Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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 "MultiColvarBase.h"
23 : #include "AtomValuePack.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/Torsion.h"
26 : #include "tools/SwitchingFunction.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 :
31 : using namespace std;
32 :
33 : namespace PLMD {
34 : namespace multicolvar {
35 :
36 : //+PLUMEDOC MCOLVAR XYTORSIONS
37 : /*
38 : Calculate the torsional angle around the x axis from the positive y direction.
39 :
40 : \par Examples
41 :
42 : The following input tells plumed to calculate the angle around the x direction between the positive y-axis and the vector connecting atom 3 to atom 5 and
43 : the angle around the x direction between the positive y axis and the vector connecting atom 1 to atom 2. The minimum of these two quantities is then output
44 : \plumedfile
45 : XYTORSIONS ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
46 : PRINT ARG=d1.min
47 : \endplumedfile
48 : (See also \ref PRINT).
49 : */
50 : //+ENDPLUMEDOC
51 :
52 : //+PLUMEDOC MCOLVAR XZTORSIONS
53 : /*
54 : Calculate the torsional angle around the x axis from the positive z direction.
55 :
56 : \par Examples
57 :
58 : The following input tells plumed to calculate the angle around the x direction between the positive z-axis and the vector connecting atom 3 to atom 5 and
59 : the angle around the x direction between the positive z direction and the vector connecting atom 1 to atom 2. The minimum of these two quantities is then output
60 : \plumedfile
61 : XZTORSIONS ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
62 : PRINT ARG=d1.min
63 : \endplumedfile
64 : (See also \ref PRINT).
65 : */
66 : //+ENDPLUMEDOC
67 :
68 : //+PLUMEDOC MCOLVAR YXTORSIONS
69 : /*
70 : Calculate the torsional angle around the y axis from the positive x direction.
71 :
72 : \par Examples
73 :
74 : The following input tells plumed to calculate the angle around the y direction between the positive x-direction and the vector connecting atom 3 to atom 5 and
75 : the angle around the y direction between the positive x axis and the vector connecting atom 1 to atom 2. The minimum of these two quantities is then output
76 : \plumedfile
77 : YXTORSIONS ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
78 : PRINT ARG=d1.min
79 : \endplumedfile
80 : (See also \ref PRINT).
81 : */
82 : //+ENDPLUMEDOC
83 :
84 : //+PLUMEDOC MCOLVAR YZTORSIONS
85 : /*
86 : Calculate the torsional angle around the y axis from the positive z direction.
87 :
88 : \par Examples
89 :
90 : The following input tells plumed to calculate the angle around the y direction between the positive z-direction and the vector connecting atom 3 to atom 5 and
91 : the angle around the y direction between the positive z direction and the vector connecting atom 1 to atom 2. The minimum of these two quantities is then output
92 : \plumedfile
93 : YZTORSIONS ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
94 : PRINT ARG=d1.min
95 : \endplumedfile
96 : (See also \ref PRINT).
97 : */
98 : //+ENDPLUMEDOC
99 :
100 : //+PLUMEDOC MCOLVAR ZXTORSIONS
101 : /*
102 : Calculate the torsional angle around the z axis from the positive x direction.
103 :
104 : \par Examples
105 :
106 : The following input tells plumed to calculate the angle around the z direction between the positive x-direction and the vector connecting atom 3 to atom 5 and
107 : the angle around the z direction between the positive x-direction and the vector connecting atom 1 to atom 2. The minimum of these two quantities is then output
108 : \plumedfile
109 : ZXTORSIONS ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
110 : PRINT ARG=d1.min
111 : \endplumedfile
112 : (See also \ref PRINT).
113 : */
114 : //+ENDPLUMEDOC
115 :
116 : //+PLUMEDOC MCOLVAR ZYTORSIONS
117 : /*
118 : Calculate the torsional angle around the z axis from the positive y direction.
119 :
120 : \par Examples
121 :
122 : The following input tells plumed to calculate the angle around the z direction between the positive y-axis and the vector connecting atom 3 to atom 5 and
123 : the angle around the z direction between the positive y axis and the vector connecting atom 1 to atom 2. The minimum of these two quantities is then output
124 : \plumedfile
125 : ZYTORSIONS ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
126 : PRINT ARG=d1.min
127 : \endplumedfile
128 : (See also \ref PRINT).
129 : */
130 : //+ENDPLUMEDOC
131 :
132 :
133 :
134 :
135 6 : class XYTorsion : public MultiColvarBase {
136 : private:
137 : bool use_sf;
138 : unsigned myc1, myc2;
139 : SwitchingFunction sf1;
140 : public:
141 : static void registerKeywords( Keywords& keys );
142 : explicit XYTorsion(const ActionOptions&);
143 : // active methods:
144 : virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
145 : double calculateWeight( const unsigned& taskCode, const double& weight, AtomValuePack& ) const ;
146 : /// Returns the number of coordinates of the field
147 3 : bool isPeriodic() { return true; }
148 2 : void retrieveDomain( std::string& min, std::string& max) { min="-pi"; max="pi"; }
149 : };
150 :
151 6452 : PLUMED_REGISTER_ACTION(XYTorsion,"XYTORSIONS")
152 6452 : PLUMED_REGISTER_ACTION(XYTorsion,"XZTORSIONS")
153 6452 : PLUMED_REGISTER_ACTION(XYTorsion,"YXTORSIONS")
154 6452 : PLUMED_REGISTER_ACTION(XYTorsion,"YZTORSIONS")
155 6454 : PLUMED_REGISTER_ACTION(XYTorsion,"ZXTORSIONS")
156 6453 : PLUMED_REGISTER_ACTION(XYTorsion,"ZYTORSIONS")
157 :
158 9 : void XYTorsion::registerKeywords( Keywords& keys ) {
159 9 : MultiColvarBase::registerKeywords( keys );
160 27 : keys.use("MAX"); keys.use("ALT_MIN");
161 27 : keys.use("MEAN"); keys.use("MIN");
162 27 : keys.use("LOWEST"); keys.use("HIGHEST");
163 36 : keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
164 36 : keys.add("numbered","ATOMS","the atoms involved in each of the torsions you wish to calculate. "
165 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one torsion will be "
166 : "calculated for each ATOM keyword you specify (all ATOM keywords should "
167 : "specify the incides of two atoms). The eventual number of quantities calculated by this "
168 : "action will depend on what functions of the distribution you choose to calculate.");
169 27 : keys.reset_style("ATOMS","atoms");
170 36 : keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
171 36 : keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
172 : "the atoms in GROUPB. This must be used in conjuction with GROUPB.");
173 36 : keys.add("atoms-2","GROUPB","Calculate the distances between all the atoms in GROUPA and all the atoms "
174 : "in GROUPB. This must be used in conjuction with GROUPA.");
175 36 : keys.add("optional","SWITCH","A switching function that ensures that only angles are only computed when atoms are within "
176 : "are within a certain fixed cutoff. The following provides information on the \\ref switchingfunction that are available.");
177 9 : }
178 :
179 3 : XYTorsion::XYTorsion(const ActionOptions&ao):
180 : Action(ao),
181 : MultiColvarBase(ao),
182 3 : use_sf(false)
183 : {
184 3 : if( getName().find("XY")!=std::string::npos) { myc1=0; myc2=1; }
185 3 : else if( getName().find("XZ")!=std::string::npos) { myc1=0; myc2=2; }
186 3 : else if( getName().find("YX")!=std::string::npos) { myc1=1; myc2=0; }
187 3 : else if( getName().find("YZ")!=std::string::npos) { myc1=1; myc2=2; }
188 3 : else if( getName().find("ZX")!=std::string::npos) { myc1=2; myc2=0; }
189 1 : else if( getName().find("ZY")!=std::string::npos) { myc1=2; myc2=1; }
190 0 : else plumed_error();
191 :
192 : // Read in switching function
193 6 : std::string sfinput, errors; parse("SWITCH",sfinput);
194 3 : if( sfinput.length()>0 ) {
195 2 : use_sf=true; weightHasDerivatives=true;
196 2 : sf1.set(sfinput,errors);
197 2 : if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
198 6 : log.printf(" only calculating angles for atoms separated by less than %s\n", sf1.description().c_str() );
199 2 : setLinkCellCutoff( sf1.get_dmax() );
200 : }
201 :
202 : // Read in the atoms
203 : std::vector<AtomNumber> all_atoms;
204 12 : readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
205 4 : if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
206 3 : setupMultiColvarBase( all_atoms );
207 : // And check everything has been read in correctly
208 3 : checkRead();
209 3 : }
210 :
211 110 : double XYTorsion::calculateWeight( const unsigned& taskCode, const double& weight, AtomValuePack& myatoms ) const {
212 110 : if(!use_sf) return 1.0;
213 :
214 200 : Vector distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
215 100 : double dw, w = sf1.calculateSqr( distance.modulo2(), dw );
216 100 : addAtomDerivatives( 0, 0, (-dw)*distance, myatoms );
217 100 : addAtomDerivatives( 0, 1, (+dw)*distance, myatoms );
218 100 : myatoms.addBoxDerivatives( 0, (-dw)*Tensor(distance,distance) );
219 100 : return w;
220 : }
221 :
222 50 : double XYTorsion::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
223 50 : Vector dd0, dd1, dd2, axis, rot, distance;
224 50 : axis.zero(); rot.zero(); rot[myc1]=1; axis[myc2]=1;
225 100 : distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
226 50 : PLMD::Torsion t; double torsion=t.compute( distance, rot, axis, dd0, dd1, dd2 );
227 :
228 50 : addAtomDerivatives( 1, 0, -dd0, myatoms );
229 50 : addAtomDerivatives( 1, 1, dd0, myatoms );
230 50 : myatoms.addBoxDerivatives( 1, -extProduct(distance,dd0) );
231 50 : return torsion;
232 : }
233 :
234 : }
235 4839 : }
236 :
|
