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