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/Angle.h"
26 : #include "tools/SwitchingFunction.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 :
31 : namespace PLMD {
32 : namespace multicolvar {
33 :
34 : //+PLUMEDOC MCOLVAR XANGLES
35 : /*
36 : Calculate the angles between the vector connecting two atoms and the x axis.
37 :
38 : \par Examples
39 :
40 : The following input tells plumed to calculate the angles between the x-axis and the vector connecting atom 3 to atom 5 and between the x-axis
41 : and the vector connecting atom 1 to atom 2. The minimum of these two quantities is then
42 : \plumedfile
43 : XANGLES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
44 : PRINT ARG=d1.min
45 : \endplumedfile
46 : (See also \ref PRINT).
47 : */
48 : //+ENDPLUMEDOC
49 :
50 : //+PLUMEDOC MCOLVAR YANGLES
51 : /*
52 : Calculate the angles between the vector connecting two atoms and the y axis.
53 :
54 : \par Examples
55 :
56 : The following input tells plumed to calculate the angles between the y-axis and the vector connecting atom 3 to atom 5 and between the y-axis
57 : and the vector connecting atom 1 to atom 2. The minimum of these two quantities is then
58 : \plumedfile
59 : YANGLES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
60 : PRINT ARG=d1.min
61 : \endplumedfile
62 : (See also \ref PRINT).
63 : */
64 : //+ENDPLUMEDOC
65 :
66 : //+PLUMEDOC MCOLVAR ZANGLES
67 : /*
68 : Calculate the angles between the vector connecting two atoms and the z axis.
69 :
70 : \par Examples
71 :
72 : The following input tells plumed to calculate the angles between the z-axis and the vector connecting atom 3 to atom 5 and between the z-axis
73 : and the vector connecting atom 1 to atom 2. The minimum of these two quantities is then
74 : \plumedfile
75 : ZANGLES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
76 : PRINT ARG=d1.min
77 : \endplumedfile
78 : (See also \ref PRINT).
79 : */
80 : //+ENDPLUMEDOC
81 :
82 :
83 :
84 : class XAngles : public MultiColvarBase {
85 : private:
86 : bool use_sf;
87 : unsigned myc;
88 : SwitchingFunction sf1;
89 : public:
90 : static void registerKeywords( Keywords& keys );
91 : explicit XAngles(const ActionOptions&);
92 : // active methods:
93 : double compute( const unsigned& tindex, AtomValuePack& myatoms ) const override;
94 : double calculateWeight( const unsigned& taskCode, const double& weight, AtomValuePack& ) const override;
95 : /// Returns the number of coordinates of the field
96 3 : bool isPeriodic() override {
97 3 : return false;
98 : }
99 : };
100 :
101 13787 : PLUMED_REGISTER_ACTION(XAngles,"XANGLES")
102 13785 : PLUMED_REGISTER_ACTION(XAngles,"YANGLES")
103 13789 : PLUMED_REGISTER_ACTION(XAngles,"ZANGLES")
104 :
105 15 : void XAngles::registerKeywords( Keywords& keys ) {
106 15 : MultiColvarBase::registerKeywords( keys );
107 15 : keys.use("MAX");
108 15 : keys.use("ALT_MIN");
109 15 : keys.use("MEAN");
110 15 : keys.use("MIN");
111 15 : keys.use("LESS_THAN");
112 15 : keys.use("LOWEST");
113 15 : keys.use("HIGHEST");
114 15 : keys.use("MORE_THAN");
115 15 : keys.use("BETWEEN");
116 15 : keys.use("HISTOGRAM");
117 15 : keys.use("MOMENTS");
118 30 : keys.add("numbered","ATOMS","the atoms involved in each of the angles you wish to calculate. "
119 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one angle will be "
120 : "calculated for each ATOM keyword you specify (all ATOM keywords should "
121 : "specify the indices of two atoms). The eventual number of quantities calculated by this "
122 : "action will depend on what functions of the distribution you choose to calculate.");
123 30 : keys.reset_style("ATOMS","atoms");
124 30 : keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
125 30 : keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
126 : "the atoms in GROUPB. This must be used in conjunction with GROUPB.");
127 30 : keys.add("atoms-2","GROUPB","Calculate the distances between all the atoms in GROUPA and all the atoms "
128 : "in GROUPB. This must be used in conjunction with GROUPA.");
129 30 : keys.add("optional","SWITCH","A switching function that ensures that only angles are only computed when atoms are within "
130 : "are within a certain fixed cutoff. The following provides information on the \\ref switchingfunction that are available.");
131 15 : }
132 :
133 3 : XAngles::XAngles(const ActionOptions&ao):
134 : Action(ao),
135 : MultiColvarBase(ao),
136 3 : use_sf(false) {
137 3 : if( getName().find("X")!=std::string::npos) {
138 1 : myc=0;
139 2 : } else if( getName().find("Y")!=std::string::npos) {
140 0 : myc=1;
141 2 : } else if( getName().find("Z")!=std::string::npos) {
142 2 : myc=2;
143 : } else {
144 0 : plumed_error();
145 : }
146 :
147 : // Read in switching function
148 : std::string sfinput, errors;
149 6 : parse("SWITCH",sfinput);
150 3 : if( sfinput.length()>0 ) {
151 2 : use_sf=true;
152 2 : weightHasDerivatives=true;
153 2 : sf1.set(sfinput,errors);
154 2 : if( errors.length()!=0 ) {
155 0 : error("problem reading SWITCH keyword : " + errors );
156 : }
157 2 : log.printf(" only calculating angles for atoms separated by less than %s\n", sf1.description().c_str() );
158 2 : setLinkCellCutoff( sf1.get_dmax() );
159 : }
160 :
161 : // Read in the atoms
162 : std::vector<AtomNumber> all_atoms;
163 6 : readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
164 3 : if( atom_lab.size()==0 ) {
165 2 : readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
166 : }
167 3 : setupMultiColvarBase( all_atoms );
168 : // And check everything has been read in correctly
169 3 : checkRead();
170 3 : }
171 :
172 110 : double XAngles::calculateWeight( const unsigned& taskCode, const double& weight, AtomValuePack& myatoms ) const {
173 110 : if(!use_sf) {
174 : return 1.0;
175 : }
176 :
177 100 : Vector distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
178 100 : double dw, w = sf1.calculateSqr( distance.modulo2(), dw );
179 100 : addAtomDerivatives( 0, 0, (-dw)*distance, myatoms );
180 100 : addAtomDerivatives( 0, 1, (+dw)*distance, myatoms );
181 100 : myatoms.addBoxDerivatives( 0, (-dw)*Tensor(distance,distance) );
182 100 : return w;
183 : }
184 :
185 50 : double XAngles::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
186 50 : Vector ddij, ddik, axis, distance;
187 50 : axis.zero();
188 50 : axis[myc]=1;
189 50 : distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
190 : PLMD::Angle a;
191 50 : double angle=a.compute( distance, axis, ddij, ddik );
192 :
193 50 : addAtomDerivatives( 1, 0, -ddij, myatoms );
194 50 : addAtomDerivatives( 1, 1, ddij, myatoms );
195 50 : myatoms.addBoxDerivatives( 1, -Tensor( distance,ddij ) );
196 50 : return angle;
197 : }
198 :
199 : }
200 : }
201 :
|
