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 one or multiple angle/s.
34 :
35 : The following input instructs PLUMED to calculate and print the angle between the vector
36 : connecting atom 2 and atom 1 and the vector connecting atom 2 and atom 3.
37 :
38 : ```plumed
39 : a1: ANGLE ATOMS=1,2,3
40 : PRINT ARG=a1 FILE=colvar
41 : ```
42 :
43 : In other words, the angle that is output by the input above is calculated as:
44 :
45 : $$
46 : \theta=\arccos\left(\frac{ r_{21}\cdot r_{23}}{
47 : | r_{21}| | r_{23}|}\right)
48 : $$
49 :
50 : Here $r_{ij}$ is the vector connecting the $i$th and $j$th atoms, which by default is evaluated
51 : in a way that takes periodic boundary conditions into account. If you wish to disregard the PBC you
52 : can use the NOPBC flag.
53 :
54 : Alternatively, we can instruct PLUMED to calculate the angle between the vectors connecting
55 : atoms 1 and atom 2 and atoms 3 and atom 4 by using the following input:
56 :
57 : ```plumed
58 : a2: ANGLE ATOMS=1,2,3,4
59 : PRINT ARG=a2 FILE=colvar
60 :
61 : The angle in this input is calculated using:
62 :
63 : $$
64 : \theta=\arccos\left(\frac{ r_{21}\cdot r_{34}}{
65 : | r_{21}| | r_{34}|}\right)
66 : $$
67 :
68 : Notice that angles defined in this way are non-periodic variables - their values must lie in between 0 and $\pi$.
69 :
70 : You can specify multiple sets of three or four atoms to calculate vectors of angles as illustrated in the following
71 : input which instructs PLUMED to calculate and output three angles:
72 :
73 : ```plumed
74 : a3: ANGLE ATOMS1=1,2,3 ATOMS2=4,5,6 ATOMS3=7,8,9
75 : PRINT ARG=a3 FILE=colvar
76 : ```
77 :
78 : It is common to assume when using this feature that all the angles being computed are indistinguishable
79 : so it makes sense to perform the same series of operations on every element of the output vector. The input
80 : file below is more approrpriate if the angles are not indistinguishable:
81 :
82 : ```plumed
83 : a4: ANGLE ATOMS=1,2,3
84 : a5: ANGLE ATOMS=4,5,6
85 : a6: ANGLE ATOMS=7,8,9
86 : PRINT ARG=a4,a5,a6 FILE=colvar
87 : ```
88 :
89 : */
90 : //+ENDPLUMEDOC
91 :
92 : //+PLUMEDOC COLVAR ANGLE_SCALAR
93 : /*
94 : Calculate an angle.
95 :
96 : \par Examples
97 :
98 : */
99 : //+ENDPLUMEDOC
100 :
101 : //+PLUMEDOC COLVAR ANGLE_VECTOR
102 : /*
103 : Calculate multiple angles.
104 :
105 : \par Examples
106 :
107 : */
108 : //+ENDPLUMEDOC
109 :
110 : class Angle : public Colvar {
111 : bool pbc;
112 : std::vector<double> value, masses, charges;
113 : std::vector<std::vector<Vector> > derivs;
114 : std::vector<Tensor> virial;
115 : public:
116 : explicit Angle(const ActionOptions&);
117 : // active methods:
118 : void calculate() override;
119 : static void registerKeywords( Keywords& keys );
120 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
121 : static unsigned getModeAndSetupValues( ActionWithValue* av );
122 : static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
123 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
124 : std::vector<Tensor>& virial, const ActionAtomistic* aa );
125 : };
126 :
127 : typedef ColvarShortcut<Angle> AngleShortcut;
128 : PLUMED_REGISTER_ACTION(AngleShortcut,"ANGLE")
129 : PLUMED_REGISTER_ACTION(Angle,"ANGLE_SCALAR")
130 : typedef MultiColvarTemplate<Angle> AngleMulti;
131 : PLUMED_REGISTER_ACTION(AngleMulti,"ANGLE_VECTOR")
132 :
133 89 : void Angle::registerKeywords( Keywords& keys ) {
134 89 : Colvar::registerKeywords(keys);
135 178 : keys.setDisplayName("ANGLE");
136 89 : keys.add("atoms","ATOMS","the list of atoms involved in this collective variable (either 3 or 4 atoms)");
137 89 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
138 178 : keys.setValueDescription("scalar/vector","the ANGLE involving these atoms");
139 89 : }
140 :
141 73 : void Angle::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
142 146 : aa->parseAtomList("ATOMS",num,atoms);
143 73 : if(atoms.size()==3) {
144 62 : aa->log.printf(" between atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial());
145 62 : atoms.resize(4);
146 62 : atoms[3]=atoms[2];
147 62 : atoms[2]=atoms[1];
148 11 : } else if(atoms.size()==4) {
149 7 : aa->log.printf(" between lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
150 4 : } else if( num<0 || atoms.size()>0 ) {
151 1 : aa->error("Number of specified atoms should be either 3 or 4");
152 : }
153 72 : }
154 :
155 3 : unsigned Angle::getModeAndSetupValues( ActionWithValue* av ) {
156 3 : av->addValueWithDerivatives();
157 3 : av->setNotPeriodic();
158 3 : return 0;
159 : }
160 :
161 23 : Angle::Angle(const ActionOptions&ao):
162 : PLUMED_COLVAR_INIT(ao),
163 23 : pbc(true),
164 23 : value(1),
165 24 : derivs(1),
166 46 : virial(1) {
167 23 : derivs[0].resize(4);
168 : std::vector<AtomNumber> atoms;
169 23 : parseAtomList( -1, atoms, this );
170 22 : bool nopbc=!pbc;
171 22 : parseFlag("NOPBC",nopbc);
172 22 : pbc=!nopbc;
173 :
174 22 : if(pbc) {
175 22 : log.printf(" using periodic boundary conditions\n");
176 : } else {
177 0 : log.printf(" without periodic boundary conditions\n");
178 : }
179 :
180 23 : addValueWithDerivatives();
181 22 : setNotPeriodic();
182 22 : requestAtoms(atoms);
183 22 : checkRead();
184 25 : }
185 :
186 : // calculator
187 319 : void Angle::calculate() {
188 :
189 319 : if(pbc) {
190 319 : makeWhole();
191 : }
192 319 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
193 319 : setValue( value[0] );
194 1595 : for(unsigned i=0; i<derivs[0].size(); ++i) {
195 1276 : setAtomsDerivatives( i, derivs[0][i] );
196 : }
197 319 : setBoxDerivatives( virial[0] );
198 319 : }
199 :
200 554 : void Angle::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
201 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
202 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
203 554 : Vector dij,dik;
204 554 : dij=delta(pos[2],pos[3]);
205 554 : dik=delta(pos[1],pos[0]);
206 554 : Vector ddij,ddik;
207 : PLMD::Angle a;
208 554 : vals[0]=a.compute(dij,dik,ddij,ddik);
209 554 : derivs[0][0]=ddik;
210 554 : derivs[0][1]=-ddik;
211 554 : derivs[0][2]=-ddij;
212 554 : derivs[0][3]=ddij;
213 554 : setBoxDerivativesNoPbc( pos, derivs, virial );
214 554 : }
215 :
216 : }
217 : }
218 :
219 :
220 :
|
