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 { return false; }
99 : };
100 :
101 10423 : PLUMED_REGISTER_ACTION(Angles,"ANGLES")
102 :
103 3 : void Angles::registerKeywords( Keywords& keys ) {
104 3 : MultiColvarBase::registerKeywords( keys );
105 6 : keys.use("MEAN"); keys.use("LESS_THAN");
106 9 : keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MORE_THAN");
107 : // Could also add Region here in theory
108 6 : keys.add("numbered","ATOMS","the atoms involved in each of the angles you wish to calculate. "
109 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one angle will be "
110 : "calculated for each ATOM keyword you specify (all ATOM keywords should "
111 : "provide the indices of three atoms). The eventual number of quantities calculated by this "
112 : "action will depend on what functions of the distribution you choose to calculate.");
113 6 : keys.reset_style("ATOMS","atoms");
114 6 : keys.add("atoms-1","GROUP","Calculate angles for each distinct set of three atoms in the group");
115 6 : keys.add("atoms-2","GROUPA","A group of central atoms about which angles should be calculated");
116 6 : keys.add("atoms-2","GROUPB","When used in conjunction with GROUPA this keyword instructs plumed "
117 : "to calculate all distinct angles involving one atom from GROUPA "
118 : "and two atoms from GROUPB. The atom from GROUPA is the central atom.");
119 6 : keys.add("atoms-3","GROUPC","This must be used in conjunction with GROUPA and GROUPB. All angles "
120 : "involving one atom from GROUPA, one atom from GROUPB and one atom from "
121 : "GROUPC are calculated. The GROUPA atoms are assumed to be the central "
122 : "atoms");
123 6 : keys.add("optional","SWITCH","A switching function that ensures that only angles between atoms that "
124 : "are within a certain fixed cutoff are calculated. The following provides "
125 : "information on the \\ref switchingfunction that are available.");
126 6 : keys.add("optional","SWITCHA","A switching function on the distance between the atoms in group A and the atoms in "
127 : "group B");
128 6 : keys.add("optional","SWITCHB","A switching function on the distance between the atoms in group A and the atoms in "
129 : "group B");
130 3 : }
131 :
132 2 : Angles::Angles(const ActionOptions&ao):
133 : Action(ao),
134 : MultiColvarBase(ao),
135 2 : use_sf(false)
136 : {
137 4 : std::string sfinput,errors; parse("SWITCH",sfinput);
138 2 : if( sfinput.length()>0 ) {
139 2 : use_sf=true;
140 2 : weightHasDerivatives=true;
141 2 : sf1.set(sfinput,errors);
142 2 : if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
143 2 : sf2.set(sfinput,errors);
144 2 : if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
145 4 : log.printf(" only calculating angles for atoms separated by less than %s\n", sf1.description().c_str() );
146 : } else {
147 0 : parse("SWITCHA",sfinput);
148 0 : if(sfinput.length()>0) {
149 0 : use_sf=true;
150 0 : weightHasDerivatives=true;
151 0 : sf1.set(sfinput,errors);
152 0 : if( errors.length()!=0 ) error("problem reading SWITCHA keyword : " + errors );
153 0 : sfinput.clear(); parse("SWITCHB",sfinput);
154 0 : if(sfinput.length()==0) error("found SWITCHA keyword without SWITCHB");
155 0 : sf2.set(sfinput,errors);
156 0 : if( errors.length()!=0 ) error("problem reading SWITCHB keyword : " + errors );
157 0 : log.printf(" only calculating angles when the distance between GROUPA and GROUPB atoms is less than %s\n", sf1.description().c_str() );
158 0 : log.printf(" only calculating angles when the distance between GROUPA and GROUPC atoms is less than %s\n", sf2.description().c_str() );
159 : }
160 : }
161 : // Read in the atoms
162 : std::vector<AtomNumber> all_atoms;
163 4 : readGroupKeywords( "GROUP", "GROUPA", "GROUPB", "GROUPC", false, true, all_atoms );
164 2 : if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 3, all_atoms );
165 2 : setupMultiColvarBase( all_atoms );
166 : // Set cutoff for link cells
167 2 : if( use_sf ) {
168 2 : setLinkCellCutoff( sf1.get_dmax() );
169 2 : rcut2_1=sf1.get_dmax()*sf1.get_dmax();
170 2 : rcut2_2=sf2.get_dmax()*sf2.get_dmax();
171 : }
172 :
173 : // And check everything has been read in correctly
174 2 : checkRead();
175 : // Setup stuff for central atom
176 2 : std::vector<bool> catom_ind(3, false); catom_ind[0]=true;
177 2 : setAtomsForCentralAtom( catom_ind );
178 2 : }
179 :
180 48510 : double Angles::calculateWeight( const unsigned& taskCode, const double& weight, AtomValuePack& myatoms ) const {
181 48510 : if(!use_sf) return 1.0;
182 48510 : Vector dij=getSeparation( myatoms.getPosition(0), myatoms.getPosition(2) );
183 48510 : Vector dik=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
184 :
185 : double w1, w2, dw1, dw2, wtot;
186 48510 : double ldij = dij.modulo2(), ldik = dik.modulo2();
187 :
188 48510 : if( use_sf ) {
189 48510 : if( ldij>rcut2_1 || ldik>rcut2_2 ) return 0.0;
190 : }
191 :
192 48510 : w1=sf1.calculateSqr( ldij, dw1 );
193 48510 : w2=sf2.calculateSqr( ldik, dw2 );
194 48510 : wtot=w1*w2; dw1*=weight*w2; dw2*=weight*w1;
195 :
196 48510 : addAtomDerivatives( 0, 1, dw2*dik, myatoms );
197 48510 : addAtomDerivatives( 0, 0, -dw1*dij - dw2*dik, myatoms );
198 48510 : addAtomDerivatives( 0, 2, dw1*dij, myatoms );
199 48510 : myatoms.addBoxDerivatives( 0, (-dw1)*Tensor(dij,dij) + (-dw2)*Tensor(dik,dik) );
200 48510 : return wtot;
201 : }
202 :
203 26964 : double Angles::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
204 26964 : Vector dij=getSeparation( myatoms.getPosition(0), myatoms.getPosition(2) );
205 26964 : Vector dik=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
206 :
207 26964 : Vector ddij,ddik; PLMD::Angle a;
208 26964 : double angle=a.compute(dij,dik,ddij,ddik);
209 :
210 : // And finish the calculation
211 26964 : addAtomDerivatives( 1, 1, ddik, myatoms );
212 26964 : addAtomDerivatives( 1, 0, - ddik - ddij, myatoms );
213 26964 : addAtomDerivatives( 1, 2, ddij, myatoms );
214 26964 : myatoms.addBoxDerivatives( 1, -(Tensor(dij,ddij)+Tensor(dik,ddik)) );
215 :
216 26964 : return angle;
217 : }
218 :
219 : }
220 : }
|
