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 { return false; }
97 : };
98 :
99 10421 : PLUMED_REGISTER_ACTION(XAngles,"XANGLES")
100 10419 : PLUMED_REGISTER_ACTION(XAngles,"YANGLES")
101 10423 : PLUMED_REGISTER_ACTION(XAngles,"ZANGLES")
102 :
103 6 : void XAngles::registerKeywords( Keywords& keys ) {
104 6 : MultiColvarBase::registerKeywords( keys );
105 12 : keys.use("MAX"); keys.use("ALT_MIN");
106 18 : keys.use("MEAN"); keys.use("MIN"); keys.use("LESS_THAN");
107 12 : keys.use("LOWEST"); keys.use("HIGHEST");
108 24 : keys.use("MORE_THAN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
109 12 : keys.add("numbered","ATOMS","the atoms involved in each of the angles you wish to calculate. "
110 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one angle will be "
111 : "calculated for each ATOM keyword you specify (all ATOM keywords should "
112 : "specify the indices of two atoms). The eventual number of quantities calculated by this "
113 : "action will depend on what functions of the distribution you choose to calculate.");
114 12 : keys.reset_style("ATOMS","atoms");
115 12 : keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
116 12 : keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
117 : "the atoms in GROUPB. This must be used in conjunction with GROUPB.");
118 12 : keys.add("atoms-2","GROUPB","Calculate the distances between all the atoms in GROUPA and all the atoms "
119 : "in GROUPB. This must be used in conjunction with GROUPA.");
120 12 : keys.add("optional","SWITCH","A switching function that ensures that only angles are only computed when atoms are within "
121 : "are within a certain fixed cutoff. The following provides information on the \\ref switchingfunction that are available.");
122 6 : }
123 :
124 3 : XAngles::XAngles(const ActionOptions&ao):
125 : Action(ao),
126 : MultiColvarBase(ao),
127 3 : use_sf(false)
128 : {
129 3 : if( getName().find("X")!=std::string::npos) myc=0;
130 2 : else if( getName().find("Y")!=std::string::npos) myc=1;
131 2 : else if( getName().find("Z")!=std::string::npos) myc=2;
132 0 : else plumed_error();
133 :
134 : // Read in switching function
135 6 : std::string sfinput, errors; parse("SWITCH",sfinput);
136 3 : if( sfinput.length()>0 ) {
137 2 : use_sf=true; weightHasDerivatives=true;
138 2 : sf1.set(sfinput,errors);
139 2 : if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
140 2 : log.printf(" only calculating angles for atoms separated by less than %s\n", sf1.description().c_str() );
141 2 : setLinkCellCutoff( sf1.get_dmax() );
142 : }
143 :
144 : // Read in the atoms
145 : std::vector<AtomNumber> all_atoms;
146 6 : readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
147 4 : if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
148 3 : setupMultiColvarBase( all_atoms );
149 : // And check everything has been read in correctly
150 3 : checkRead();
151 3 : }
152 :
153 110 : double XAngles::calculateWeight( const unsigned& taskCode, const double& weight, AtomValuePack& myatoms ) const {
154 110 : if(!use_sf) return 1.0;
155 :
156 100 : Vector distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
157 100 : double dw, w = sf1.calculateSqr( distance.modulo2(), dw );
158 100 : addAtomDerivatives( 0, 0, (-dw)*distance, myatoms );
159 100 : addAtomDerivatives( 0, 1, (+dw)*distance, myatoms );
160 100 : myatoms.addBoxDerivatives( 0, (-dw)*Tensor(distance,distance) );
161 100 : return w;
162 : }
163 :
164 50 : double XAngles::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
165 50 : Vector ddij, ddik, axis, distance; axis.zero(); axis[myc]=1;
166 50 : distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
167 50 : PLMD::Angle a; double angle=a.compute( distance, axis, ddij, ddik );
168 :
169 50 : addAtomDerivatives( 1, 0, -ddij, myatoms );
170 50 : addAtomDerivatives( 1, 1, ddij, myatoms );
171 50 : myatoms.addBoxDerivatives( 1, -Tensor( distance,ddij ) );
172 50 : return angle;
173 : }
174 :
175 : }
176 : }
177 :
|
