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 "core/ActionRegister.h"
25 : #include "tools/SwitchingFunction.h"
26 :
27 : //+PLUMEDOC MCOLVARF LOCAL_AVERAGE
28 : /*
29 : Calculate averages over spherical regions centered on atoms
30 :
31 : As is explained in <a href="http://www.youtube.com/watch?v=iDvZmbWE5ps"> this video </a> certain multicolvars
32 : calculate one scalar quantity or one vector for each of the atoms in the system. For example
33 : \ref COORDINATIONNUMBER measures the coordination number of each of the atoms in the system and \ref Q4 measures
34 : the fourth order Steinhardt parameter for each of the atoms in the system. These quantities provide tell us something about
35 : the disposition of the atoms in the first coordination sphere of each of the atoms of interest. Lechner and Dellago \cite dellago-q6
36 : have suggested that one can probe local order in a system by taking the average value of such symmetry functions over
37 : the atoms within a spherical cutoff of each of these atoms in the systems. When this is done with Steinhardt parameters
38 : they claim this gives a coordinate that is better able to distinguish solid and liquid configurations of Lennard-Jones atoms.
39 :
40 : You can calculate such locally averaged quantities within plumed by using the LOCAL_AVERAGE command. This command calculates
41 : the following atom-centered quantities:
42 :
43 : \f[
44 : s_i = \frac{ c_i + \sum_j \sigma(r_{ij})c_j }{ 1 + \sum_j \sigma(r_{ij}) }
45 : \f]
46 :
47 : 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
48 : multicolvars. The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
49 : 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
50 : when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
51 :
52 : The \f$s_i\f$ quantities calculated using the above command can be again thought of as atom-centered symmetry functions. They
53 : thus operate much like multicolvars. You can thus calculate properties of the distribution of \f$s_i\f$ values using MEAN, LESS_THAN, HISTOGRAM
54 : 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
55 : \ref AROUND command.
56 :
57 : \par Examples
58 :
59 : This example input calculates the coordination numbers for all the atoms in the system. These coordination numbers are then averaged over
60 : spherical regions. The number of averaged coordination numbers that are greater than 4 is then output to a file.
61 :
62 : \plumedfile
63 : COORDINATIONNUMBER SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=d1
64 : LOCAL_AVERAGE SPECIES=d1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MORE_THAN={RATIONAL R_0=4} LABEL=la
65 : PRINT ARG=la.* FILE=colvar
66 : \endplumedfile
67 :
68 : 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
69 : component by component over a spherical region. The average value for this quantity is then output to a file. This calculates the
70 : quantities that were used in the paper by Lechner and Dellago \cite dellago-q6
71 :
72 : \plumedfile
73 : Q4 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2} LABEL=q4
74 : LOCAL_AVERAGE SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=la
75 : PRINT ARG=la.* FILE=colvar
76 : \endplumedfile
77 :
78 : */
79 : //+ENDPLUMEDOC
80 :
81 : namespace PLMD {
82 : namespace multicolvar {
83 :
84 : class LocalAverage : public MultiColvarBase {
85 : private:
86 : /// Cutoff
87 : double rcut2;
88 : /// The switching function that tells us if atoms are close enough together
89 : SwitchingFunction switchingFunction;
90 : public:
91 : static void registerKeywords( Keywords& keys );
92 : explicit LocalAverage(const ActionOptions&);
93 : /// We have to overwrite this here
94 : unsigned getNumberOfQuantities() const override;
95 : /// Actually do the calculation
96 : double compute( const unsigned& tindex, AtomValuePack& myatoms ) const override;
97 : /// We overwrite this in order to have dumpmulticolvar working for local average
98 0 : void normalizeVector( std::vector<double>& vals ) const override {}
99 : /// Is the variable periodic
100 3 : bool isPeriodic() override { return false; }
101 : };
102 :
103 10427 : PLUMED_REGISTER_ACTION(LocalAverage,"LOCAL_AVERAGE")
104 :
105 5 : void LocalAverage::registerKeywords( Keywords& keys ) {
106 5 : MultiColvarBase::registerKeywords( keys );
107 10 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
108 10 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
109 10 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
110 10 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
111 10 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
112 : "The following provides information on the \\ref switchingfunction that are available. "
113 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
114 : // Use actionWithDistributionKeywords
115 15 : keys.use("SPECIES"); keys.use("SPECIESA"); keys.use("SPECIESB");
116 20 : keys.remove("LOWMEM"); keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("LESS_THAN");
117 15 : keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
118 10 : keys.addFlag("LOWMEM",false,"lower the memory requirements");
119 15 : if( keys.reserved("VMEAN") ) keys.use("VMEAN");
120 15 : if( keys.reserved("VSUM") ) keys.use("VSUM");
121 5 : }
122 :
123 4 : LocalAverage::LocalAverage(const ActionOptions& ao):
124 : Action(ao),
125 4 : MultiColvarBase(ao)
126 : {
127 4 : if( getNumberOfBaseMultiColvars()>1 ) error("local average with more than one base colvar makes no sense");
128 : // Read in the switching function
129 8 : std::string sw, errors; parse("SWITCH",sw);
130 4 : if(sw.length()>0) {
131 4 : switchingFunction.set(sw,errors);
132 : } else {
133 0 : double r_0=-1.0, d_0; int nn, mm;
134 0 : parse("NN",nn); parse("MM",mm);
135 0 : parse("R_0",r_0); parse("D_0",d_0);
136 0 : if( r_0<0.0 ) error("you must set a value for R_0");
137 0 : switchingFunction.set(nn,mm,r_0,d_0);
138 : }
139 4 : log.printf(" averaging over central molecule and those within %s\n",( switchingFunction.description() ).c_str() );
140 4 : rcut2 = switchingFunction.get_dmax()*switchingFunction.get_dmax();
141 4 : setLinkCellCutoff( switchingFunction.get_dmax() );
142 4 : std::vector<AtomNumber> all_atoms; setupMultiColvarBase( all_atoms );
143 4 : }
144 :
145 71 : unsigned LocalAverage::getNumberOfQuantities() const {
146 71 : return getBaseMultiColvar(0)->getNumberOfQuantities();
147 : }
148 :
149 1004 : double LocalAverage::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
150 : double sw, dfunc; MultiValue& myvals = myatoms.getUnderlyingMultiValue();
151 1004 : std::vector<double> values( getBaseMultiColvar(0)->getNumberOfQuantities() );
152 :
153 1004 : getInputData( 0, false, myatoms, values );
154 : myvals.addTemporyValue( values[0] );
155 1004 : if( values.size()>2 ) {
156 27108 : for(unsigned j=2; j<values.size(); ++j) myatoms.addValue( j, values[0]*values[j] );
157 : } else {
158 0 : myatoms.addValue( 1, values[0]*values[1] );
159 : }
160 :
161 1004 : if( !doNotCalculateDerivatives() ) {
162 1004 : MultiValue& myder=getInputDerivatives( 0, false, myatoms );
163 :
164 : // Convert input atom to local index
165 : unsigned katom = myatoms.getIndex( 0 ); plumed_dbg_assert( katom<atom_lab.size() ); plumed_dbg_assert( atom_lab[katom].first>0 );
166 : // Find base colvar
167 1004 : unsigned mmc=atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
168 : // Get start of indices for this atom
169 1004 : unsigned basen=0; for(unsigned i=0; i<mmc; ++i) basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
170 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
171 :
172 1004 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
173 1004 : if( values.size()>2 ) {
174 49505 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
175 : unsigned jder=myder.getActiveIndex(j);
176 48501 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
177 39465 : unsigned kder=basen+jder;
178 1065555 : for(unsigned k=2; k<values.size(); ++k) {
179 1026090 : myatoms.addDerivative( k, kder, values[0]*myder.getDerivative(k,jder) );
180 1026090 : myatoms.addDerivative( k, kder, values[k]*myder.getDerivative(0,jder) );
181 : }
182 : } else {
183 9036 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
184 243972 : for(unsigned k=2; k<values.size(); ++k) {
185 234936 : myatoms.addDerivative( k, kder, values[0]*myder.getDerivative(k,jder) );
186 234936 : myatoms.addDerivative( k, kder, values[k]*myder.getDerivative(0,jder) );
187 : }
188 : }
189 : }
190 : } else {
191 0 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
192 : unsigned jder=myder.getActiveIndex(j);
193 0 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
194 0 : unsigned kder=basen+jder;
195 0 : myatoms.addDerivative( 1, kder, values[0]*myder.getDerivative(1,jder) );
196 0 : myatoms.addDerivative( 1, kder, values[1]*myder.getDerivative(0,jder) );
197 : } else {
198 0 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
199 0 : myatoms.addDerivative( 1, kder, values[0]*myder.getDerivative(1,jder) );
200 0 : myatoms.addDerivative( 1, kder, values[1]*myder.getDerivative(0,jder) );
201 : }
202 : }
203 : }
204 49505 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
205 : unsigned jder=myder.getActiveIndex(j);
206 48501 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
207 39465 : unsigned kder=basen+jder;
208 39465 : myvals.addTemporyDerivative( kder, myder.getDerivative(0, jder) );
209 : } else {
210 9036 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
211 9036 : myvals.addTemporyDerivative( kder, myder.getDerivative(0, jder) );
212 : }
213 : }
214 1004 : myder.clearAll();
215 : }
216 :
217 274074 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
218 : Vector& distance=myatoms.getPosition(i); // getSeparation( myatoms.getPosition(0), myatoms.getPosition(i) );
219 : double d2;
220 363595 : if ( (d2=distance[0]*distance[0])<rcut2 &&
221 90525 : (d2+=distance[1]*distance[1])<rcut2 &&
222 305390 : (d2+=distance[2]*distance[2])<rcut2 &&
223 : d2>epsilon) {
224 :
225 14639 : sw = switchingFunction.calculateSqr( d2, dfunc );
226 :
227 14639 : getInputData( i, false, myatoms, values );
228 14639 : if( values.size()>2 ) {
229 395253 : for(unsigned j=2; j<values.size(); ++j) myatoms.addValue( j, sw*values[0]*values[j] );
230 : } else {
231 0 : myatoms.addValue( 1, sw*values[0]*values[1] );
232 : }
233 : myvals.addTemporyValue(sw);
234 :
235 14639 : if( !doNotCalculateDerivatives() ) {
236 14639 : Tensor vir(distance,distance);
237 14639 : MultiValue& myder=getInputDerivatives( i, false, myatoms );
238 :
239 : // Convert input atom to local index
240 : unsigned katom = myatoms.getIndex( i ); plumed_dbg_assert( katom<atom_lab.size() ); plumed_dbg_assert( atom_lab[katom].first>0 );
241 : // Find base colvar
242 14639 : unsigned mmc=atom_lab[katom].first - 1; plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
243 : // Get start of indices for this atom
244 19559 : unsigned basen=0; for(unsigned j=0; j<mmc; ++j) basen+=mybasemulticolvars[j]->getNumberOfDerivatives() - 9;
245 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
246 :
247 14639 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
248 14639 : if( values.size()>2 ) {
249 1424993 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
250 : unsigned jder=myder.getActiveIndex(j);
251 1410354 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
252 1278603 : unsigned kder=basen+jder;
253 34522281 : for(unsigned k=2; k<values.size(); ++k) {
254 33243678 : myatoms.addDerivative( k, kder, sw*values[0]*myder.getDerivative(k,jder) );
255 33243678 : myatoms.addDerivative( k, kder, sw*values[k]*myder.getDerivative(0,jder) );
256 : }
257 : } else {
258 131751 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
259 3557277 : for(unsigned k=2; k<values.size(); ++k) {
260 3425526 : myatoms.addDerivative( k, kder, sw*values[0]*myder.getDerivative(k,jder) );
261 3425526 : myatoms.addDerivative( k, kder, sw*values[k]*myder.getDerivative(0,jder) );
262 : }
263 : }
264 : }
265 395253 : for(unsigned k=2; k<values.size(); ++k) {
266 380614 : addAtomDerivatives( k, 0, (-dfunc)*values[0]*values[k]*distance, myatoms );
267 380614 : addAtomDerivatives( k, i, (+dfunc)*values[0]*values[k]*distance, myatoms );
268 380614 : myatoms.addBoxDerivatives( k, (-dfunc)*values[0]*values[k]*vir );
269 : }
270 : } else {
271 0 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
272 : unsigned jder=myder.getActiveIndex(j);
273 0 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
274 0 : unsigned kder=basen+jder;
275 0 : myatoms.addDerivative( 1, kder, sw*values[0]*myder.getDerivative(1,jder) );
276 0 : myatoms.addDerivative( 1, kder, sw*values[1]*myder.getDerivative(0,jder) );
277 : } else {
278 0 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
279 0 : myatoms.addDerivative( 1, kder, sw*values[0]*myder.getDerivative(1,jder) );
280 0 : myatoms.addDerivative( 1, kder, sw*values[1]*myder.getDerivative(0,jder) );
281 : }
282 : }
283 0 : addAtomDerivatives( 1, 0, (-dfunc)*values[0]*values[1]*distance, myatoms );
284 0 : addAtomDerivatives( 1, i, (+dfunc)*values[0]*values[1]*distance, myatoms );
285 0 : myatoms.addBoxDerivatives( 1, (-dfunc)*values[0]*values[1]*vir );
286 : }
287 : // And the bit we use to average the vector
288 14639 : addAtomDerivatives( -1, 0, (-dfunc)*values[0]*distance, myatoms );
289 14639 : addAtomDerivatives( -1, i, (+dfunc)*values[0]*distance, myatoms );
290 1424993 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
291 : unsigned jder=myder.getActiveIndex(j);
292 1410354 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
293 1278603 : unsigned kder=basen+jder;
294 1278603 : myvals.addTemporyDerivative( kder, sw*myder.getDerivative(0, jder) );
295 : } else {
296 131751 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
297 131751 : myvals.addTemporyDerivative( kder, sw*myder.getDerivative(0, jder) );
298 : }
299 : }
300 14639 : myatoms.addTemporyBoxDerivatives( (-dfunc)*values[0]*vir );
301 14639 : myder.clearAll();
302 : }
303 : }
304 : }
305 :
306 : // Set the tempory weight
307 1004 : updateActiveAtoms( myatoms );
308 1004 : if( values.size()>2) {
309 : double norm=0;
310 27108 : for(unsigned i=2; i<values.size(); ++i) {
311 26104 : myvals.quotientRule( i, i );
312 : // Calculate length of vector
313 26104 : norm+=myvals.get(i)*myvals.get(i);
314 : }
315 1004 : norm=sqrt(norm); myatoms.setValue(1, norm); double inorm = 1.0 / norm;
316 132755 : for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
317 131751 : unsigned jder=myvals.getActiveIndex(j);
318 3557277 : for(unsigned i=2; i<values.size(); ++i) {
319 3425526 : myvals.addDerivative( 1, jder, myvals.get(i)*inorm*myvals.getDerivative(i,jder) );
320 : }
321 : }
322 : } else {
323 0 : myvals.quotientRule( 1, 1 );
324 : }
325 :
326 1004 : return myatoms.getValue(1);
327 : }
328 :
329 : }
330 : }
|
