Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "ActionWithVirtualAtom.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include <cmath>
26 : #include <limits>
27 :
28 : namespace PLMD {
29 : namespace vatom {
30 :
31 : //+PLUMEDOC VATOM CENTER
32 : /*
33 : Calculate the center for a group of atoms, with arbitrary weights.
34 :
35 : The position of the center ${r}_{\rm COM}$ given by:
36 :
37 : $$
38 : {r}_{\rm COM}=\frac{\sum_{i=1}^{n} {r}_i\ w_i }{\sum_i^{n} w_i}
39 : $$
40 :
41 : In these expressions $r_i$ indicates the position of atom $i$ and $w_i$ is the weight for atom $i$. The following input
42 : shows how you can calculate the expressions for a set of atoms by using PLUMED:
43 :
44 : ```plumed
45 : # This action calculates the position of the point on the line connecting atoms 1 and 10 that is
46 : # twice as far atom 10 as it is from atom 1
47 : c1: CENTER ATOMS=1,1,10
48 : # this is another way of calculating this position
49 : c1bis: CENTER ATOMS=1,10 WEIGHTS=2,1
50 :
51 : DUMPATOMS ATOMS=c1,c1bis FILE=atoms.xyz
52 : ```
53 :
54 : Notice that center's position is stored as [a virtual atom](specifying_atoms.md). The positions of the centers in the above input
55 : used in the DUMPATOMS command by using the labels for the CENTER actions. Notice, furthermore, that the mass and charge of this new center
56 : are equal to the sums of the masses and charges of the input atoms.
57 :
58 : The input below shows how you can use the CENTER action in place of the [COM](COM.md) action to calculate the center of mass for a group of atoms.
59 :
60 : ```plumed
61 : c: CENTER ATOMS=1-5 MASS
62 : ```
63 :
64 : Center is more powerful than COM because you can use arbitrary vectors of weights as in the first example above or vector of weights that are calculated by
65 : another action as has been done in the input below.
66 :
67 : ```plumed
68 : fcc: FCCUBIC SPECIES=1-1000 SWITCH={RATIONAL D_0=3.0 R_0=1.5}
69 : sfcc: MORE_THAN ARG=fcc SWITCH={RATIONAL R_0=0.5}
70 : c: CENTER ATOMS=1-1000 WEIGHTS=sfcc
71 : DUMPATOMS ATOMS=c FILE=atom.xyz
72 : ```
73 :
74 : This input assumes you have a cluster of solid atoms in a liquid. The actions with labels `fcc` and `sfcc`
75 : are used to differentiate between atoms in solid-like and liquid-like atoms. `sfcc` is thus a vector with
76 : one element for each atom. These elements are equal to one if the environment around the corresponding atom
77 : are solid like and zero if the environment around the atom is liquid like.
78 :
79 : ## A note on periodic boundary conditions
80 :
81 : If you run with periodic boundary conditions
82 : these are taken into account automatically when computing the center of mass. The way this is
83 : handled is akin to the way molecules are rebuilt in the [WHOLEMOLECULES](WHOLEMOLECULES.md) command.
84 : However, at variance to [WHOLEMOLECULES](WHOLEMOLECULES.md), the copies of the atomic positions in this
85 : action are modified. The global positions (i.e. those that are used in all other actions)
86 : are not changed when the alignment is performed.
87 :
88 : If you believe that PBC should not be applied when calculating the position fo the center of mass you can use the
89 : NOPBC flag as shown below:
90 :
91 : ```plumed
92 : c: CENTER ATOMS=1-100 NOPBC
93 : ```
94 :
95 : An additional way of managing periodic boundary conditions is offered in CENTER by using the PHASES keyword as shown
96 : in the example input below
97 :
98 : ```plumed
99 : c: CENTER ATOMS=1-100 PHASES
100 : ```
101 :
102 : The scaled value for the $x$ component of the position of the center is calculated from the scaled components of the input atoms, $x_i$,
103 : using the following expression when the PHASES option is employed
104 :
105 : $$
106 : x_\textrm{com} = \frac{1}{2\pi} \arctan\left( \frac{\sum_{i=1}^n w_i \sin(2 \pi x_i ) }{ \sum_{i=1}^n w_i \cos(2 \pi x_i ) } \right)
107 : $$
108 :
109 : Similar, expressions are used to calculae the values of the scaled $y$ and $z$ components. The final cartesian coordinates of the center are then computed
110 : by multiplying these scaled components by the cell vectors. Notice that by construction this center position is
111 : not invariant with respect to rotations of the atoms at fixed cell lattice.
112 : In addition, for symmetric Bravais lattices, it is not invariant with respect
113 : to special symmetries. E.g., if you have an hexagonal cell, the center will
114 : not be invariant with respect to rotations of 120 degrees.
115 : On the other hand, it might make the treatment of PBC easier in difficult cases.
116 :
117 : As an alternative to PHASES you can use the SAFE_PHASES flag as shown below:
118 :
119 : ```plumed
120 : c: CENTER ATOMS=1-100 SAFE_PHASES
121 : ```
122 :
123 : This option will use the method described above for the PHASES keyword if the cell coordinates are set. If, however, the cell coordinates are not set the
124 : position of the center will be calculated in the usual way. We wouldn't recommend using this option in place of PHASES as additional computational overhead
125 : is introduced. Furthermore, you normally know if the cell parameters are not passed to PLUMED in advance of doing the calculation
126 :
127 : */
128 : //+ENDPLUMEDOC
129 :
130 : //+PLUMEDOC VATOM COM
131 : /*
132 : Calculate the center of mass for a group of atoms.
133 :
134 : The following example input shows how you can use this shortcut to calculate the center of
135 : mass for atoms 1,2,3,4,5,6,7 and for atoms 15,20. You can then compute the distance between
136 : the two center of masses.
137 :
138 : ```plumed
139 : c1: COM ATOMS=1-7
140 : c2: COM ATOMS=15,20
141 : d1: DISTANCE ATOMS=c1,c2
142 : PRINT ARG=d1
143 : ```
144 :
145 : The center of mass is stored as a [virtual atom](specifying_atoms.md). As you can see from the
146 : above to refer to the position of the center of mass when specifying the atoms that should be used
147 : when calculating some other action you use the lable of the COM action that computes the center of mass
148 : of interest.
149 :
150 : The COM command is a shortcut because it is a wrapper to [CENTER](CENTER.md). CENTER is more powerful than
151 : comm as it allows you to use arbitrary weights in place of the masses.
152 :
153 : ## Periodic boundary conditions
154 :
155 : If you run with periodic boundary conditions
156 : these are taken into account automatically when computing the center of mass. By default the way this is
157 : handled is akin to the way molecules are rebuilt in the [WHOLEMOLECULES](WHOLEMOLECULES.md) command.
158 : However, at variance to [WHOLEMOLECULES](WHOLEMOLECULES.md), the copies of the atomic positions in this
159 : action are modified. The global positions (i.e. those that are used in all other actions)
160 : are not changed when the alignment is performed.
161 :
162 : If you believe that PBC should not be applied when calculating the position fo the center of mass you can use the
163 : NOPBC flag as shown below:
164 :
165 : ```plumed
166 : c1: COM ATOMS=1-7 NOPBC
167 : ```
168 :
169 : Alternatively, if you would like to calculate the position of the center of mass using the PHASES method that is dicussed
170 : in the documentation for [CENTER](CENTER.md) you can add the keyword PHASES as shwo below:
171 :
172 : ```plumed
173 : c1: COM ATOMS=1-7 PHASES
174 : ```
175 :
176 : ## The mass and charge of the vatom
177 :
178 : By default the mass and charge of the COM are calculated by adding the masses and charges of the input atoms together. In other
179 : words, if your COM is calculated from the positions of $N$ atoms, the mass of the vatom is calculated as:
180 :
181 : $$
182 : m_\textrm{com} = \sum_{i=1}^N m_i
183 : $$
184 :
185 : where $m_i$ is the mass of the $i$th input atom. The charge is then calculated as:
186 :
187 : $$
188 : q_\textrm{com} = \sum_{i=1}^N q_i
189 : $$
190 :
191 : where $q_i$ is the charge of the $i$th atom. If for any reason you don't want the mass and charge of the VATOM to be set equal to these
192 : values you can set them manually by using the `SET_CHARGE` and `SET_MASS` keywords as shown below:
193 :
194 : ```plumed
195 : c1: COM ATOMS=1-7 SET_MASS=1 SET_CHARGE=-1
196 : ```
197 :
198 : */
199 : //+ENDPLUMEDOC
200 :
201 :
202 : class Center:
203 : public ActionWithVirtualAtom {
204 : std::vector<double> weights;
205 : std::vector<Tensor> dcenter_sin;
206 : std::vector<Tensor> dcenter_cos;
207 : std::vector<Tensor> deriv;
208 : bool isChargeSet_;
209 : bool isMassSet_;
210 : bool weight_mass;
211 : bool nopbc;
212 : bool first;
213 : bool phases;
214 : // true if the masses and charges of all the atoms are declared constant (e.g. they come from the MD engine)
215 : bool constantMassCharge;
216 : // computes mass and charge of the virtual atom and, if needed, the mass-based weights
217 : void updateMassAndCharge();
218 : public:
219 : explicit Center(const ActionOptions&ao);
220 : void calculate() override;
221 : static void registerKeywords( Keywords& keys );
222 : };
223 :
224 : PLUMED_REGISTER_ACTION(Center,"CENTER_FAST")
225 : PLUMED_REGISTER_ACTION(Center,"COM")
226 :
227 17240 : void Center::registerKeywords(Keywords& keys) {
228 17240 : ActionWithVirtualAtom::registerKeywords(keys);
229 17240 : if( keys.getDisplayName()!="COM" ) {
230 30794 : keys.setDisplayName("CENTER");
231 15397 : keys.add("optional","WEIGHTS","Center is computed as a weighted average.");
232 15397 : keys.addFlag("MASS",false,"If set center is mass weighted");
233 : }
234 17240 : keys.add("optional","SET_CHARGE","Set the charge of the virtual atom to a given value.");
235 17240 : keys.add("optional","SET_MASS","Set the mass of the virtual atom to a given value.");
236 17240 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
237 17240 : keys.addFlag("PHASES",false,"Compute center using trigonometric phases");
238 17240 : }
239 :
240 9536 : Center::Center(const ActionOptions&ao):
241 : Action(ao),
242 : ActionWithVirtualAtom(ao),
243 9536 : isChargeSet_(false),
244 9536 : isMassSet_(false),
245 9536 : weight_mass(false),
246 9536 : nopbc(false),
247 9536 : first(true),
248 9536 : phases(false),
249 9536 : constantMassCharge(true) {
250 9536 : std::vector<AtomNumber> atoms;
251 19072 : parseAtomList("ATOMS",atoms);
252 9536 : if(atoms.size()==0) {
253 0 : error("at least one atom should be specified");
254 : }
255 9536 : if( getName()!="COM") {
256 7697 : parseVector("WEIGHTS",weights);
257 15394 : parseFlag("MASS",weight_mass);
258 : }
259 9536 : parseFlag("NOPBC",nopbc);
260 9536 : parseFlag("PHASES",phases);
261 9536 : double charge_=std::numeric_limits<double>::lowest();
262 9536 : parse("SET_CHARGE",charge_);
263 9536 : setCharge(charge_);
264 9536 : if(charge_!=std::numeric_limits<double>::lowest()) {
265 3 : isChargeSet_=true;
266 : }
267 : const double massNotSet=std::numeric_limits<double>::lowest();
268 9536 : double mass_=massNotSet;
269 9536 : parse("SET_MASS",mass_);
270 9536 : if(mass_!=massNotSet) {
271 3 : if(!(mass_>0.)) {
272 0 : error("SET_MASS must be greater than 0");
273 : }
274 3 : isMassSet_=true;
275 3 : setMass(mass_);
276 : } else {
277 9533 : setMass(-1);
278 : }
279 9536 : if( getName()=="COM") {
280 1839 : weight_mass=true;
281 : }
282 9536 : checkRead();
283 9536 : log.printf(" of atoms:");
284 89781 : for(unsigned i=0; i<atoms.size(); ++i) {
285 80245 : if(i%25==0) {
286 11352 : log<<"\n";
287 : }
288 80245 : log.printf(" %d",atoms[i].serial());
289 : }
290 9536 : log<<"\n";
291 9536 : if(weight_mass) {
292 1844 : log<<" mass weighted\n";
293 1844 : if(weights.size()!=0) {
294 0 : error("WEIGHTS and MASS keywords cannot not be used simultaneously");
295 : }
296 : } else {
297 7692 : if( weights.size()==0) {
298 648 : log<<" using the geometric center\n";
299 648 : weights.resize( atoms.size() );
300 4275 : for(unsigned i=0; i<atoms.size(); i++) {
301 3627 : weights[i] = 1.;
302 : }
303 : } else {
304 7044 : log<<" with weights:";
305 7044 : if( weights.size()!=atoms.size() ) {
306 1 : error("number of elements in weight vector does not match the number of atoms");
307 : }
308 35278 : for(unsigned i=0; i<weights.size(); ++i) {
309 28235 : if(i%25==0) {
310 7043 : log<<"\n";
311 : }
312 28235 : log.printf(" %f",weights[i]);
313 : }
314 7043 : log.printf("\n");
315 : }
316 : // weights are constant: they are normalized once here
317 : double wtot=0.0;
318 39553 : for(unsigned i=0; i<weights.size(); i++) {
319 31862 : wtot+=weights[i];
320 : }
321 7691 : if(wtot==0.0) {
322 0 : error("the sum of the weights is zero");
323 : }
324 39553 : for(unsigned i=0; i<weights.size(); i++) {
325 31862 : weights[i]/=wtot;
326 : }
327 : }
328 9535 : if(phases) {
329 9 : log<<" Phases will be used to take into account PBC\n";
330 9 : if(nopbc) {
331 1 : log<<" WARNING: NOPBC has no effect when PHASES is used\n";
332 : }
333 9526 : } else if(nopbc) {
334 50 : log<<" PBC will be ignored\n";
335 : } else {
336 9476 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
337 : }
338 : // masses and charges coming from the MD engine are declared constant, and in that case
339 : // they are processed only once. Otherwise (e.g. center of virtual atoms) at every step.
340 9535 : std::vector<Value*> xv, yv, zv, mv, cv;
341 9535 : getAtomValuesFromPlumedObject( plumed, xv, yv, zv, mv, cv );
342 89776 : for(unsigned i=0; i<atoms.size(); ++i) {
343 80241 : const std::pair<std::size_t,std::size_t> p = getValueIndices( atoms[i] );
344 80241 : if( !mv[p.first]->isConstant() || !cv[p.first]->isConstant() ) {
345 22 : constantMassCharge=false;
346 : }
347 : }
348 9535 : requestAtoms(atoms);
349 9542 : }
350 :
351 9406 : void Center::updateMassAndCharge() {
352 : const std::vector<double>& m=getMasses();
353 : double mass(0.0);
354 89187 : for(unsigned i=0; i<m.size(); i++) {
355 79781 : mass+=m[i];
356 : }
357 9406 : if( weight_mass ) {
358 : // a NaN in any mass propagates to the sum
359 1838 : if(std::isnan(mass)) {
360 0 : error(
361 : "You are trying to compute a CENTER or COM but masses are not known.\n"
362 : " If you are using plumed driver, please use the --mc option"
363 : );
364 : }
365 1838 : if(mass==0.0) {
366 0 : error("You are trying to compute a CENTER or COM but the total mass of the atoms is zero");
367 : }
368 1838 : weights.resize( m.size() );
369 50140 : for(unsigned i=0; i<weights.size(); i++) {
370 48302 : weights[i] = m[i] / mass;
371 : }
372 : }
373 9406 : if( !isMassSet_ ) {
374 9403 : setMass(mass);
375 : }
376 9406 : if( !isChargeSet_ ) {
377 9403 : if( chargesWereSet ) {
378 9399 : const std::vector<double>& q=getCharges();
379 : double charge(0.0);
380 89135 : for(unsigned i=0; i<q.size(); i++) {
381 79736 : charge+=q[i];
382 : }
383 9399 : setCharge(charge);
384 : } else {
385 4 : setCharge(0.0);
386 : }
387 : }
388 9406 : }
389 :
390 37057 : void Center::calculate() {
391 37057 : Vector pos;
392 37057 : const bool dophases=(getPbc().isSet() ? phases : false);
393 :
394 37057 : if(!nopbc && !dophases) {
395 24707 : makeWhole();
396 : }
397 :
398 37057 : if( first ) {
399 9402 : first=false;
400 9402 : if( !weight_mass && !isMassSet_ && !massesWereSet ) {
401 0 : log.printf("WARNING: masses were not passed to plumed\n");
402 : }
403 9402 : updateMassAndCharge();
404 27655 : } else if( !constantMassCharge ) {
405 4 : updateMassAndCharge();
406 : }
407 :
408 37057 : deriv.resize(getNumberOfAtoms());
409 :
410 37057 : if(dophases) {
411 256 : dcenter_sin.resize(getNumberOfAtoms());
412 256 : dcenter_cos.resize(getNumberOfAtoms());
413 256 : Vector center_sin;
414 256 : Vector center_cos;
415 256 : Tensor invbox2pi=2*pi*getPbc().getInvBox();
416 256 : Tensor box2pi=getPbc().getBox() / (2*pi);
417 1024 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
418 768 : double w=weights[i];
419 :
420 : // real to scaled
421 768 : const Vector scaled=matmul(getPosition(i),invbox2pi);
422 : const Vector ccos(
423 768 : w*std::cos(scaled[0]),
424 768 : w*std::cos(scaled[1]),
425 768 : w*std::cos(scaled[2])
426 768 : );
427 : const Vector csin(
428 768 : w*std::sin(scaled[0]),
429 768 : w*std::sin(scaled[1]),
430 768 : w*std::sin(scaled[2])
431 768 : );
432 : center_cos+=ccos;
433 : center_sin+=csin;
434 3072 : for(unsigned l=0; l<3; l++)
435 9216 : for(unsigned k=0; k<3; k++) {
436 : // k over real coordinates
437 : // l over scaled coordinates
438 6912 : dcenter_sin[i][l][k]=ccos[l]*invbox2pi[k][l];
439 6912 : dcenter_cos[i][l][k]=-csin[l]*invbox2pi[k][l];
440 : }
441 : }
442 : const Vector c(
443 : std::atan2(center_sin[0],center_cos[0]),
444 : std::atan2(center_sin[1],center_cos[1]),
445 : std::atan2(center_sin[2],center_cos[2])
446 256 : );
447 :
448 : // normalization is convenient for doing derivatives later
449 1024 : for(unsigned l=0; l<3; l++) {
450 768 : const double den=center_sin[l]*center_sin[l]+center_cos[l]*center_cos[l];
451 768 : if(den==0.0) {
452 0 : error("the center cannot be computed with PHASES: the weighted sum of sines and cosines vanishes along one of the cell vectors");
453 : }
454 768 : double norm=1.0/den;
455 768 : center_sin[l]*=norm;
456 768 : center_cos[l]*=norm;
457 : }
458 :
459 1024 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
460 768 : Tensor dd;
461 3072 : for(unsigned l=0; l<3; l++)
462 9216 : for(unsigned k=0; k<3; k++) {
463 : // k over real coordinates
464 : // l over scaled coordinates
465 6912 : dd[l][k]= (center_cos[l]*dcenter_sin[i][l][k] - center_sin[l]*dcenter_cos[i][l][k]);
466 : }
467 : // dd[l][k] is the derivative of the scaled component l with respect to the real coordinate k,
468 : // whereas deriv[i][k][m] must be the derivative of the real component m with respect to the
469 : // real coordinate k: dd has to be transposed before going from scaled to real
470 : // (this matters only for non-orthorhombic cells)
471 768 : deriv[i]=matmul(dd.transpose(),box2pi);
472 : }
473 256 : setAtomsDerivatives(deriv);
474 : // scaled to real
475 256 : setPosition(matmul(c,box2pi));
476 : } else {
477 438647 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
478 401846 : double w=weights[i];
479 401846 : pos+=w*getPosition(i);
480 401846 : deriv[i]=w*Tensor::identity();
481 : }
482 36801 : setPosition(pos);
483 36801 : setAtomsDerivatives(deriv);
484 : }
485 37057 : }
486 :
487 : }
488 : }
|
