Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "tools/Angle.h"
25 : #include "tools/SwitchingFunction.h"
26 : #include "core/ActionRegister.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 :
31 : namespace PLMD {
32 : namespace multicolvar {
33 :
34 : //+PLUMEDOC MCOLVAR ANGLES
35 : /*
36 : Calculate functions of the distribution of angles .
37 :
38 : You can use this command to calculate functions such as:
39 :
40 : \f[
41 : f(x) = \sum_{ijk} g( \theta_{ijk} )
42 : \f]
43 :
44 : Alternatively you can use this command to calculate functions such as:
45 :
46 : \f[
47 : f(x) = \sum_{ijk} s(r_{ij})s(r_{jk}) g(\theta_{ijk})
48 : \f]
49 :
50 : where \f$s(r)\f$ is a \ref switchingfunction. This second form means that you can
51 : use this to calculate functions of the angles in the first coordination sphere of
52 : an atom / molecule \cite lj-recon.
53 :
54 : \par Examples
55 :
56 : The following example instructs plumed to find the average of two angles and to
57 : print it to a file
58 :
59 : \plumedfile
60 : ANGLES ATOMS1=1,2,3 ATOMS2=4,5,6 MEAN LABEL=a1
61 : PRINT ARG=a1.mean FILE=colvar
62 : \endplumedfile
63 :
64 : The following example tells plumed to calculate all angles involving
65 : at least one atom from GROUPA and two atoms from GROUPB in which the distances
66 : are less than 1.0. The number of angles between \f$\frac{\pi}{4}\f$ and
67 : \f$\frac{3\pi}{4}\f$ is then output
68 :
69 : \plumedfile
70 : ANGLES GROUPA=1-10 GROUPB=11-100 BETWEEN={GAUSSIAN LOWER=0.25pi UPPER=0.75pi} SWITCH={GAUSSIAN R_0=1.0} LABEL=a1
71 : PRINT ARG=a1.between FILE=colvar
72 : \endplumedfile
73 :
74 : This final example instructs plumed to calculate all the angles in the first coordination
75 : spheres of the atoms. The bins for a normalized histogram of the distribution is then output
76 :
77 : \plumedfile
78 : ANGLES GROUP=1-38 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=pi NBINS=20} SWITCH={GAUSSIAN R_0=1.0} LABEL=a1
79 : PRINT ARG=a1.* FILE=colvar
80 : \endplumedfile
81 :
82 : */
83 : //+ENDPLUMEDOC
84 :
85 : class Angles : public MultiColvarBase {
86 : private:
87 : bool use_sf;
88 : double rcut2_1, rcut2_2;
89 : SwitchingFunction sf1;
90 : SwitchingFunction sf2;
91 : public:
92 : static void registerKeywords( Keywords& keys );
93 : explicit Angles(const ActionOptions&);
94 : /// Updates neighbor list
95 : double compute( const unsigned& tindex, AtomValuePack& ) const override;
96 : /// Returns the number of coordinates of the field
97 : double calculateWeight( const unsigned& taskCode, const double& weight, AtomValuePack& ) const override;
98 6 : bool isPeriodic() override {
99 6 : return false;
100 : }
101 : };
102 :
103 13789 : PLUMED_REGISTER_ACTION(Angles,"ANGLES")
104 :
105 6 : void Angles::registerKeywords( Keywords& keys ) {
106 6 : MultiColvarBase::registerKeywords( keys );
107 6 : keys.use("MEAN");
108 6 : keys.use("LESS_THAN");
109 6 : keys.use("BETWEEN");
110 6 : keys.use("HISTOGRAM");
111 6 : keys.use("MORE_THAN");
112 : // Could also add Region here in theory
113 12 : keys.add("numbered","ATOMS","the atoms involved in each of the angles you wish to calculate. "
114 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one angle will be "
115 : "calculated for each ATOM keyword you specify (all ATOM keywords should "
116 : "provide the indices of three atoms). The eventual number of quantities calculated by this "
117 : "action will depend on what functions of the distribution you choose to calculate.");
118 12 : keys.reset_style("ATOMS","atoms");
119 12 : keys.add("atoms-1","GROUP","Calculate angles for each distinct set of three atoms in the group");
120 12 : keys.add("atoms-2","GROUPA","A group of central atoms about which angles should be calculated");
121 12 : keys.add("atoms-2","GROUPB","When used in conjunction with GROUPA this keyword instructs plumed "
122 : "to calculate all distinct angles involving one atom from GROUPA "
123 : "and two atoms from GROUPB. The atom from GROUPA is the central atom.");
124 12 : keys.add("atoms-3","GROUPC","This must be used in conjunction with GROUPA and GROUPB. All angles "
125 : "involving one atom from GROUPA, one atom from GROUPB and one atom from "
126 : "GROUPC are calculated. The GROUPA atoms are assumed to be the central "
127 : "atoms");
128 12 : keys.add("optional","SWITCH","A switching function that ensures that only angles between atoms that "
129 : "are within a certain fixed cutoff are calculated. The following provides "
130 : "information on the \\ref switchingfunction that are available.");
131 12 : keys.add("optional","SWITCHA","A switching function on the distance between the atoms in group A and the atoms in "
132 : "group B");
133 12 : keys.add("optional","SWITCHB","A switching function on the distance between the atoms in group A and the atoms in "
134 : "group B");
135 6 : }
136 :
137 2 : Angles::Angles(const ActionOptions&ao):
138 : Action(ao),
139 : MultiColvarBase(ao),
140 2 : use_sf(false) {
141 : std::string sfinput,errors;
142 4 : parse("SWITCH",sfinput);
143 2 : if( sfinput.length()>0 ) {
144 2 : use_sf=true;
145 2 : weightHasDerivatives=true;
146 2 : sf1.set(sfinput,errors);
147 2 : if( errors.length()!=0 ) {
148 0 : error("problem reading SWITCH keyword : " + errors );
149 : }
150 2 : sf2.set(sfinput,errors);
151 2 : if( errors.length()!=0 ) {
152 0 : error("problem reading SWITCH keyword : " + errors );
153 : }
154 4 : log.printf(" only calculating angles for atoms separated by less than %s\n", sf1.description().c_str() );
155 : } else {
156 0 : parse("SWITCHA",sfinput);
157 0 : if(sfinput.length()>0) {
158 0 : use_sf=true;
159 0 : weightHasDerivatives=true;
160 0 : sf1.set(sfinput,errors);
161 0 : if( errors.length()!=0 ) {
162 0 : error("problem reading SWITCHA keyword : " + errors );
163 : }
164 : sfinput.clear();
165 0 : parse("SWITCHB",sfinput);
166 0 : if(sfinput.length()==0) {
167 0 : error("found SWITCHA keyword without SWITCHB");
168 : }
169 0 : sf2.set(sfinput,errors);
170 0 : if( errors.length()!=0 ) {
171 0 : error("problem reading SWITCHB keyword : " + errors );
172 : }
173 0 : log.printf(" only calculating angles when the distance between GROUPA and GROUPB atoms is less than %s\n", sf1.description().c_str() );
174 0 : log.printf(" only calculating angles when the distance between GROUPA and GROUPC atoms is less than %s\n", sf2.description().c_str() );
175 : }
176 : }
177 : // Read in the atoms
178 : std::vector<AtomNumber> all_atoms;
179 4 : readGroupKeywords( "GROUP", "GROUPA", "GROUPB", "GROUPC", false, true, all_atoms );
180 2 : if( atom_lab.size()==0 ) {
181 0 : readAtomsLikeKeyword( "ATOMS", 3, all_atoms );
182 : }
183 2 : setupMultiColvarBase( all_atoms );
184 : // Set cutoff for link cells
185 2 : if( use_sf ) {
186 2 : setLinkCellCutoff( sf1.get_dmax() );
187 2 : rcut2_1=sf1.get_dmax()*sf1.get_dmax();
188 2 : rcut2_2=sf2.get_dmax()*sf2.get_dmax();
189 : }
190 :
191 : // And check everything has been read in correctly
192 2 : checkRead();
193 : // Setup stuff for central atom
194 2 : std::vector<bool> catom_ind(3, false);
195 : catom_ind[0]=true;
196 2 : setAtomsForCentralAtom( catom_ind );
197 2 : }
198 :
199 48510 : double Angles::calculateWeight( const unsigned& taskCode, const double& weight, AtomValuePack& myatoms ) const {
200 48510 : if(!use_sf) {
201 : return 1.0;
202 : }
203 48510 : Vector dij=getSeparation( myatoms.getPosition(0), myatoms.getPosition(2) );
204 48510 : Vector dik=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
205 :
206 : double w1, w2, dw1, dw2, wtot;
207 48510 : double ldij = dij.modulo2(), ldik = dik.modulo2();
208 :
209 48510 : if( use_sf ) {
210 48510 : if( ldij>rcut2_1 || ldik>rcut2_2 ) {
211 : return 0.0;
212 : }
213 : }
214 :
215 48510 : w1=sf1.calculateSqr( ldij, dw1 );
216 48510 : w2=sf2.calculateSqr( ldik, dw2 );
217 48510 : wtot=w1*w2;
218 48510 : dw1*=weight*w2;
219 48510 : dw2*=weight*w1;
220 :
221 48510 : addAtomDerivatives( 0, 1, dw2*dik, myatoms );
222 48510 : addAtomDerivatives( 0, 0, -dw1*dij - dw2*dik, myatoms );
223 48510 : addAtomDerivatives( 0, 2, dw1*dij, myatoms );
224 48510 : myatoms.addBoxDerivatives( 0, (-dw1)*Tensor(dij,dij) + (-dw2)*Tensor(dik,dik) );
225 48510 : return wtot;
226 : }
227 :
228 26964 : double Angles::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
229 26964 : Vector dij=getSeparation( myatoms.getPosition(0), myatoms.getPosition(2) );
230 26964 : Vector dik=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
231 :
232 26964 : Vector ddij,ddik;
233 : PLMD::Angle a;
234 26964 : double angle=a.compute(dij,dik,ddij,ddik);
235 :
236 : // And finish the calculation
237 26964 : addAtomDerivatives( 1, 1, ddik, myatoms );
238 26964 : addAtomDerivatives( 1, 0, - ddik - ddij, myatoms );
239 26964 : addAtomDerivatives( 1, 2, ddij, myatoms );
240 26964 : myatoms.addBoxDerivatives( 1, -(Tensor(dij,ddij)+Tensor(dik,ddik)) );
241 :
242 26964 : return angle;
243 : }
244 :
245 : }
246 : }
|
