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