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