Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "Colvar.h"
23 : #include "ColvarShortcut.h"
24 : #include "MultiColvarTemplate.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/Angle.h"
27 :
28 : namespace PLMD {
29 : namespace colvar {
30 :
31 : //+PLUMEDOC COLVAR ANGLE
32 : /*
33 : Calculate an angle.
34 :
35 : This command can be used to compute the angle between three atoms. Alternatively
36 : if four atoms appear in the atom
37 : specification it calculates the angle between
38 : two vectors identified by two pairs of atoms.
39 :
40 : If _three_ atoms are given, the angle is defined as:
41 : \f[
42 : \theta=\arccos\left(\frac{ {\bf r}_{21}\cdot {\bf r}_{23}}{
43 : |{\bf r}_{21}| |{\bf r}_{23}|}\right)
44 : \f]
45 : Here \f$ {\bf r}_{ij}\f$ is the distance vector among the
46 : \f$i\f$th and the \f$j\f$th listed atom.
47 :
48 : If _four_ atoms are given, the angle is defined as:
49 : \f[
50 : \theta=\arccos\left(\frac{ {\bf r}_{21}\cdot {\bf r}_{34}}{
51 : |{\bf r}_{21}| |{\bf r}_{34}|}\right)
52 : \f]
53 :
54 : Notice that angles defined in this way are non-periodic variables and
55 : their value is limited by definition between 0 and \f$\pi\f$.
56 :
57 : The vectors \f$ {\bf r}_{ij}\f$ are by default evaluated taking
58 : periodic boundary conditions into account.
59 : This behavior can be changed with the NOPBC flag.
60 :
61 : \par Examples
62 :
63 : This command tells plumed to calculate the angle between the vector connecting atom 1 to atom 2 and
64 : the vector connecting atom 2 to atom 3 and to print it on file COLVAR1. At the same time,
65 : the angle between vector connecting atom 1 to atom 2 and the vector connecting atom 3 to atom 4 is printed
66 : on file COLVAR2.
67 : \plumedfile
68 :
69 : a: ANGLE ATOMS=1,2,3
70 : # equivalently one could state:
71 : # a: ANGLE ATOMS=1,2,2,3
72 :
73 : b: ANGLE ATOMS=1,2,3,4
74 :
75 : PRINT ARG=a FILE=COLVAR1
76 : PRINT ARG=b FILE=COLVAR2
77 : \endplumedfile
78 :
79 :
80 : */
81 : //+ENDPLUMEDOC
82 :
83 : //+PLUMEDOC COLVAR ANGLE_SCALAR
84 : /*
85 : Calculate an angle.
86 :
87 : \par Examples
88 :
89 : */
90 : //+ENDPLUMEDOC
91 :
92 : //+PLUMEDOC COLVAR ANGLE_VECTOR
93 : /*
94 : Calculate multiple angles.
95 :
96 : \par Examples
97 :
98 : */
99 : //+ENDPLUMEDOC
100 :
101 : class Angle : public Colvar {
102 : bool pbc;
103 : std::vector<double> value, masses, charges;
104 : std::vector<std::vector<Vector> > derivs;
105 : std::vector<Tensor> virial;
106 : public:
107 : explicit Angle(const ActionOptions&);
108 : // active methods:
109 : void calculate() override;
110 : static void registerKeywords( Keywords& keys );
111 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
112 : static unsigned getModeAndSetupValues( ActionWithValue* av );
113 : static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
114 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
115 : std::vector<Tensor>& virial, const ActionAtomistic* aa );
116 : };
117 :
118 : typedef ColvarShortcut<Angle> AngleShortcut;
119 : PLUMED_REGISTER_ACTION(AngleShortcut,"ANGLE")
120 : PLUMED_REGISTER_ACTION(Angle,"ANGLE_SCALAR")
121 : typedef MultiColvarTemplate<Angle> AngleMulti;
122 : PLUMED_REGISTER_ACTION(AngleMulti,"ANGLE_VECTOR")
123 :
124 89 : void Angle::registerKeywords( Keywords& keys ) {
125 89 : Colvar::registerKeywords(keys);
126 89 : keys.setDisplayName("ANGLE");
127 178 : keys.add("atoms","ATOMS","the list of atoms involved in this collective variable (either 3 or 4 atoms)");
128 178 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
129 89 : keys.setValueDescription("the ANGLE involving these atoms");
130 89 : }
131 :
132 73 : void Angle::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
133 146 : aa->parseAtomList("ATOMS",num,atoms);
134 73 : if(atoms.size()==3) {
135 62 : aa->log.printf(" between atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial());
136 62 : atoms.resize(4);
137 62 : atoms[3]=atoms[2];
138 62 : atoms[2]=atoms[1];
139 11 : } else if(atoms.size()==4) {
140 7 : aa->log.printf(" between lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
141 4 : } else if( num<0 || atoms.size()>0 ) {
142 1 : aa->error("Number of specified atoms should be either 3 or 4");
143 : }
144 72 : }
145 :
146 3 : unsigned Angle::getModeAndSetupValues( ActionWithValue* av ) {
147 3 : av->addValueWithDerivatives();
148 3 : av->setNotPeriodic();
149 3 : return 0;
150 : }
151 :
152 23 : Angle::Angle(const ActionOptions&ao):
153 : PLUMED_COLVAR_INIT(ao),
154 23 : pbc(true),
155 23 : value(1),
156 24 : derivs(1),
157 46 : virial(1) {
158 23 : derivs[0].resize(4);
159 : std::vector<AtomNumber> atoms;
160 23 : parseAtomList( -1, atoms, this );
161 22 : bool nopbc=!pbc;
162 22 : parseFlag("NOPBC",nopbc);
163 22 : pbc=!nopbc;
164 :
165 22 : if(pbc) {
166 22 : log.printf(" using periodic boundary conditions\n");
167 : } else {
168 0 : log.printf(" without periodic boundary conditions\n");
169 : }
170 :
171 23 : addValueWithDerivatives();
172 22 : setNotPeriodic();
173 22 : requestAtoms(atoms);
174 22 : checkRead();
175 25 : }
176 :
177 : // calculator
178 319 : void Angle::calculate() {
179 :
180 319 : if(pbc) {
181 319 : makeWhole();
182 : }
183 319 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
184 319 : setValue( value[0] );
185 1595 : for(unsigned i=0; i<derivs[0].size(); ++i) {
186 1276 : setAtomsDerivatives( i, derivs[0][i] );
187 : }
188 319 : setBoxDerivatives( virial[0] );
189 319 : }
190 :
191 554 : void Angle::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
192 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
193 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
194 554 : Vector dij,dik;
195 554 : dij=delta(pos[2],pos[3]);
196 554 : dik=delta(pos[1],pos[0]);
197 554 : Vector ddij,ddik;
198 : PLMD::Angle a;
199 554 : vals[0]=a.compute(dij,dik,ddij,ddik);
200 554 : derivs[0][0]=ddik;
201 554 : derivs[0][1]=-ddik;
202 554 : derivs[0][2]=-ddij;
203 554 : derivs[0][3]=ddij;
204 554 : setBoxDerivativesNoPbc( pos, derivs, virial );
205 554 : }
206 :
207 : }
208 : }
209 :
210 :
211 :
|
