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 {
101 3 : return false;
102 : }
103 : };
104 :
105 13793 : PLUMED_REGISTER_ACTION(LocalAverage,"LOCAL_AVERAGE")
106 :
107 8 : void LocalAverage::registerKeywords( Keywords& keys ) {
108 8 : MultiColvarBase::registerKeywords( keys );
109 16 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
110 16 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
111 16 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
112 16 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
113 16 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
114 : "The following provides information on the \\ref switchingfunction that are available. "
115 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
116 : // Use actionWithDistributionKeywords
117 8 : keys.use("SPECIES");
118 8 : keys.use("SPECIESA");
119 8 : keys.use("SPECIESB");
120 8 : keys.remove("LOWMEM");
121 8 : keys.use("MEAN");
122 8 : keys.use("MORE_THAN");
123 8 : keys.use("LESS_THAN");
124 8 : keys.use("BETWEEN");
125 8 : keys.use("HISTOGRAM");
126 8 : keys.use("MOMENTS");
127 16 : keys.addFlag("LOWMEM",false,"lower the memory requirements");
128 16 : if( keys.reserved("VMEAN") ) {
129 16 : keys.use("VMEAN");
130 : }
131 16 : if( keys.reserved("VSUM") ) {
132 16 : keys.use("VSUM");
133 : }
134 8 : }
135 :
136 4 : LocalAverage::LocalAverage(const ActionOptions& ao):
137 : Action(ao),
138 4 : MultiColvarBase(ao) {
139 4 : if( getNumberOfBaseMultiColvars()>1 ) {
140 0 : error("local average with more than one base colvar makes no sense");
141 : }
142 : // Read in the switching function
143 : std::string sw, errors;
144 8 : parse("SWITCH",sw);
145 4 : if(sw.length()>0) {
146 4 : switchingFunction.set(sw,errors);
147 : } else {
148 0 : double r_0=-1.0, d_0;
149 : int nn, mm;
150 0 : parse("NN",nn);
151 0 : parse("MM",mm);
152 0 : parse("R_0",r_0);
153 0 : parse("D_0",d_0);
154 0 : if( r_0<0.0 ) {
155 0 : error("you must set a value for R_0");
156 : }
157 0 : switchingFunction.set(nn,mm,r_0,d_0);
158 : }
159 4 : log.printf(" averaging over central molecule and those within %s\n",( switchingFunction.description() ).c_str() );
160 4 : rcut2 = switchingFunction.get_dmax()*switchingFunction.get_dmax();
161 4 : setLinkCellCutoff( switchingFunction.get_dmax() );
162 : std::vector<AtomNumber> all_atoms;
163 4 : setupMultiColvarBase( all_atoms );
164 4 : }
165 :
166 71 : unsigned LocalAverage::getNumberOfQuantities() const {
167 71 : return getBaseMultiColvar(0)->getNumberOfQuantities();
168 : }
169 :
170 1004 : double LocalAverage::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
171 : double sw, dfunc;
172 : MultiValue& myvals = myatoms.getUnderlyingMultiValue();
173 1004 : std::vector<double> values( getBaseMultiColvar(0)->getNumberOfQuantities() );
174 :
175 1004 : getInputData( 0, false, myatoms, values );
176 : myvals.addTemporyValue( values[0] );
177 1004 : if( values.size()>2 ) {
178 27108 : for(unsigned j=2; j<values.size(); ++j) {
179 26104 : myatoms.addValue( j, values[0]*values[j] );
180 : }
181 : } else {
182 0 : myatoms.addValue( 1, values[0]*values[1] );
183 : }
184 :
185 1004 : if( !doNotCalculateDerivatives() ) {
186 1004 : MultiValue& myder=getInputDerivatives( 0, false, myatoms );
187 :
188 : // Convert input atom to local index
189 : unsigned katom = myatoms.getIndex( 0 );
190 : plumed_dbg_assert( katom<atom_lab.size() );
191 : plumed_dbg_assert( atom_lab[katom].first>0 );
192 : // Find base colvar
193 1004 : unsigned mmc=atom_lab[katom].first - 1;
194 : plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
195 : // Get start of indices for this atom
196 : unsigned basen=0;
197 1004 : for(unsigned i=0; i<mmc; ++i) {
198 0 : basen+=mybasemulticolvars[i]->getNumberOfDerivatives() - 9;
199 : }
200 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
201 :
202 1004 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
203 1004 : if( values.size()>2 ) {
204 53663 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
205 : unsigned jder=myder.getActiveIndex(j);
206 52659 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
207 43623 : unsigned kder=basen+jder;
208 1177821 : for(unsigned k=2; k<values.size(); ++k) {
209 1134198 : myatoms.addDerivative( k, kder, values[0]*myder.getDerivative(k,jder) );
210 1134198 : myatoms.addDerivative( k, kder, values[k]*myder.getDerivative(0,jder) );
211 : }
212 : } else {
213 9036 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
214 243972 : for(unsigned k=2; k<values.size(); ++k) {
215 234936 : myatoms.addDerivative( k, kder, values[0]*myder.getDerivative(k,jder) );
216 234936 : myatoms.addDerivative( k, kder, values[k]*myder.getDerivative(0,jder) );
217 : }
218 : }
219 : }
220 : } else {
221 0 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
222 : unsigned jder=myder.getActiveIndex(j);
223 0 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
224 0 : unsigned kder=basen+jder;
225 0 : myatoms.addDerivative( 1, kder, values[0]*myder.getDerivative(1,jder) );
226 0 : myatoms.addDerivative( 1, kder, values[1]*myder.getDerivative(0,jder) );
227 : } else {
228 0 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
229 0 : myatoms.addDerivative( 1, kder, values[0]*myder.getDerivative(1,jder) );
230 0 : myatoms.addDerivative( 1, kder, values[1]*myder.getDerivative(0,jder) );
231 : }
232 : }
233 : }
234 53663 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
235 : unsigned jder=myder.getActiveIndex(j);
236 52659 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
237 43623 : unsigned kder=basen+jder;
238 43623 : myvals.addTemporyDerivative( kder, myder.getDerivative(0, jder) );
239 : } else {
240 9036 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
241 9036 : myvals.addTemporyDerivative( kder, myder.getDerivative(0, jder) );
242 : }
243 : }
244 1004 : myder.clearAll();
245 : }
246 :
247 78226 : for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
248 : Vector& distance=myatoms.getPosition(i); // getSeparation( myatoms.getPosition(0), myatoms.getPosition(i) );
249 : double d2;
250 123529 : if ( (d2=distance[0]*distance[0])<rcut2 &&
251 46307 : (d2+=distance[1]*distance[1])<rcut2 &&
252 103834 : (d2+=distance[2]*distance[2])<rcut2 &&
253 : d2>epsilon) {
254 :
255 16025 : sw = switchingFunction.calculateSqr( d2, dfunc );
256 :
257 16025 : getInputData( i, false, myatoms, values );
258 16025 : if( values.size()>2 ) {
259 432675 : for(unsigned j=2; j<values.size(); ++j) {
260 416650 : myatoms.addValue( j, sw*values[0]*values[j] );
261 : }
262 : } else {
263 0 : myatoms.addValue( 1, sw*values[0]*values[1] );
264 : }
265 : myvals.addTemporyValue(sw);
266 :
267 16025 : if( !doNotCalculateDerivatives() ) {
268 16025 : Tensor vir(distance,distance);
269 16025 : MultiValue& myder=getInputDerivatives( i, false, myatoms );
270 :
271 : // Convert input atom to local index
272 : unsigned katom = myatoms.getIndex( i );
273 : plumed_dbg_assert( katom<atom_lab.size() );
274 : plumed_dbg_assert( atom_lab[katom].first>0 );
275 : // Find base colvar
276 16025 : unsigned mmc=atom_lab[katom].first - 1;
277 : plumed_dbg_assert( mybasemulticolvars[mmc]->taskIsCurrentlyActive( atom_lab[katom].second ) );
278 : // Get start of indices for this atom
279 : unsigned basen=0;
280 22331 : for(unsigned j=0; j<mmc; ++j) {
281 6306 : basen+=mybasemulticolvars[j]->getNumberOfDerivatives() - 9;
282 : }
283 : plumed_dbg_assert( basen%3==0 ); // Check the number of atoms is consistent with input derivatives
284 :
285 16025 : unsigned virbas = myvals.getNumberOfDerivatives()-9;
286 16025 : if( values.size()>2 ) {
287 1523879 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
288 : unsigned jder=myder.getActiveIndex(j);
289 1507854 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
290 1363629 : unsigned kder=basen+jder;
291 36817983 : for(unsigned k=2; k<values.size(); ++k) {
292 35454354 : myatoms.addDerivative( k, kder, sw*values[0]*myder.getDerivative(k,jder) );
293 35454354 : myatoms.addDerivative( k, kder, sw*values[k]*myder.getDerivative(0,jder) );
294 : }
295 : } else {
296 144225 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
297 3894075 : for(unsigned k=2; k<values.size(); ++k) {
298 3749850 : myatoms.addDerivative( k, kder, sw*values[0]*myder.getDerivative(k,jder) );
299 3749850 : myatoms.addDerivative( k, kder, sw*values[k]*myder.getDerivative(0,jder) );
300 : }
301 : }
302 : }
303 432675 : for(unsigned k=2; k<values.size(); ++k) {
304 416650 : addAtomDerivatives( k, 0, (-dfunc)*values[0]*values[k]*distance, myatoms );
305 416650 : addAtomDerivatives( k, i, (+dfunc)*values[0]*values[k]*distance, myatoms );
306 416650 : myatoms.addBoxDerivatives( k, (-dfunc)*values[0]*values[k]*vir );
307 : }
308 : } else {
309 0 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
310 : unsigned jder=myder.getActiveIndex(j);
311 0 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
312 0 : unsigned kder=basen+jder;
313 0 : myatoms.addDerivative( 1, kder, sw*values[0]*myder.getDerivative(1,jder) );
314 0 : myatoms.addDerivative( 1, kder, sw*values[1]*myder.getDerivative(0,jder) );
315 : } else {
316 0 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
317 0 : myatoms.addDerivative( 1, kder, sw*values[0]*myder.getDerivative(1,jder) );
318 0 : myatoms.addDerivative( 1, kder, sw*values[1]*myder.getDerivative(0,jder) );
319 : }
320 : }
321 0 : addAtomDerivatives( 1, 0, (-dfunc)*values[0]*values[1]*distance, myatoms );
322 0 : addAtomDerivatives( 1, i, (+dfunc)*values[0]*values[1]*distance, myatoms );
323 0 : myatoms.addBoxDerivatives( 1, (-dfunc)*values[0]*values[1]*vir );
324 : }
325 : // And the bit we use to average the vector
326 16025 : addAtomDerivatives( -1, 0, (-dfunc)*values[0]*distance, myatoms );
327 16025 : addAtomDerivatives( -1, i, (+dfunc)*values[0]*distance, myatoms );
328 1523879 : for(unsigned j=0; j<myder.getNumberActive(); ++j) {
329 : unsigned jder=myder.getActiveIndex(j);
330 1507854 : if( jder<mybasemulticolvars[mmc]->getNumberOfDerivatives()-9 ) {
331 1363629 : unsigned kder=basen+jder;
332 1363629 : myvals.addTemporyDerivative( kder, sw*myder.getDerivative(0, jder) );
333 : } else {
334 144225 : unsigned kder=virbas + (jder - mybasemulticolvars[mmc]->getNumberOfDerivatives() + 9);
335 144225 : myvals.addTemporyDerivative( kder, sw*myder.getDerivative(0, jder) );
336 : }
337 : }
338 16025 : myatoms.addTemporyBoxDerivatives( (-dfunc)*values[0]*vir );
339 16025 : myder.clearAll();
340 : }
341 : }
342 : }
343 :
344 : // Set the tempory weight
345 1004 : updateActiveAtoms( myatoms );
346 1004 : if( values.size()>2) {
347 : double norm=0;
348 27108 : for(unsigned i=2; i<values.size(); ++i) {
349 26104 : myvals.quotientRule( i, i );
350 : // Calculate length of vector
351 26104 : norm+=myvals.get(i)*myvals.get(i);
352 : }
353 1004 : norm=sqrt(norm);
354 : myatoms.setValue(1, norm);
355 1004 : double inorm = 1.0 / norm;
356 165719 : for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
357 164715 : unsigned jder=myvals.getActiveIndex(j);
358 4447305 : for(unsigned i=2; i<values.size(); ++i) {
359 4282590 : myvals.addDerivative( 1, jder, myvals.get(i)*inorm*myvals.getDerivative(i,jder) );
360 : }
361 : }
362 : } else {
363 0 : myvals.quotientRule( 1, 1 );
364 : }
365 :
366 1004 : return myatoms.getValue(1);
367 : }
368 :
369 : }
370 : }
|
