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/Torsion.h"
27 :
28 : namespace PLMD {
29 : namespace colvar {
30 :
31 : //+PLUMEDOC COLVAR TORSION
32 : /*
33 : Calculate one or multiple torsional angles.
34 :
35 : This command can be used to compute the torsion between four atoms as shown by the input below:
36 :
37 : ```plumed
38 : t: TORSION ATOMS=1,2,3,4
39 : PRINT ARG=t FILE=COLVAR
40 : ```
41 :
42 : Alternatively you can use this action to calculate the angle between two vectors projected on the plane
43 : orthogonal to an axis. The example below uses this syntax and computes the same torsion as the first example
44 : input above.
45 :
46 : ```plumed
47 : t: TORSION VECTOR1=2,1 AXIS=2,3 VECTOR2=3,4
48 : PRINT ARG=t FILE=COLVAR
49 : ```
50 :
51 : If you combine this sytax with the functionality in [FIXEDATOM](FIXEDATOM.md) you can see how we can calculate the
52 : torsional angle between two bond vectors around the z-axis as shown below:
53 :
54 : ```plumed
55 : a0: FIXEDATOM AT=0,0,0
56 : az: FIXEDATOM AT=0,0,1
57 : t1: TORSION VECTOR1=1,2 AXIS=a0,az VECTOR2=5,6
58 : PRINT ARG=t1 FILE=colvar STRIDE=20
59 : ```
60 :
61 : If you are working with a protein you can specify the special named torsion angles $\phi$, $\psi$, $\omega$ and $\chi_1$
62 : by using TORSION in combination with the [MOLINFO](MOLINFO.md) command. This can be done by using the following
63 : syntax.
64 :
65 : ```plumed
66 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
67 : MOLINFO MOLTYPE=protein STRUCTURE=regtest/basic/rt32/helix.pdb
68 : t1: TORSION ATOMS=@phi-3
69 : t2: TORSION ATOMS=@psi-4
70 : PRINT ARG=t1,t2 FILE=colvar STRIDE=10
71 : ```
72 :
73 : Here, `@phi-3` tells plumed that you would like to calculate the $\phi$ angle in the third residue of the protein.
74 : Similarly `@psi-4` tells plumed that you want to calculate the $\psi$ angle of the fourth residue of the protein.
75 :
76 : If you would like to calculate multiple torsion angles at the same time you can use a command like the one shown below:
77 :
78 : ```plumed
79 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
80 : MOLINFO MOLTYPE=protein STRUCTURE=regtest/basic/rt32/helix.pdb
81 : t1: TORSION ATOMS1=@phi-3 ATOMS2=@phi-4 ATOMS3=@phi-5 ATOMS4=@phi-6 ATOMS5=@phi-7
82 : PRINT ARG=t1 FILE=colvar STRIDE=10
83 : ```
84 :
85 : This input tells PLUMED to calculate the $\phi$ angles in residues 3-7 of the protein. The output, `t1`, is a 5 dimensional vector.
86 :
87 :
88 : */
89 : //+ENDPLUMEDOC
90 :
91 : //+PLUMEDOC COLVAR TORSION_SCALAR
92 : /*
93 : Calculate a torsional angle.
94 :
95 : \par Examples
96 :
97 : */
98 : //+ENDPLUMEDOC
99 :
100 : //+PLUMEDOC MCOLVAR TORSION_VECTOR
101 : /*
102 : Calculate multiple torsional angles.
103 :
104 : \par Examples
105 :
106 : */
107 : //+ENDPLUMEDOC
108 :
109 : class Torsion : public Colvar {
110 : bool pbc;
111 : bool do_cosine;
112 :
113 : std::vector<double> value, masses, charges;
114 : std::vector<std::vector<Vector> > derivs;
115 : std::vector<Tensor> virial;
116 : public:
117 : explicit Torsion(const ActionOptions&);
118 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
119 : static unsigned getModeAndSetupValues( ActionWithValue* av );
120 : // active methods:
121 : void calculate() override;
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 : static void registerKeywords(Keywords& keys);
126 : };
127 :
128 : typedef ColvarShortcut<Torsion> TorsionShortcut;
129 : PLUMED_REGISTER_ACTION(TorsionShortcut,"TORSION")
130 : PLUMED_REGISTER_ACTION(Torsion,"TORSION_SCALAR")
131 : typedef MultiColvarTemplate<Torsion> TorsionMulti;
132 : PLUMED_REGISTER_ACTION(TorsionMulti,"TORSION_VECTOR")
133 :
134 2130 : void Torsion::registerKeywords(Keywords& keys) {
135 2130 : Colvar::registerKeywords( keys );
136 4260 : keys.setDisplayName("TORSION");
137 2130 : keys.add("atoms-1","ATOMS","the four atoms involved in the torsional angle");
138 2130 : keys.add("atoms-2","AXIS","two atoms that define an axis. You can use this to find the angle in the plane perpendicular to the axis between the vectors specified using the VECTORA and VECTORB keywords.");
139 2130 : keys.add("atoms-2","VECTORA","two atoms that define a vector. You can use this in combination with VECTOR2 and AXIS");
140 2130 : keys.add("atoms-2","VECTORB","two atoms that define a vector. You can use this in combination with VECTOR1 and AXIS");
141 2130 : keys.add("atoms-3","VECTOR1","two atoms that define a vector. You can use this in combination with VECTOR2 and AXIS");
142 2130 : keys.add("atoms-3","VECTOR2","two atoms that define a vector. You can use this in combination with VECTOR1 and AXIS");
143 2130 : keys.addFlag("COSINE",false,"calculate cosine instead of dihedral");
144 2130 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
145 4260 : keys.setValueDescription("scalar/vector","the TORSION involving these atoms");
146 2130 : }
147 :
148 691 : Torsion::Torsion(const ActionOptions&ao):
149 : PLUMED_COLVAR_INIT(ao),
150 691 : pbc(true),
151 691 : do_cosine(false),
152 691 : value(1),
153 693 : derivs(1),
154 1382 : virial(1) {
155 691 : derivs[0].resize(6);
156 : std::vector<AtomNumber> atoms;
157 : std::vector<AtomNumber> v1;
158 1382 : ActionAtomistic::parseAtomList("VECTOR1",v1);
159 691 : if( v1.size()>0 ) {
160 : std::vector<AtomNumber> v2;
161 4 : ActionAtomistic::parseAtomList("VECTOR2",v2);
162 : std::vector<AtomNumber> axis;
163 4 : ActionAtomistic::parseAtomList("AXIS",axis);
164 2 : if( !(v1.size()==2 && v2.size()==2 && axis.size()==2)) {
165 0 : error("VECTOR1, VECTOR2 and AXIS should specify 2 atoms each");
166 : }
167 2 : atoms.resize(6);
168 2 : atoms[0]=v1[1];
169 2 : atoms[1]=v1[0];
170 2 : atoms[2]=axis[0];
171 2 : atoms[3]=axis[1];
172 2 : atoms[4]=v2[0];
173 2 : atoms[5]=v2[1];
174 2 : log.printf(" between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n",
175 : v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial());
176 : } else {
177 689 : parseAtomList(-1,atoms,this);
178 : }
179 690 : unsigned mode=getModeAndSetupValues(this);
180 690 : if( mode==1 ) {
181 3 : do_cosine=true;
182 : }
183 :
184 690 : bool nopbc=!pbc;
185 692 : parseFlag("NOPBC",nopbc);
186 690 : pbc=!nopbc;
187 690 : checkRead();
188 :
189 689 : if(pbc) {
190 577 : log.printf(" using periodic boundary conditions\n");
191 : } else {
192 112 : log.printf(" without periodic boundary conditions\n");
193 : }
194 689 : requestAtoms(atoms);
195 695 : }
196 :
197 794 : void Torsion::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
198 : std::vector<AtomNumber> v1,v2,axis;
199 794 : aa->parseAtomList("ATOMS",num,t);
200 794 : aa->parseAtomList("VECTORA",num,v1);
201 794 : aa->parseAtomList("VECTORB",num,v2);
202 1588 : aa->parseAtomList("AXIS",num,axis);
203 :
204 794 : if(t.size()==4) {
205 747 : if(!(v1.empty() && v2.empty() && axis.empty())) {
206 0 : aa->error("ATOMS keyword is not compatible with VECTORA, VECTORB and AXIS keywords");
207 : }
208 747 : aa->log.printf(" between atoms %d %d %d %d\n",t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial());
209 747 : t.resize(6);
210 747 : t[5]=t[3];
211 747 : t[4]=t[2];
212 747 : t[3]=t[2];
213 747 : t[2]=t[1];
214 47 : } else if(t.empty()) {
215 46 : if( num>0 && v1.empty() && v2.empty() && axis.empty() ) {
216 : return;
217 : }
218 32 : if(!(v1.size()==2 && v2.size()==2 && axis.size()==2)) {
219 0 : aa->error("VECTORA, VECTORB and AXIS should specify 2 atoms each");
220 : }
221 32 : aa->log.printf(" between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n",
222 : v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial());
223 32 : t.resize(6);
224 32 : t[0]=v1[1];
225 32 : t[1]=v1[0];
226 32 : t[2]=axis[0];
227 32 : t[3]=axis[1];
228 32 : t[4]=v2[0];
229 32 : t[5]=v2[1];
230 : } else if( t.size()!=4 ) {
231 2 : aa->error("ATOMS should specify 4 atoms");
232 : }
233 : }
234 :
235 704 : unsigned Torsion::getModeAndSetupValues( ActionWithValue* av ) {
236 : bool do_cos;
237 704 : av->parseFlag("COSINE",do_cos);
238 704 : if(do_cos) {
239 4 : av->log.printf(" calculating cosine instead of torsion\n");
240 : }
241 :
242 704 : av->addValueWithDerivatives();
243 704 : if(!do_cos) {
244 1400 : av->setPeriodic("-pi","pi");
245 700 : return 0;
246 : }
247 4 : av->setNotPeriodic();
248 : return 1;
249 : }
250 :
251 : // calculator
252 37137 : void Torsion::calculate() {
253 37137 : if(pbc) {
254 34865 : makeWhole();
255 : }
256 37137 : if(do_cosine) {
257 15 : calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this );
258 : } else {
259 37122 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
260 : }
261 259959 : for(unsigned i=0; i<6; ++i) {
262 222822 : setAtomsDerivatives(i,derivs[0][i] );
263 : }
264 37137 : setValue(value[0]);
265 37137 : setBoxDerivatives( virial[0] );
266 37137 : }
267 :
268 37916 : void Torsion::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
269 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
270 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
271 37916 : Vector d0=delta(pos[1],pos[0]);
272 37916 : Vector d1=delta(pos[3],pos[2]);
273 37916 : Vector d2=delta(pos[5],pos[4]);
274 37916 : Vector dd0,dd1,dd2;
275 : PLMD::Torsion t;
276 37916 : vals[0] = t.compute(d0,d1,d2,dd0,dd1,dd2);
277 37916 : if(mode==1) {
278 30 : dd0 *= -std::sin(vals[0]);
279 30 : dd1 *= -std::sin(vals[0]);
280 30 : dd2 *= -std::sin(vals[0]);
281 30 : vals[0] = std::cos(vals[0]);
282 : }
283 37916 : derivs[0][0] = dd0;
284 37916 : derivs[0][1] = -dd0;
285 37916 : derivs[0][2] = dd1;
286 37916 : derivs[0][3] = -dd1;
287 37916 : derivs[0][4] = dd2;
288 37916 : derivs[0][5] = -dd2;
289 37916 : setBoxDerivativesNoPbc( pos, derivs, virial );
290 37916 : }
291 :
292 : }
293 : }
294 :
295 :
296 :
|
