Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2019 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 "ActionRegister.h"
24 : #include "tools/Torsion.h"
25 :
26 : #include <string>
27 : #include <cmath>
28 :
29 : using namespace std;
30 :
31 : namespace PLMD {
32 : namespace colvar {
33 :
34 : //+PLUMEDOC COLVAR TORSION
35 : /*
36 : Calculate a torsional angle.
37 :
38 : This command can be used to compute the torsion between four atoms or alternatively
39 : to calculate the angle between two vectors projected on the plane
40 : orthogonal to an axis.
41 :
42 : \par Examples
43 :
44 : This input tells plumed to print the torsional angle between atoms 1, 2, 3 and 4
45 : on file COLVAR.
46 : \plumedfile
47 : t: TORSION ATOMS=1,2,3,4
48 : # this is an alternative, equivalent, definition:
49 : # t: TORSION VECTOR1=2,1 AXIS=2,3 VECTOR2=3,4
50 : PRINT ARG=t FILE=COLVAR
51 : \endplumedfile
52 :
53 : 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$
54 : by using TORSION in combination with the \ref MOLINFO command. This can be done by using the following
55 : syntax.
56 :
57 : \plumedfile
58 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
59 : MOLINFO MOLTYPE=protein STRUCTURE=myprotein.pdb
60 : t1: TORSION ATOMS=@phi-3
61 : t2: TORSION ATOMS=@psi-4
62 : PRINT ARG=t1,t2 FILE=colvar STRIDE=10
63 : \endplumedfile
64 :
65 : Here, \@phi-3 tells plumed that you would like to calculate the \f$\phi\f$ angle in the third residue of the protein.
66 : Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle of the fourth residue of the protein.
67 :
68 : Both of the previous examples specify that the torsion angle should be calculated based on the position of four atoms.
69 : For the first example in particular the assumption when the torsion is specified in this way is that there are chemical
70 : bonds between atoms 1 and 2, atoms 2 and 3 and atoms 3 and 4. In general, however, a torsional angle measures the angle
71 : between two planes, which have at least one vector in common. As shown below, there is thus an alternate, more general, way
72 : through which we can define a torsional angle:
73 :
74 : \plumedfile
75 : t1: TORSION VECTOR1=1,2 AXIS=3,4 VECTOR2=5,6
76 : PRINT ARG=t1 FILE=colvar STRIDE=20
77 : \endplumedfile
78 :
79 : This input instructs PLUMED to calculate the angle between the plane containing the vector connecting atoms 1 and 2 and the vector
80 : connecting atoms 3 and 4 and the plane containing this second vector and the vector connecting atoms 5 and 6. We can even use
81 : PLUMED to calculate the torsional angle between two bond vectors around the z-axis as shown below:
82 :
83 : \plumedfile
84 : a0: FIXEDATOM AT=0,0,0
85 : az: FIXEDATOM AT=0,0,1
86 : t1: TORSION VECTOR1=1,2 AXIS=a0,az VECTOR2=5,6
87 : PRINT ARG=t1 FILE=colvar STRIDE=20
88 : \endplumedfile
89 :
90 :
91 : */
92 : //+ENDPLUMEDOC
93 :
94 1014 : class Torsion : public Colvar {
95 : bool pbc;
96 : bool do_cosine;
97 :
98 : public:
99 : explicit Torsion(const ActionOptions&);
100 : // active methods:
101 : void calculate() override;
102 : static void registerKeywords(Keywords& keys);
103 : };
104 :
105 8848 : PLUMED_REGISTER_ACTION(Torsion,"TORSION")
106 :
107 510 : void Torsion::registerKeywords(Keywords& keys) {
108 510 : Colvar::registerKeywords( keys );
109 2040 : keys.add("atoms-1","ATOMS","the four atoms involved in the torsional angle");
110 2040 : 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 VECTOR1 and VECTOR2 keywords.");
111 2040 : keys.add("atoms-2","VECTOR1","two atoms that define a vector. You can use this in combination with VECTOR2 and AXIS");
112 2040 : keys.add("atoms-2","VECTOR2","two atoms that define a vector. You can use this in combination with VECTOR1 and AXIS");
113 1530 : keys.addFlag("COSINE",false,"calculate cosine instead of dihedral");
114 510 : }
115 :
116 509 : Torsion::Torsion(const ActionOptions&ao):
117 : PLUMED_COLVAR_INIT(ao),
118 : pbc(true),
119 511 : do_cosine(false)
120 : {
121 : vector<AtomNumber> atoms,v1,v2,axis;
122 1020 : parseAtomList("ATOMS",atoms);
123 1020 : parseAtomList("VECTOR1",v1);
124 1020 : parseAtomList("VECTOR2",v2);
125 1020 : parseAtomList("AXIS",axis);
126 :
127 1020 : parseFlag("COSINE",do_cosine);
128 :
129 509 : bool nopbc=!pbc;
130 1020 : parseFlag("NOPBC",nopbc);
131 509 : pbc=!nopbc;
132 509 : checkRead();
133 :
134 509 : if(atoms.size()==4) {
135 1510 : if(!(v1.empty() && v2.empty() && axis.empty()))
136 4 : error("ATOMS keyword is not compatible with VECTOR1, VECTOR2 and AXIS keywords");
137 503 : log.printf(" between atoms %d %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
138 503 : atoms.resize(6);
139 503 : atoms[5]=atoms[3];
140 503 : atoms[4]=atoms[2];
141 503 : atoms[3]=atoms[2];
142 503 : atoms[2]=atoms[1];
143 5 : } else if(atoms.empty()) {
144 12 : if(!(v1.size()==2 && v2.size()==2 && axis.size()==2))
145 2 : error("VECTOR1, VECTOR2 and AXIS should specify 2 atoms each");
146 : log.printf(" between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n",
147 4 : v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial());
148 4 : atoms.resize(6);
149 4 : atoms[0]=v1[1];
150 4 : atoms[1]=v1[0];
151 4 : atoms[2]=axis[0];
152 4 : atoms[3]=axis[1];
153 4 : atoms[4]=v2[0];
154 4 : atoms[5]=v2[1];
155 4 : } else error("ATOMS should specify 4 atoms");
156 :
157 507 : if(pbc) log.printf(" using periodic boundary conditions\n");
158 111 : else log.printf(" without periodic boundary conditions\n");
159 :
160 507 : if(do_cosine) log.printf(" calculating cosine instead of torsion\n");
161 :
162 507 : addValueWithDerivatives();
163 1521 : if(!do_cosine) setPeriodic("-pi","pi");
164 0 : else setNotPeriodic();
165 507 : requestAtoms(atoms);
166 507 : }
167 :
168 : // calculator
169 29580 : void Torsion::calculate() {
170 :
171 29580 : Vector d0,d1,d2;
172 29580 : if(pbc) makeWhole();
173 29580 : d0=delta(getPosition(1),getPosition(0));
174 29580 : d1=delta(getPosition(3),getPosition(2));
175 29580 : d2=delta(getPosition(5),getPosition(4));
176 29580 : Vector dd0,dd1,dd2;
177 : PLMD::Torsion t;
178 29580 : double torsion=t.compute(d0,d1,d2,dd0,dd1,dd2);
179 29580 : if(do_cosine) {
180 0 : dd0 *= -sin(torsion);
181 0 : dd1 *= -sin(torsion);
182 0 : dd2 *= -sin(torsion);
183 0 : torsion = cos(torsion);
184 : }
185 29580 : setAtomsDerivatives(0,dd0);
186 29580 : setAtomsDerivatives(1,-dd0);
187 29580 : setAtomsDerivatives(2,dd1);
188 29580 : setAtomsDerivatives(3,-dd1);
189 29580 : setAtomsDerivatives(4,dd2);
190 29580 : setAtomsDerivatives(5,-dd2);
191 :
192 29580 : setValue (torsion);
193 29580 : setBoxDerivativesNoPbc();
194 29580 : }
195 :
196 : }
197 5874 : }
198 :
199 :
200 :
|
