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 a torsional angle.
34 :
35 : This command can be used to compute the torsion between four atoms or alternatively
36 : to calculate the angle between two vectors projected on the plane
37 : orthogonal to an axis.
38 :
39 : \par Examples
40 :
41 : This input tells plumed to print the torsional angle between atoms 1, 2, 3 and 4
42 : on file COLVAR.
43 : \plumedfile
44 : t: TORSION ATOMS=1,2,3,4
45 : # this is an alternative, equivalent, definition:
46 : # t: TORSION VECTOR1=2,1 AXIS=2,3 VECTOR2=3,4
47 : PRINT ARG=t FILE=COLVAR
48 : \endplumedfile
49 :
50 : If you are working with a protein you can specify the special named torsion angles \f$\phi\f$, \f$\psi\f$, \f$\omega\f$ and \f$\chi_1\f$
51 : by using TORSION in combination with the \ref MOLINFO command. This can be done by using the following
52 : syntax.
53 :
54 : \plumedfile
55 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
56 : MOLINFO MOLTYPE=protein STRUCTURE=myprotein.pdb
57 : t1: TORSION ATOMS=@phi-3
58 : t2: TORSION ATOMS=@psi-4
59 : PRINT ARG=t1,t2 FILE=colvar STRIDE=10
60 : \endplumedfile
61 :
62 : Here, \@phi-3 tells plumed that you would like to calculate the \f$\phi\f$ angle in the third residue of the protein.
63 : Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle of the fourth residue of the protein.
64 :
65 : Both of the previous examples specify that the torsion angle should be calculated based on the position of four atoms.
66 : For the first example in particular the assumption when the torsion is specified in this way is that there are chemical
67 : bonds between atoms 1 and 2, atoms 2 and 3 and atoms 3 and 4. In general, however, a torsional angle measures the angle
68 : between two planes, which have at least one vector in common. As shown below, there is thus an alternate, more general, way
69 : through which we can define a torsional angle:
70 :
71 : \plumedfile
72 : t1: TORSION VECTOR1=1,2 AXIS=3,4 VECTOR2=5,6
73 : PRINT ARG=t1 FILE=colvar STRIDE=20
74 : \endplumedfile
75 :
76 : This input instructs PLUMED to calculate the angle between the plane containing the vector connecting atoms 1 and 2 and the vector
77 : connecting atoms 3 and 4 and the plane containing this second vector and the vector connecting atoms 5 and 6. We can even use
78 : PLUMED to calculate the torsional angle between two bond vectors around the z-axis as shown below:
79 :
80 : \plumedfile
81 : a0: FIXEDATOM AT=0,0,0
82 : az: FIXEDATOM AT=0,0,1
83 : t1: TORSION VECTOR1=1,2 AXIS=a0,az VECTOR2=5,6
84 : PRINT ARG=t1 FILE=colvar STRIDE=20
85 : \endplumedfile
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 2187 : void Torsion::registerKeywords(Keywords& keys) {
135 2187 : Colvar::registerKeywords( keys );
136 2187 : keys.setDisplayName("TORSION");
137 4374 : keys.add("atoms-1","ATOMS","the four atoms involved in the torsional angle");
138 4374 : 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 4374 : keys.add("atoms-2","VECTORA","two atoms that define a vector. You can use this in combination with VECTOR2 and AXIS");
140 4374 : keys.add("atoms-2","VECTORB","two atoms that define a vector. You can use this in combination with VECTOR1 and AXIS");
141 4374 : keys.add("atoms-3","VECTOR1","two atoms that define a vector. You can use this in combination with VECTOR2 and AXIS");
142 4374 : keys.add("atoms-3","VECTOR2","two atoms that define a vector. You can use this in combination with VECTOR1 and AXIS");
143 4374 : keys.addFlag("COSINE",false,"calculate cosine instead of dihedral");
144 4374 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
145 2187 : keys.setValueDescription("the TORSION involving these atoms");
146 2187 : }
147 :
148 707 : Torsion::Torsion(const ActionOptions&ao):
149 : PLUMED_COLVAR_INIT(ao),
150 707 : pbc(true),
151 707 : do_cosine(false),
152 707 : value(1),
153 709 : derivs(1),
154 1414 : virial(1) {
155 707 : derivs[0].resize(6);
156 : std::vector<AtomNumber> atoms;
157 : std::vector<AtomNumber> v1;
158 1414 : ActionAtomistic::parseAtomList("VECTOR1",v1);
159 707 : 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 705 : parseAtomList(-1,atoms,this);
178 : }
179 706 : unsigned mode=getModeAndSetupValues(this);
180 706 : if( mode==1 ) {
181 3 : do_cosine=true;
182 : }
183 :
184 706 : bool nopbc=!pbc;
185 708 : parseFlag("NOPBC",nopbc);
186 706 : pbc=!nopbc;
187 706 : checkRead();
188 :
189 705 : if(pbc) {
190 589 : log.printf(" using periodic boundary conditions\n");
191 : } else {
192 116 : log.printf(" without periodic boundary conditions\n");
193 : }
194 705 : requestAtoms(atoms);
195 711 : }
196 :
197 813 : void Torsion::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
198 : std::vector<AtomNumber> v1,v2,axis;
199 813 : aa->parseAtomList("ATOMS",num,t);
200 813 : aa->parseAtomList("VECTORA",num,v1);
201 813 : aa->parseAtomList("VECTORB",num,v2);
202 1626 : aa->parseAtomList("AXIS",num,axis);
203 :
204 813 : if(t.size()==4) {
205 765 : if(!(v1.empty() && v2.empty() && axis.empty())) {
206 0 : aa->error("ATOMS keyword is not compatible with VECTORA, VECTORB and AXIS keywords");
207 : }
208 765 : aa->log.printf(" between atoms %d %d %d %d\n",t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial());
209 765 : t.resize(6);
210 765 : t[5]=t[3];
211 765 : t[4]=t[2];
212 765 : t[3]=t[2];
213 765 : t[2]=t[1];
214 48 : } else if(t.empty()) {
215 47 : 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 721 : unsigned Torsion::getModeAndSetupValues( ActionWithValue* av ) {
236 : bool do_cos;
237 721 : av->parseFlag("COSINE",do_cos);
238 721 : if(do_cos) {
239 4 : av->log.printf(" calculating cosine instead of torsion\n");
240 : }
241 :
242 721 : av->addValueWithDerivatives();
243 721 : if(!do_cos) {
244 1434 : av->setPeriodic("-pi","pi");
245 717 : return 0;
246 : }
247 4 : av->setNotPeriodic();
248 : return 1;
249 : }
250 :
251 : // calculator
252 37738 : void Torsion::calculate() {
253 37738 : if(pbc) {
254 35338 : makeWhole();
255 : }
256 37738 : if(do_cosine) {
257 15 : calculateCV( 1, masses, charges, getPositions(), value, derivs, virial, this );
258 : } else {
259 37723 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
260 : }
261 264166 : for(unsigned i=0; i<6; ++i) {
262 226428 : setAtomsDerivatives(i,derivs[0][i] );
263 : }
264 37738 : setValue(value[0]);
265 37738 : setBoxDerivatives( virial[0] );
266 37738 : }
267 :
268 38581 : 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 38581 : Vector d0=delta(pos[1],pos[0]);
272 38581 : Vector d1=delta(pos[3],pos[2]);
273 38581 : Vector d2=delta(pos[5],pos[4]);
274 38581 : Vector dd0,dd1,dd2;
275 : PLMD::Torsion t;
276 38581 : vals[0] = t.compute(d0,d1,d2,dd0,dd1,dd2);
277 38581 : 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 38581 : derivs[0][0] = dd0;
284 38581 : derivs[0][1] = -dd0;
285 38581 : derivs[0][2] = dd1;
286 38581 : derivs[0][3] = -dd1;
287 38581 : derivs[0][4] = dd2;
288 38581 : derivs[0][5] = -dd2;
289 38581 : setBoxDerivativesNoPbc( pos, derivs, virial );
290 38581 : }
291 :
292 : }
293 : }
294 :
295 :
296 :
|
