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 computed
45 : center is stored as a virtual atom that can be accessed in
46 : an atom list through the label for the CENTER action that creates it.
47 : Notice that the generated virtual atom has charge equal to the sum of the
48 : charges and mass equal to the sum of the masses. If used with the MASS flag,
49 : then it provides a result identical to \ref COM.
50 :
51 : When running with periodic boundary conditions, the atoms should be
52 : in the proper periodic image. This is done automatically since PLUMED 2.2,
53 : by considering the ordered list of atoms and rebuilding the molecule using a procedure
54 : that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
55 : rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
56 : which actually modifies the coordinates stored in PLUMED.
57 :
58 : In case you want to recover the old behavior you should use the NOPBC flag.
59 : In that case you need to take care that atoms are in the correct
60 : periodic image.
61 :
62 : \note As an experimental feature, CENTER also supports a keyword PHASES.
63 : This keyword finds the center of mass for sets of atoms that have been split by the period boundaries by computing scaled coordinates and average
64 : trigonometric functions, similarly to \ref CENTER_OF_MULTICOLVAR.
65 : Notice that by construction this center position is
66 : not invariant with respect to rotations of the atoms at fixed cell lattice.
67 : In addition, for symmetric Bravais lattices, it is not invariant with respect
68 : to special symmetries. E.g., if you have an hexagonal cell, the center will
69 : not be invariant with respect to rotations of 120 degrees.
70 : On the other hand, it might make the treatment of PBC easier in difficult cases.
71 :
72 : \par Examples
73 :
74 : \plumedfile
75 : # a point which is on the line connecting atoms 1 and 10, so that its distance
76 : # from 10 is twice its distance from 1:
77 : c1: CENTER ATOMS=1,1,10
78 : # this is another way of stating the same:
79 : c1bis: CENTER ATOMS=1,10 WEIGHTS=2,1
80 :
81 : # center of mass among these atoms:
82 : c2: CENTER ATOMS=2,3,4,5 MASS
83 :
84 : d1: DISTANCE ATOMS=c1,c2
85 :
86 : PRINT ARG=d1
87 : \endplumedfile
88 :
89 : */
90 : //+ENDPLUMEDOC
91 :
92 : //+PLUMEDOC VATOM COM
93 : /*
94 : Calculate the center of mass for a group of atoms.
95 :
96 : The computed
97 : center of mass is stored as a virtual atom that can be accessed in
98 : an atom list through the label for the COM action that creates it.
99 :
100 : For arbitrary weights (e.g. geometric center) see \ref CENTER.
101 :
102 : When running with periodic boundary conditions, the atoms should be
103 : in the proper periodic image. This is done automatically since PLUMED 2.2,
104 : by considering the ordered list of atoms and rebuilding the molecule using a procedure
105 : that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
106 : rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
107 : which actually modifies the coordinates stored in PLUMED.
108 :
109 : In case you want to recover the old behavior you should use the NOPBC flag.
110 : In that case you need to take care that atoms are in the correct
111 : periodic image.
112 :
113 : \par Examples
114 :
115 : The following input instructs plumed to print the distance between the
116 : center of mass for atoms 1,2,3,4,5,6,7 and that for atoms 15,20:
117 : \plumedfile
118 : c1: COM ATOMS=1-7
119 : c2: COM ATOMS=15,20
120 : d1: DISTANCE ATOMS=c1,c2
121 : PRINT ARG=d1
122 : \endplumedfile
123 :
124 : */
125 : //+ENDPLUMEDOC
126 :
127 :
128 : class Center:
129 : public ActionWithVirtualAtom {
130 : std::vector<double> weights;
131 : std::vector<Tensor> dcenter_sin;
132 : std::vector<Tensor> dcenter_cos;
133 : std::vector<Tensor> deriv;
134 : bool isChargeSet_;
135 : bool isMassSet_;
136 : bool weight_mass;
137 : bool nopbc;
138 : bool first;
139 : bool phases;
140 : public:
141 : explicit Center(const ActionOptions&ao);
142 : void calculate() override;
143 : static void registerKeywords( Keywords& keys );
144 : };
145 :
146 : PLUMED_REGISTER_ACTION(Center,"CENTER_FAST")
147 : PLUMED_REGISTER_ACTION(Center,"COM")
148 :
149 16735 : void Center::registerKeywords(Keywords& keys) {
150 16735 : ActionWithVirtualAtom::registerKeywords(keys);
151 33470 : if( keys.getDisplayName()!="COM" ) {
152 29790 : keys.setDisplayName("CENTER");
153 : }
154 33470 : keys.add("optional","WEIGHTS","Center is computed as a weighted average.");
155 33470 : keys.add("optional","SET_CHARGE","Set the charge of the virtual atom to a given value.");
156 33470 : keys.add("optional","SET_MASS","Set the mass of the virtual atom to a given value.");
157 33470 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
158 33470 : keys.addFlag("MASS",false,"If set center is mass weighted");
159 33470 : keys.addFlag("PHASES",false,"Compute center using trigonometric phases");
160 16735 : }
161 :
162 9285 : Center::Center(const ActionOptions&ao):
163 : Action(ao),
164 : ActionWithVirtualAtom(ao),
165 9285 : isChargeSet_(false),
166 9285 : isMassSet_(false),
167 9285 : weight_mass(false),
168 9285 : nopbc(false),
169 9285 : first(true),
170 9285 : phases(false) {
171 : std::vector<AtomNumber> atoms;
172 18570 : parseAtomList("ATOMS",atoms);
173 9285 : if(atoms.size()==0) {
174 0 : error("at least one atom should be specified");
175 : }
176 9285 : parseVector("WEIGHTS",weights);
177 9285 : parseFlag("MASS",weight_mass);
178 9285 : parseFlag("NOPBC",nopbc);
179 9285 : parseFlag("PHASES",phases);
180 9285 : double charge_=std::numeric_limits<double>::lowest();
181 9285 : parse("SET_CHARGE",charge_);
182 9285 : setCharge(charge_);
183 9285 : if(charge_!=std::numeric_limits<double>::lowest()) {
184 2 : isChargeSet_=true;
185 : }
186 9285 : double mass_=-1;
187 9285 : parse("SET_MASS",mass_);
188 9285 : setMass(mass_);
189 9285 : if(mass_>0.) {
190 2 : isMassSet_=true;
191 : }
192 9285 : if(mass_==0.) {
193 0 : error("SETMASS must be greater than 0");
194 : }
195 9285 : if( getName()=="COM") {
196 1838 : weight_mass=true;
197 : }
198 9285 : checkRead();
199 9285 : log.printf(" of atoms:");
200 87728 : for(unsigned i=0; i<atoms.size(); ++i) {
201 78443 : if(i%25==0) {
202 11056 : log<<"\n";
203 : }
204 78443 : log.printf(" %d",atoms[i].serial());
205 : }
206 9285 : log<<"\n";
207 9285 : if(weight_mass) {
208 1840 : log<<" mass weighted\n";
209 1840 : if(weights.size()!=0) {
210 0 : error("WEIGHTS and MASS keywords cannot not be used simultaneously");
211 : }
212 : } else {
213 7445 : if( weights.size()==0) {
214 402 : log<<" using the geometric center\n";
215 402 : weights.resize( atoms.size() );
216 2242 : for(unsigned i=0; i<atoms.size(); i++) {
217 1840 : weights[i] = 1.;
218 : }
219 : } else {
220 7043 : log<<" with weights:";
221 7043 : if( weights.size()!=atoms.size() ) {
222 2 : error("number of elements in weight vector does not match the number of atoms");
223 : }
224 35274 : for(unsigned i=0; i<weights.size(); ++i) {
225 28232 : if(i%25==0) {
226 7042 : log<<"\n";
227 : }
228 28232 : log.printf(" %f",weights[i]);
229 : }
230 7042 : log.printf("\n");
231 : }
232 : }
233 9284 : if(phases) {
234 4 : log<<" Phases will be used to take into account PBC\n";
235 9280 : } else if(nopbc) {
236 47 : log<<" PBC will be ignored\n";
237 : } else {
238 9233 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
239 : }
240 9284 : requestAtoms(atoms);
241 9286 : }
242 :
243 32946 : void Center::calculate() {
244 32946 : Vector pos;
245 32946 : const bool dophases=(getPbc().isSet() ? phases : false);
246 :
247 32946 : if(!nopbc && !dophases) {
248 20638 : makeWhole();
249 : }
250 :
251 32946 : if( first ) {
252 15376 : if( weight_mass ) {
253 172257 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
254 164230 : if(std::isnan(getMass(i))) {
255 0 : error(
256 : "You are trying to compute a CENTER or COM but masses are not known.\n"
257 : " If you are using plumed driver, please use the --mc option"
258 : );
259 : }
260 : }
261 : }
262 : double mass(0.0);
263 209364 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
264 193988 : mass+=getMass(i);
265 : }
266 15376 : if( chargesWereSet && !isChargeSet_) {
267 : double charge(0.0);
268 192732 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
269 179414 : charge+=getCharge(i);
270 : }
271 13318 : setCharge(charge);
272 2058 : } else if( !isChargeSet_ ) {
273 2056 : setCharge(0.0);
274 : }
275 15376 : if(!isMassSet_) {
276 15374 : setMass(mass);
277 : }
278 :
279 15376 : if( weight_mass ) {
280 8027 : weights.resize( getNumberOfAtoms() );
281 172257 : for(unsigned i=0; i<weights.size(); i++) {
282 164230 : weights[i] = getMass(i) / mass;
283 : }
284 : } else {
285 : double wtot=0.0;
286 37107 : for(unsigned i=0; i<weights.size(); i++) {
287 29758 : wtot+=weights[i];
288 : }
289 37107 : for(unsigned i=0; i<weights.size(); i++) {
290 29758 : weights[i]=weights[i]/wtot;
291 : }
292 7349 : first=false;
293 : }
294 : }
295 :
296 32946 : deriv.resize(getNumberOfAtoms());
297 :
298 32946 : if(dophases) {
299 241 : dcenter_sin.resize(getNumberOfAtoms());
300 241 : dcenter_cos.resize(getNumberOfAtoms());
301 241 : Vector center_sin;
302 241 : Vector center_cos;
303 241 : Tensor invbox2pi=2*pi*getPbc().getInvBox();
304 241 : Tensor box2pi=getPbc().getBox() / (2*pi);
305 964 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
306 723 : double w=weights[i];
307 :
308 : // real to scaled
309 723 : const Vector scaled=matmul(getPosition(i),invbox2pi);
310 : const Vector ccos(
311 723 : w*std::cos(scaled[0]),
312 723 : w*std::cos(scaled[1]),
313 723 : w*std::cos(scaled[2])
314 723 : );
315 : const Vector csin(
316 723 : w*std::sin(scaled[0]),
317 723 : w*std::sin(scaled[1]),
318 723 : w*std::sin(scaled[2])
319 723 : );
320 723 : center_cos+=ccos;
321 723 : center_sin+=csin;
322 2892 : for(unsigned l=0; l<3; l++)
323 8676 : for(unsigned k=0; k<3; k++) {
324 : // k over real coordinates
325 : // l over scaled coordinates
326 6507 : dcenter_sin[i][l][k]=ccos[l]*invbox2pi[k][l];
327 6507 : dcenter_cos[i][l][k]=-csin[l]*invbox2pi[k][l];
328 : }
329 : }
330 : const Vector c(
331 241 : std::atan2(center_sin[0],center_cos[0]),
332 241 : std::atan2(center_sin[1],center_cos[1]),
333 241 : std::atan2(center_sin[2],center_cos[2])
334 241 : );
335 :
336 : // normalization is convenient for doing derivatives later
337 964 : for(unsigned l=0; l<3; l++) {
338 723 : double norm=1.0/(center_sin[l]*center_sin[l]+center_cos[l]*center_cos[l]);
339 723 : center_sin[l]*=norm;
340 723 : center_cos[l]*=norm;
341 : }
342 :
343 964 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
344 723 : Tensor dd;
345 2892 : for(unsigned l=0; l<3; l++)
346 8676 : for(unsigned k=0; k<3; k++) {
347 : // k over real coordinates
348 : // l over scaled coordinates
349 6507 : dd[l][k]= (center_cos[l]*dcenter_sin[i][l][k] - center_sin[l]*dcenter_cos[i][l][k]);
350 : }
351 : // scaled to real
352 723 : deriv[i]=matmul(dd,box2pi);
353 : }
354 241 : setAtomsDerivatives(deriv);
355 : // scaled to real
356 241 : setPosition(matmul(c,box2pi));
357 : } else {
358 412212 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
359 379507 : double w=weights[i];
360 379507 : pos+=w*getPosition(i);
361 379507 : deriv[i]=w*Tensor::identity();
362 : }
363 32705 : setPosition(pos);
364 32705 : setAtomsDerivatives(deriv);
365 : }
366 32946 : }
367 :
368 : }
369 : }
|
