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