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_FAST
32 : /*
33 : Calculate the center for a group of atoms, with arbitrary weights.
34 :
35 : \par Examples
36 :
37 : */
38 : //+ENDPLUMEDOC
39 :
40 : //+PLUMEDOC VATOM CENTER
41 : /*
42 : Calculate the center for a group of atoms, with arbitrary weights.
43 :
44 : The position of the center ${r}_{\rm COM}$ given by:
45 :
46 : $$
47 : {r}_{\rm COM}=\frac{\sum_{i=1}^{n} {r}_i\ w_i }{\sum_i^{n} w_i}
48 : $$
49 :
50 : In these expressions $r_i$ indicates the position of atom $i$ and $w_i$ is the weight for atom $i$. The following input
51 : shows how you can calculate the expressions for a set of atoms by using PLUMED:
52 :
53 : ```plumed
54 : # This action calculates the position of the point on the line connecting atoms 1 and 10 that is
55 : # twice as far atom 10 as it is from atom 1
56 : c1: CENTER ATOMS=1,1,10
57 : # this is another way of calculating this position
58 : c1bis: CENTER ATOMS=1,10 WEIGHTS=2,1
59 :
60 : DUMPATOMS ATOMS=c1,c1bis FILE=atoms.xyz
61 : ```
62 :
63 : Notice that center's position is stored as [a virtual atom](specifying_atoms.md). The positions of the centers in the above input
64 : used in the DUMPATOMS command by using the labels for the CENTER actions. Notice, furthermore, that the mass and charge of this new center
65 : are equal to the sums of the masses and charges of the input atoms.
66 :
67 : 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.
68 :
69 : ```plumed
70 : c: CENTER ATOMS=1-5 MASS
71 : ```
72 :
73 : 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
74 : another action as has been done in the input below.
75 :
76 : ```plumed
77 : fcc: FCCUBIC SPECIES=1-1000 SWITCH={RATIONAL D_0=3.0 R_0=1.5}
78 : sfcc: MORE_THAN ARG=fcc SWITCH={RATIONAL R_0=0.5}
79 : c: CENTER ATOMS=1-1000 WEIGHTS=sfcc
80 : DUMPATOMS ATOMS=c FILE=atom.xyz
81 : ```
82 :
83 : This input assumes you have a cluster of solid atoms in a liquid. The actions with labels `fcc` and `sfcc`
84 : are used to differentiate between atoms in solid-like and liquid-like atoms. `sfcc` is thus a vector with
85 : one element for each atom. These elements are equal to one if the environment around the corresponding atom
86 : are solid like and zero if the environment around the atom is liquid like.
87 :
88 : ## A note on periodic boundary conditions
89 :
90 : If you run with periodic boundary conditions
91 : these are taken into account automatically when computing the center of mass. The way this is
92 : handled is akin to the way molecules are rebuilt in the [WHOLEMOLECULES](WHOLEMOLECULES.md) command.
93 : However, at variance to [WHOLEMOLECULES](WHOLEMOLECULES.md), the copies of the atomic positions in this
94 : action are modified. The global positions (i.e. those that are used in all other actions)
95 : are not changed when the alignment is performed.
96 :
97 : If you believe that PBC should not be applied when calculating the position fo the center of mass you can use the
98 : NOPBC flag.
99 :
100 : An additional way of managing periodic boundary conditions is offered in CENTER by using the PHASES keyword as shown
101 : in the example input below
102 :
103 : ```plumed
104 : c: CENTER ATOMS=1-100 PHASES
105 : ```
106 :
107 : 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$,
108 : using the following expression when the PHASES option is employed
109 :
110 : $$
111 : x_\textrm{com} = \frac{1}{2\pi} \atan\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)
112 : $$
113 :
114 : 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
115 : by multiplying these scaled components by the cell vectors. Notice that by construction this center position is
116 : not invariant with respect to rotations of the atoms at fixed cell lattice.
117 : In addition, for symmetric Bravais lattices, it is not invariant with respect
118 : to special symmetries. E.g., if you have an hexagonal cell, the center will
119 : not be invariant with respect to rotations of 120 degrees.
120 : On the other hand, it might make the treatment of PBC easier in difficult cases.
121 :
122 : */
123 : //+ENDPLUMEDOC
124 :
125 : //+PLUMEDOC VATOM COM
126 : /*
127 : Calculate the center of mass for a group of atoms.
128 :
129 : The following example input shows how you can use this shortcut to calculate the center of
130 : mass for atoms 1,2,3,4,5,6,7 and for atoms 15,20. You can then compute the distance between
131 : the two center of masses.
132 :
133 : ```plumed
134 : c1: COM ATOMS=1-7
135 : c2: COM ATOMS=15,20
136 : d1: DISTANCE ATOMS=c1,c2
137 : PRINT ARG=d1
138 : ```
139 :
140 : The center of mass is stored as a [virtual atom](specifying_atoms.md). As you can see from the
141 : above to refer to the position of the center of mass when specifying the atoms that should be used
142 : when calculating some other action you use the lable of the COM action that computes the center of mass
143 : of interest.
144 :
145 : The COM command is a shortcut because it is a wrapper to [CENTER](CENTER.md). CENTER is more powerful than
146 : comm as it allows you to use arbitrary weights in place of the masses.
147 :
148 : If you run with periodic boundary conditions
149 : these are taken into account automatically when computing the center of mass. The way this is
150 : handled is akin to the way molecules are rebuilt in the [WHOLEMOLECULES](WHOLEMOLECULES.md) command.
151 : However, at variance to [WHOLEMOLECULES](WHOLEMOLECULES.md), the copies of the atomic positions in this
152 : action are modified. The global positions (i.e. those that are used in all other actions)
153 : are not changed when the alignment is performed.
154 :
155 : If you believe that PBC should not be applied when calculating the position fo the center of mass you can use the
156 : NOPBC flag.
157 :
158 : */
159 : //+ENDPLUMEDOC
160 :
161 :
162 : class Center:
163 : public ActionWithVirtualAtom {
164 : std::vector<double> weights;
165 : std::vector<Tensor> dcenter_sin;
166 : std::vector<Tensor> dcenter_cos;
167 : std::vector<Tensor> deriv;
168 : bool isChargeSet_;
169 : bool isMassSet_;
170 : bool weight_mass;
171 : bool nopbc;
172 : bool first;
173 : bool phases;
174 : public:
175 : explicit Center(const ActionOptions&ao);
176 : void calculate() override;
177 : static void registerKeywords( Keywords& keys );
178 : };
179 :
180 : PLUMED_REGISTER_ACTION(Center,"CENTER_FAST")
181 : PLUMED_REGISTER_ACTION(Center,"COM")
182 :
183 16735 : void Center::registerKeywords(Keywords& keys) {
184 16735 : ActionWithVirtualAtom::registerKeywords(keys);
185 33470 : if( keys.getDisplayName()!="COM" ) {
186 29790 : keys.setDisplayName("CENTER");
187 : }
188 16735 : keys.add("optional","WEIGHTS","Center is computed as a weighted average.");
189 16735 : keys.add("optional","SET_CHARGE","Set the charge of the virtual atom to a given value.");
190 16735 : keys.add("optional","SET_MASS","Set the mass of the virtual atom to a given value.");
191 16735 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
192 16735 : keys.addFlag("MASS",false,"If set center is mass weighted");
193 16735 : keys.addFlag("PHASES",false,"Compute center using trigonometric phases");
194 16735 : }
195 :
196 9285 : Center::Center(const ActionOptions&ao):
197 : Action(ao),
198 : ActionWithVirtualAtom(ao),
199 9285 : isChargeSet_(false),
200 9285 : isMassSet_(false),
201 9285 : weight_mass(false),
202 9285 : nopbc(false),
203 9285 : first(true),
204 9285 : phases(false) {
205 : std::vector<AtomNumber> atoms;
206 18570 : parseAtomList("ATOMS",atoms);
207 9285 : if(atoms.size()==0) {
208 0 : error("at least one atom should be specified");
209 : }
210 9285 : parseVector("WEIGHTS",weights);
211 9285 : parseFlag("MASS",weight_mass);
212 9285 : parseFlag("NOPBC",nopbc);
213 9285 : parseFlag("PHASES",phases);
214 9285 : double charge_=std::numeric_limits<double>::lowest();
215 9285 : parse("SET_CHARGE",charge_);
216 9285 : setCharge(charge_);
217 9285 : if(charge_!=std::numeric_limits<double>::lowest()) {
218 2 : isChargeSet_=true;
219 : }
220 9285 : double mass_=-1;
221 9285 : parse("SET_MASS",mass_);
222 9285 : setMass(mass_);
223 9285 : if(mass_>0.) {
224 2 : isMassSet_=true;
225 : }
226 9285 : if(mass_==0.) {
227 0 : error("SETMASS must be greater than 0");
228 : }
229 9285 : if( getName()=="COM") {
230 1838 : weight_mass=true;
231 : }
232 9285 : checkRead();
233 9285 : log.printf(" of atoms:");
234 87728 : for(unsigned i=0; i<atoms.size(); ++i) {
235 78443 : if(i%25==0) {
236 11056 : log<<"\n";
237 : }
238 78443 : log.printf(" %d",atoms[i].serial());
239 : }
240 9285 : log<<"\n";
241 9285 : if(weight_mass) {
242 1840 : log<<" mass weighted\n";
243 1840 : if(weights.size()!=0) {
244 0 : error("WEIGHTS and MASS keywords cannot not be used simultaneously");
245 : }
246 : } else {
247 7445 : if( weights.size()==0) {
248 402 : log<<" using the geometric center\n";
249 402 : weights.resize( atoms.size() );
250 2242 : for(unsigned i=0; i<atoms.size(); i++) {
251 1840 : weights[i] = 1.;
252 : }
253 : } else {
254 7043 : log<<" with weights:";
255 7043 : if( weights.size()!=atoms.size() ) {
256 2 : error("number of elements in weight vector does not match the number of atoms");
257 : }
258 35274 : for(unsigned i=0; i<weights.size(); ++i) {
259 28232 : if(i%25==0) {
260 7042 : log<<"\n";
261 : }
262 28232 : log.printf(" %f",weights[i]);
263 : }
264 7042 : log.printf("\n");
265 : }
266 : }
267 9284 : if(phases) {
268 4 : log<<" Phases will be used to take into account PBC\n";
269 9280 : } else if(nopbc) {
270 47 : log<<" PBC will be ignored\n";
271 : } else {
272 9233 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
273 : }
274 9284 : requestAtoms(atoms);
275 9286 : }
276 :
277 32946 : void Center::calculate() {
278 32946 : Vector pos;
279 32946 : const bool dophases=(getPbc().isSet() ? phases : false);
280 :
281 32946 : if(!nopbc && !dophases) {
282 20638 : makeWhole();
283 : }
284 :
285 32946 : if( first ) {
286 15376 : if( weight_mass ) {
287 172257 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
288 164230 : if(std::isnan(getMass(i))) {
289 0 : error(
290 : "You are trying to compute a CENTER or COM but masses are not known.\n"
291 : " If you are using plumed driver, please use the --mc option"
292 : );
293 : }
294 : }
295 : }
296 : double mass(0.0);
297 209364 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
298 193988 : mass+=getMass(i);
299 : }
300 15376 : if( chargesWereSet && !isChargeSet_) {
301 : double charge(0.0);
302 192732 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
303 179414 : charge+=getCharge(i);
304 : }
305 13318 : setCharge(charge);
306 2058 : } else if( !isChargeSet_ ) {
307 2056 : setCharge(0.0);
308 : }
309 15376 : if(!isMassSet_) {
310 15374 : setMass(mass);
311 : }
312 :
313 15376 : if( weight_mass ) {
314 8027 : weights.resize( getNumberOfAtoms() );
315 172257 : for(unsigned i=0; i<weights.size(); i++) {
316 164230 : weights[i] = getMass(i) / mass;
317 : }
318 : } else {
319 : double wtot=0.0;
320 37107 : for(unsigned i=0; i<weights.size(); i++) {
321 29758 : wtot+=weights[i];
322 : }
323 37107 : for(unsigned i=0; i<weights.size(); i++) {
324 29758 : weights[i]=weights[i]/wtot;
325 : }
326 7349 : first=false;
327 : }
328 : }
329 :
330 32946 : deriv.resize(getNumberOfAtoms());
331 :
332 32946 : if(dophases) {
333 241 : dcenter_sin.resize(getNumberOfAtoms());
334 241 : dcenter_cos.resize(getNumberOfAtoms());
335 241 : Vector center_sin;
336 241 : Vector center_cos;
337 241 : Tensor invbox2pi=2*pi*getPbc().getInvBox();
338 241 : Tensor box2pi=getPbc().getBox() / (2*pi);
339 964 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
340 723 : double w=weights[i];
341 :
342 : // real to scaled
343 723 : const Vector scaled=matmul(getPosition(i),invbox2pi);
344 : const Vector ccos(
345 723 : w*std::cos(scaled[0]),
346 723 : w*std::cos(scaled[1]),
347 723 : w*std::cos(scaled[2])
348 723 : );
349 : const Vector csin(
350 723 : w*std::sin(scaled[0]),
351 723 : w*std::sin(scaled[1]),
352 723 : w*std::sin(scaled[2])
353 723 : );
354 723 : center_cos+=ccos;
355 723 : center_sin+=csin;
356 2892 : for(unsigned l=0; l<3; l++)
357 8676 : for(unsigned k=0; k<3; k++) {
358 : // k over real coordinates
359 : // l over scaled coordinates
360 6507 : dcenter_sin[i][l][k]=ccos[l]*invbox2pi[k][l];
361 6507 : dcenter_cos[i][l][k]=-csin[l]*invbox2pi[k][l];
362 : }
363 : }
364 : const Vector c(
365 241 : std::atan2(center_sin[0],center_cos[0]),
366 241 : std::atan2(center_sin[1],center_cos[1]),
367 241 : std::atan2(center_sin[2],center_cos[2])
368 241 : );
369 :
370 : // normalization is convenient for doing derivatives later
371 964 : for(unsigned l=0; l<3; l++) {
372 723 : double norm=1.0/(center_sin[l]*center_sin[l]+center_cos[l]*center_cos[l]);
373 723 : center_sin[l]*=norm;
374 723 : center_cos[l]*=norm;
375 : }
376 :
377 964 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
378 723 : Tensor dd;
379 2892 : for(unsigned l=0; l<3; l++)
380 8676 : for(unsigned k=0; k<3; k++) {
381 : // k over real coordinates
382 : // l over scaled coordinates
383 6507 : dd[l][k]= (center_cos[l]*dcenter_sin[i][l][k] - center_sin[l]*dcenter_cos[i][l][k]);
384 : }
385 : // scaled to real
386 723 : deriv[i]=matmul(dd,box2pi);
387 : }
388 241 : setAtomsDerivatives(deriv);
389 : // scaled to real
390 241 : setPosition(matmul(c,box2pi));
391 : } else {
392 412212 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
393 379507 : double w=weights[i];
394 379507 : pos+=w*getPosition(i);
395 379507 : deriv[i]=w*Tensor::identity();
396 : }
397 32705 : setPosition(pos);
398 32705 : setAtomsDerivatives(deriv);
399 : }
400 32946 : }
401 :
402 : }
403 : }
|
