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 "MultiColvarFunction.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/SwitchingFunction.h"
25 :
26 : //+PLUMEDOC MCOLVARF LOCAL_AVERAGE
27 : /*
28 : Calculate averages over spherical regions centered on atoms
29 :
30 : As is explained in <a href="http://www.youtube.com/watch?v=iDvZmbWE5ps"> this video </a> certain multicolvars
31 : calculate one scalar quantity or one vector for each of the atoms in the system. For example
32 : \ref COORDINATIONNUMBER measures the coordination number of each of the atoms in the system and \ref Q4 measures
33 : the 4th order Steinhardt parameter for each of the atoms in the system. These quantities provide tell us something about
34 : the disposition of the atoms in the first coordination sphere of each of the atoms of interest. Lechner and Dellago \cite dellago-q6
35 : have suggested that one can probe local order in a system by taking the average value of such symmetry functions over
36 : the atoms within a spherical cutoff of each of these atoms in the systems. When this is done with Steinhardt parameters
37 : they claim this gives a coordinate that is better able to distinguish solid and liquid configurations of Lennard-Jones atoms.
38 :
39 : You can calculate such locally averaged quantities within plumed by using the LOCAL_AVERAGE command. This command calculates
40 : the following atom-centered quantities:
41 :
42 : \f[
43 : s_i = \frac{ c_i + \sum_j \sigma(r_{ij})c_j }{ 1 + \sum_j \sigma(r_{ij}) }
44 : \f]
45 :
46 : where the \f$c_i\f$ and \f$c_j\f$ values can be for any one of the symmetry functions that can be calculated using plumed
47 : multicolvars. The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
48 : atoms \f$i\f$ and \f$j\f$. Lechner and Dellago suggest that the parameters of this function should be set so that it the function is equal to one
49 : when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
50 :
51 : The \f$s_i\f$ quantities calculated using the above command can be again thought of as atom-centred symmetry functions. They
52 : thus operate much like multicolvars. You can thus calculate properties of the distribution of \f$s_i\f$ values using MEAN, LESS_THAN, HISTOGRAM
53 : and so on. You can also probe the value of these averaged variables in regions of the box by using the command in tandem with the
54 : \ref AROUND command.
55 :
56 : \par Examples
57 :
58 : This example input calculates the coordination numbers for all the atoms in the system. These coordination numbers are then averaged over
59 : spherical regions. The number of averaged coordination numbers that are greater than 4 is then output to a file.
60 :
61 : \verbatim
62 : COORDINATIONNUMBER SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=d1
63 : LOCAL_AVERAGE ARG=d1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MORE_THAN={RATIONAL R_0=4} LABEL=la
64 : PRINT ARG=la.* FILE=colvar
65 : \endverbatim
66 :
67 : This example input calculates the \f$q_4\f$ (see \ref Q4) vectors for each of the atoms in the system. These vectors are then averaged
68 : component by component over a spherical region. The average value for this quantity is then outputeed to a file. This calculates the
69 : quantities that were used in the paper by Lechner and Dellago \cite dellago-q6
70 :
71 : \verbatim
72 : Q4 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2} LABEL=q4
73 : LOCAL_AVERAGE ARG=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=la
74 : PRINT ARG=la.* FILE=colvar
75 : \endverbatim
76 :
77 : */
78 : //+ENDPLUMEDOC
79 :
80 : namespace PLMD {
81 : namespace multicolvar {
82 :
83 4 : class LocalAverage : public MultiColvarFunction {
84 : private:
85 : /// Cutoff
86 : double rcut2;
87 : /// The switching function that tells us if atoms are close enough together
88 : SwitchingFunction switchingFunction;
89 : public:
90 : static void registerKeywords( Keywords& keys );
91 : explicit LocalAverage(const ActionOptions&);
92 : /// We have to overwrite this here
93 : unsigned getNumberOfQuantities() const ;
94 : /// Actually do the calculation
95 : double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
96 : /// We overwrite this in order to have dumpmulticolvar working for local average
97 0 : void normalizeVector( std::vector<double>& vals ) const {}
98 : /// Is the variable periodic
99 2 : bool isPeriodic() { return false; }
100 : };
101 :
102 2525 : PLUMED_REGISTER_ACTION(LocalAverage,"LOCAL_AVERAGE")
103 :
104 3 : void LocalAverage::registerKeywords( Keywords& keys ) {
105 3 : MultiColvarFunction::registerKeywords( keys );
106 3 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
107 3 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
108 3 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
109 3 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
110 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
111 : "The following provides information on the \\ref switchingfunction that are available. "
112 3 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
113 : // Use actionWithDistributionKeywords
114 3 : keys.use("SPECIES"); keys.use("SPECIESA"); keys.use("SPECIESB");
115 3 : keys.remove("LOWMEM"); keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("LESS_THAN");
116 3 : keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS"); keys.remove("DATA");
117 3 : keys.addFlag("LOWMEM",false,"lower the memory requirements");
118 3 : if( keys.reserved("VMEAN") ) keys.use("VMEAN");
119 3 : if( keys.reserved("VSUM") ) keys.use("VSUM");
120 3 : }
121 :
122 2 : LocalAverage::LocalAverage(const ActionOptions& ao):
123 : Action(ao),
124 2 : MultiColvarFunction(ao)
125 : {
126 2 : if( getNumberOfBaseMultiColvars()>1 ) error("local average with more than one base colvar makes no sense");
127 : // Read in the switching function
128 4 : std::string sw, errors; parse("SWITCH",sw);
129 2 : if(sw.length()>0) {
130 2 : switchingFunction.set(sw,errors);
131 : } else {
132 0 : double r_0=-1.0, d_0; int nn, mm;
133 0 : parse("NN",nn); parse("MM",mm);
134 0 : parse("R_0",r_0); parse("D_0",d_0);
135 0 : if( r_0<0.0 ) error("you must set a value for R_0");
136 0 : switchingFunction.set(nn,mm,r_0,d_0);
137 : }
138 2 : log.printf(" averaging over central molecule and those within %s\n",( switchingFunction.description() ).c_str() );
139 2 : rcut2 = switchingFunction.get_dmax()*switchingFunction.get_dmax();
140 2 : setLinkCellCutoff( switchingFunction.get_dmax() );
141 4 : std::vector<AtomNumber> all_atoms; setupMultiColvarBase( all_atoms );
142 2 : }
143 :
144 4 : unsigned LocalAverage::getNumberOfQuantities() const {
145 4 : return getBaseMultiColvar(0)->getNumberOfQuantities();
146 : }
147 :
148 128 : double LocalAverage::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
149 256 : double d2, sw, dfunc; CatomPack atom0, atom1; MultiValue& myvals = myatoms.getUnderlyingMultiValue();
150 256 : std::vector<double> values( getBaseMultiColvar(0)->getNumberOfQuantities() );
151 :
152 128 : getInputData( 0, false, myatoms, values );
153 128 : myvals.addTemporyValue( values[0] );
154 128 : if( values.size()>2 ) {
155 128 : for(unsigned j=2; j<values.size(); ++j) myatoms.addValue( j, values[0]*values[j] );
156 : } else {
157 0 : myatoms.addValue( 1, values[0]*values[1] );
158 : }
159 :
160 128 : if( !doNotCalculateDerivatives() ) {
161 128 : atom0=getCentralAtomPackFromInput( myatoms.getIndex(0) );
162 128 : MultiValue& myder=getInputDerivatives( 0, false, myatoms );
163 128 : if( values.size()>2 ) {
164 18392 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
165 18264 : unsigned jder=myder.getActiveIndex(j);
166 493128 : for(unsigned k=2; k<values.size(); ++k) {
167 474864 : myatoms.addDerivative( k, jder, values[0]*myder.getDerivative(k,jder) );
168 474864 : myatoms.addDerivative( k, jder, values[k]*myder.getDerivative(0,jder) );
169 : }
170 : }
171 : } else {
172 0 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
173 0 : unsigned jder=myder.getActiveIndex(j);
174 0 : myatoms.addDerivative( 1, jder, values[0]*myder.getDerivative(1,jder) );
175 0 : myatoms.addDerivative( 1, jder, values[1]*myder.getDerivative(0,jder) );
176 : }
177 : }
178 18392 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
179 18264 : unsigned jder=myder.getActiveIndex(j); myvals.addTemporyDerivative( jder, myder.getDerivative(0, jder) );
180 : }
181 128 : myder.clearAll();
182 : }
183 :
184 8192 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
185 8064 : Vector& distance=myatoms.getPosition(i); // getSeparation( myatoms.getPosition(0), myatoms.getPosition(i) );
186 40320 : if ( (d2=distance[0]*distance[0])<rcut2 &&
187 32256 : (d2+=distance[1]*distance[1])<rcut2 &&
188 40320 : (d2+=distance[2]*distance[2])<rcut2 &&
189 8064 : d2>epsilon) {
190 :
191 8064 : sw = switchingFunction.calculateSqr( d2, dfunc );
192 :
193 8064 : getInputData( i, false, myatoms, values );
194 8064 : if( values.size()>2 ) {
195 8064 : for(unsigned j=2; j<values.size(); ++j) myatoms.addValue( j, sw*values[0]*values[j] );
196 : } else {
197 0 : myatoms.addValue( 1, sw*values[0]*values[1] );
198 : }
199 8064 : myvals.addTemporyValue(sw);
200 :
201 8064 : if( !doNotCalculateDerivatives() ) {
202 8064 : Tensor vir(distance,distance);
203 8064 : MultiValue& myder=getInputDerivatives( i, false, myatoms );
204 8064 : atom1=getCentralAtomPackFromInput( myatoms.getIndex(i) );
205 8064 : if( values.size()>2 ) {
206 1158696 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
207 1150632 : unsigned jder=myder.getActiveIndex(j);
208 31067064 : for(unsigned k=2; k<values.size(); ++k) {
209 29916432 : myatoms.addDerivative( k, jder, sw*values[0]*myder.getDerivative(k,jder) );
210 29916432 : myatoms.addDerivative( k, jder, sw*values[k]*myder.getDerivative(0,jder) );
211 : }
212 : }
213 217728 : for(unsigned k=2; k<values.size(); ++k) {
214 209664 : myatoms.addComDerivatives( k, (-dfunc)*values[0]*values[k]*distance, atom0 );
215 209664 : myatoms.addComDerivatives( k, (+dfunc)*values[0]*values[k]*distance, atom1 );
216 209664 : myatoms.addBoxDerivatives( k, (-dfunc)*values[0]*values[k]*vir );
217 : }
218 : } else {
219 0 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
220 0 : unsigned jder=myder.getActiveIndex(j);
221 0 : myatoms.addDerivative( 1, jder, sw*values[0]*myder.getDerivative(1,jder) );
222 0 : myatoms.addDerivative( 1, jder, sw*values[1]*myder.getDerivative(0,jder) );
223 : }
224 0 : myatoms.addComDerivatives( 1, (-dfunc)*values[0]*values[1]*distance, atom0 );
225 0 : myatoms.addComDerivatives( 1, (+dfunc)*values[0]*values[1]*distance, atom1 );
226 0 : myatoms.addBoxDerivatives( 1, (-dfunc)*values[0]*values[1]*vir );
227 : }
228 : // And the bit we use to average the vector
229 8064 : myatoms.addComDerivatives( -1, (-dfunc)*values[0]*distance, atom0 );
230 8064 : myatoms.addComDerivatives( -1, (+dfunc)*values[0]*distance, atom1 );
231 1158696 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
232 1150632 : unsigned jder=myder.getActiveIndex(j); myvals.addTemporyDerivative( jder, sw*myder.getDerivative(0, jder) );
233 : }
234 8064 : myatoms.addTemporyBoxDerivatives( (-dfunc)*values[0]*vir );
235 8064 : myder.clearAll();
236 : }
237 : }
238 : }
239 :
240 : // Set the tempory weight
241 128 : updateActiveAtoms( myatoms );
242 128 : if( values.size()>2) {
243 128 : double norm=0;
244 3456 : for(unsigned i=2; i<values.size(); ++i) {
245 3328 : myvals.quotientRule( i, i );
246 : // Calculate length of vector
247 3328 : norm+=myvals.get(i)*myvals.get(i);
248 : }
249 128 : norm=sqrt(norm); myatoms.setValue(1, norm); double inorm = 1.0 / norm;
250 25856 : for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
251 25728 : unsigned jder=myvals.getActiveIndex(j);
252 694656 : for(unsigned i=2; i<values.size(); ++i) {
253 668928 : myvals.addDerivative( 1, jder, myvals.get(i)*inorm*myvals.getDerivative(i,jder) );
254 : }
255 : }
256 : } else {
257 0 : myvals.quotientRule( 1, 1 );
258 : }
259 :
260 256 : return myatoms.getValue(1);
261 : }
262 :
263 : }
264 2523 : }
|
