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