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 11721 : PLUMED_REGISTER_ACTION(Torsion,"TORSION")
101 :
102 653 : void Torsion::registerKeywords(Keywords& keys) {
103 653 : Colvar::registerKeywords( keys );
104 1306 : keys.add("atoms-1","ATOMS","the four atoms involved in the torsional angle");
105 1306 : 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 1306 : keys.add("atoms-2","VECTOR1","two atoms that define a vector. You can use this in combination with VECTOR2 and AXIS");
107 1306 : keys.add("atoms-2","VECTOR2","two atoms that define a vector. You can use this in combination with VECTOR1 and AXIS");
108 1306 : keys.addFlag("COSINE",false,"calculate cosine instead of dihedral");
109 653 : }
110 :
111 652 : Torsion::Torsion(const ActionOptions&ao):
112 : PLUMED_COLVAR_INIT(ao),
113 652 : pbc(true),
114 652 : do_cosine(false)
115 : {
116 : std::vector<AtomNumber> atoms,v1,v2,axis;
117 652 : parseAtomList("ATOMS",atoms);
118 652 : parseAtomList("VECTOR1",v1);
119 652 : parseAtomList("VECTOR2",v2);
120 652 : parseAtomList("AXIS",axis);
121 :
122 652 : parseFlag("COSINE",do_cosine);
123 :
124 652 : bool nopbc=!pbc;
125 652 : parseFlag("NOPBC",nopbc);
126 652 : pbc=!nopbc;
127 652 : checkRead();
128 :
129 652 : if(atoms.size()==4) {
130 647 : if(!(v1.empty() && v2.empty() && axis.empty()))
131 1 : error("ATOMS keyword is not compatible with VECTOR1, VECTOR2 and AXIS keywords");
132 646 : log.printf(" between atoms %d %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
133 646 : atoms.resize(6);
134 646 : atoms[5]=atoms[3];
135 646 : atoms[4]=atoms[2];
136 646 : atoms[3]=atoms[2];
137 646 : 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 4 : log.printf(" between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n",
142 : v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial());
143 4 : atoms.resize(6);
144 4 : atoms[0]=v1[1];
145 4 : atoms[1]=v1[0];
146 4 : atoms[2]=axis[0];
147 4 : atoms[3]=axis[1];
148 4 : atoms[4]=v2[0];
149 4 : atoms[5]=v2[1];
150 1 : } else error("ATOMS should specify 4 atoms");
151 :
152 650 : if(pbc) log.printf(" using periodic boundary conditions\n");
153 112 : else log.printf(" without periodic boundary conditions\n");
154 :
155 650 : if(do_cosine) log.printf(" calculating cosine instead of torsion\n");
156 :
157 650 : addValueWithDerivatives();
158 1302 : if(!do_cosine) setPeriodic("-pi","pi");
159 0 : else setNotPeriodic();
160 650 : requestAtoms(atoms);
161 654 : }
162 :
163 : // calculator
164 34536 : void Torsion::calculate() {
165 :
166 34536 : Vector d0,d1,d2;
167 34536 : if(pbc) makeWhole();
168 34536 : d0=delta(getPosition(1),getPosition(0));
169 34536 : d1=delta(getPosition(3),getPosition(2));
170 34536 : d2=delta(getPosition(5),getPosition(4));
171 34536 : Vector dd0,dd1,dd2;
172 : PLMD::Torsion t;
173 34536 : double torsion=t.compute(d0,d1,d2,dd0,dd1,dd2);
174 34536 : if(do_cosine) {
175 0 : dd0 *= -std::sin(torsion);
176 0 : dd1 *= -std::sin(torsion);
177 0 : dd2 *= -std::sin(torsion);
178 0 : torsion = std::cos(torsion);
179 : }
180 34536 : setAtomsDerivatives(0,dd0);
181 34536 : setAtomsDerivatives(1,-dd0);
182 34536 : setAtomsDerivatives(2,dd1);
183 34536 : setAtomsDerivatives(3,-dd1);
184 34536 : setAtomsDerivatives(4,dd2);
185 34536 : setAtomsDerivatives(5,-dd2);
186 :
187 34536 : setValue (torsion);
188 34536 : setBoxDerivativesNoPbc();
189 34536 : }
190 :
191 : }
192 : }
193 :
194 :
195 :
|
