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