Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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 : \verbatim
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 : \endverbatim
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 : \verbatim
58 : MOLINFO MOLTYPE=protein STRUCTURE=myprotein.pdb
59 : t1: TORSION ATOMS=@phi-3
60 : t2: TORSION ATOMS=@psi-4
61 : PRINT ARG=t1,t2 FILE=colvar STRIDE=10
62 : \endverbatim
63 :
64 : Here, \@phi-3 tells plumed that you would like to calculate the \f$\phi\f$ angle in the third residue of the protein.
65 : Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle of the 4th residue of the protein.
66 : */
67 : //+ENDPLUMEDOC
68 :
69 142 : class Torsion : public Colvar {
70 : bool pbc;
71 : bool do_cosine;
72 :
73 : public:
74 : explicit Torsion(const ActionOptions&);
75 : // active methods:
76 : virtual void calculate();
77 : static void registerKeywords(Keywords& keys);
78 : };
79 :
80 2594 : PLUMED_REGISTER_ACTION(Torsion,"TORSION")
81 :
82 72 : void Torsion::registerKeywords(Keywords& keys) {
83 72 : Colvar::registerKeywords( keys );
84 72 : keys.add("atoms-1","ATOMS","the four atoms involved in the torsional angle");
85 72 : 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.");
86 72 : keys.add("atoms-2","VECTOR1","two atoms that define a vector. You can use this in combination with VECTOR2 and AXIS");
87 72 : keys.add("atoms-2","VECTOR2","two atoms that define a vector. You can use this in combination with VECTOR1 and AXIS");
88 72 : keys.addFlag("COSINE",false,"calculate cosine instead of dihedral");
89 72 : }
90 :
91 71 : Torsion::Torsion(const ActionOptions&ao):
92 : PLUMED_COLVAR_INIT(ao),
93 : pbc(true),
94 71 : do_cosine(false)
95 : {
96 142 : vector<AtomNumber> atoms,v1,v2,axis;
97 71 : parseAtomList("ATOMS",atoms);
98 71 : parseAtomList("VECTOR1",v1);
99 71 : parseAtomList("VECTOR2",v2);
100 71 : parseAtomList("AXIS",axis);
101 :
102 71 : parseFlag("COSINE",do_cosine);
103 :
104 71 : bool nopbc=!pbc;
105 71 : parseFlag("NOPBC",nopbc);
106 71 : pbc=!nopbc;
107 71 : checkRead();
108 :
109 71 : if(atoms.size()==4) {
110 67 : if(!(v1.empty() && v2.empty() && axis.empty()))
111 0 : error("ATOMS keyword is not compatible with VECTOR1, VECTOR2 and AXIS keywords");
112 67 : log.printf(" between atoms %d %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
113 67 : atoms.resize(6);
114 67 : atoms[5]=atoms[3];
115 67 : atoms[4]=atoms[2];
116 67 : atoms[3]=atoms[2];
117 67 : atoms[2]=atoms[1];
118 4 : } else if(atoms.empty()) {
119 4 : if(!(v1.size()==2 && v2.size()==2 && axis.size()==2))
120 0 : error("VECTOR1, VECTOR2 and AXIS should specify 2 atoms each");
121 : log.printf(" between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n",
122 4 : v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial());
123 4 : atoms.resize(6);
124 4 : atoms[0]=v1[1];
125 4 : atoms[1]=v1[0];
126 4 : atoms[2]=axis[0];
127 4 : atoms[3]=axis[1];
128 4 : atoms[4]=v2[0];
129 4 : atoms[5]=v2[1];
130 0 : } else error("ATOMS should specify 4 atoms");
131 :
132 71 : if(pbc) log.printf(" using periodic boundary conditions\n");
133 10 : else log.printf(" without periodic boundary conditions\n");
134 :
135 71 : if(do_cosine) log.printf(" calculating cosine instead of torsion\n");
136 :
137 71 : addValueWithDerivatives();
138 71 : if(!do_cosine) setPeriodic("-pi","pi");
139 0 : else setNotPeriodic();
140 142 : requestAtoms(atoms);
141 71 : }
142 :
143 : // calculator
144 3117 : void Torsion::calculate() {
145 :
146 3117 : Vector d0,d1,d2;
147 3117 : if(pbc) makeWhole();
148 3117 : d0=delta(getPosition(1),getPosition(0));
149 3117 : d1=delta(getPosition(3),getPosition(2));
150 3117 : d2=delta(getPosition(5),getPosition(4));
151 3117 : Vector dd0,dd1,dd2;
152 : PLMD::Torsion t;
153 3117 : double torsion=t.compute(d0,d1,d2,dd0,dd1,dd2);
154 3117 : if(do_cosine) {
155 0 : dd0 *= -sin(torsion);
156 0 : dd1 *= -sin(torsion);
157 0 : dd2 *= -sin(torsion);
158 0 : torsion = cos(torsion);
159 : }
160 3117 : setAtomsDerivatives(0,dd0);
161 3117 : setAtomsDerivatives(1,-dd0);
162 3117 : setAtomsDerivatives(2,dd1);
163 3117 : setAtomsDerivatives(3,-dd1);
164 3117 : setAtomsDerivatives(4,dd2);
165 3117 : setAtomsDerivatives(5,-dd2);
166 :
167 3117 : setValue (torsion);
168 3117 : setBoxDerivativesNoPbc();
169 3117 : }
170 :
171 : }
172 2523 : }
173 :
174 :
175 :
|
