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 : {
131 : std::vector<double> weights;
132 : std::vector<Tensor> dcenter_sin;
133 : std::vector<Tensor> dcenter_cos;
134 : std::vector<Tensor> deriv;
135 : bool isChargeSet_;
136 : bool isMassSet_;
137 : bool weight_mass;
138 : bool nopbc;
139 : bool first;
140 : bool phases;
141 : public:
142 : explicit Center(const ActionOptions&ao);
143 : void calculate() override;
144 : static void registerKeywords( Keywords& keys );
145 : };
146 :
147 : PLUMED_REGISTER_ACTION(Center,"CENTER_FAST")
148 : PLUMED_REGISTER_ACTION(Center,"COM")
149 :
150 17317 : void Center::registerKeywords(Keywords& keys) {
151 17317 : ActionWithVirtualAtom::registerKeywords(keys);
152 50111 : if( keys.getDisplayName()!="COM" ) keys.setDisplayName("CENTER");
153 34634 : keys.add("optional","WEIGHTS","Center is computed as a weighted average.");
154 34634 : keys.add("optional","SET_CHARGE","Set the charge of the virtual atom to a given value.");
155 34634 : keys.add("optional","SET_MASS","Set the mass of the virtual atom to a given value.");
156 34634 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
157 34634 : keys.addFlag("MASS",false,"If set center is mass weighted");
158 34634 : keys.addFlag("PHASES",false,"Compute center using trigonometric phases");
159 17317 : }
160 :
161 9576 : Center::Center(const ActionOptions&ao):
162 : Action(ao),
163 : ActionWithVirtualAtom(ao),
164 9576 : isChargeSet_(false),
165 9576 : isMassSet_(false),
166 9576 : weight_mass(false),
167 9576 : nopbc(false),
168 9576 : first(true),
169 9576 : phases(false)
170 : {
171 : std::vector<AtomNumber> atoms;
172 19152 : parseAtomList("ATOMS",atoms);
173 9576 : if(atoms.size()==0) error("at least one atom should be specified");
174 9576 : parseVector("WEIGHTS",weights);
175 9576 : parseFlag("MASS",weight_mass);
176 9576 : parseFlag("NOPBC",nopbc);
177 9576 : parseFlag("PHASES",phases);
178 19152 : double charge_=std::numeric_limits<double>::lowest(); parse("SET_CHARGE",charge_); setCharge(charge_);
179 9576 : if(charge_!=std::numeric_limits<double>::lowest()) isChargeSet_=true;
180 19152 : double mass_=-1; parse("SET_MASS",mass_); setMass(mass_);
181 9576 : if(mass_>0.) isMassSet_=true;
182 9576 : if(mass_==0.) error("SETMASS must be greater than 0");
183 9576 : if( getName()=="COM") weight_mass=true;
184 9576 : checkRead();
185 9576 : log.printf(" of atoms:");
186 88601 : for(unsigned i=0; i<atoms.size(); ++i) {
187 79025 : if(i%25==0) log<<"\n";
188 79025 : log.printf(" %d",atoms[i].serial());
189 : }
190 9576 : log<<"\n";
191 9576 : if(weight_mass) {
192 1840 : log<<" mass weighted\n";
193 1840 : if(weights.size()!=0) error("WEIGHTS and MASS keywords cannot not be used simultaneously");
194 : } else {
195 7736 : if( weights.size()==0) {
196 693 : log<<" using the geometric center\n";
197 693 : weights.resize( atoms.size() );
198 3115 : for(unsigned i=0; i<atoms.size(); i++) weights[i] = 1.;
199 : } else {
200 7043 : log<<" with weights:";
201 7044 : if( weights.size()!=atoms.size() ) error("number of elements in weight vector does not match the number of atoms");
202 35274 : for(unsigned i=0; i<weights.size(); ++i) {
203 28232 : if(i%25==0) log<<"\n";
204 28232 : log.printf(" %f",weights[i]);
205 : }
206 7042 : log.printf("\n");
207 : }
208 : }
209 9575 : if(phases) {
210 4 : log<<" Phases will be used to take into account PBC\n";
211 9571 : } else if(nopbc) {
212 47 : log<<" PBC will be ignored\n";
213 : } else {
214 9524 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
215 : }
216 9575 : requestAtoms(atoms);
217 9577 : }
218 :
219 32946 : void Center::calculate() {
220 32946 : Vector pos;
221 32946 : const bool dophases=(getPbc().isSet() ? phases : false);
222 :
223 32946 : if(!nopbc && !dophases) makeWhole();
224 :
225 32946 : if( first ) {
226 15376 : if( weight_mass ) {
227 172257 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
228 164230 : if(std::isnan(getMass(i))) {
229 0 : error(
230 : "You are trying to compute a CENTER or COM but masses are not known.\n"
231 : " If you are using plumed driver, please use the --mc option"
232 : );
233 : }
234 : }
235 : }
236 : double mass(0.0);
237 209364 : for(unsigned i=0; i<getNumberOfAtoms(); i++) mass+=getMass(i);
238 15376 : if( chargesWereSet && !isChargeSet_) {
239 : double charge(0.0);
240 192732 : for(unsigned i=0; i<getNumberOfAtoms(); i++) charge+=getCharge(i);
241 13318 : setCharge(charge);
242 2058 : } else if( !isChargeSet_ ) {
243 2056 : setCharge(0.0);
244 : }
245 15376 : if(!isMassSet_) setMass(mass);
246 :
247 15376 : if( weight_mass ) {
248 8027 : weights.resize( getNumberOfAtoms() );
249 172257 : for(unsigned i=0; i<weights.size(); i++) weights[i] = getMass(i) / mass;
250 : } else {
251 : double wtot=0.0;
252 37107 : for(unsigned i=0; i<weights.size(); i++) wtot+=weights[i];
253 37107 : for(unsigned i=0; i<weights.size(); i++) weights[i]=weights[i]/wtot;
254 7349 : first=false;
255 : }
256 : }
257 :
258 32946 : deriv.resize(getNumberOfAtoms());
259 :
260 32946 : if(dophases) {
261 241 : dcenter_sin.resize(getNumberOfAtoms());
262 241 : dcenter_cos.resize(getNumberOfAtoms());
263 241 : Vector center_sin;
264 241 : Vector center_cos;
265 241 : Tensor invbox2pi=2*pi*getPbc().getInvBox();
266 241 : Tensor box2pi=getPbc().getBox() / (2*pi);
267 964 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
268 723 : double w=weights[i];
269 :
270 : // real to scaled
271 723 : const Vector scaled=matmul(getPosition(i),invbox2pi);
272 : const Vector ccos(
273 723 : w*std::cos(scaled[0]),
274 723 : w*std::cos(scaled[1]),
275 723 : w*std::cos(scaled[2])
276 723 : );
277 : const Vector csin(
278 723 : w*std::sin(scaled[0]),
279 723 : w*std::sin(scaled[1]),
280 723 : w*std::sin(scaled[2])
281 723 : );
282 723 : center_cos+=ccos;
283 723 : center_sin+=csin;
284 9399 : for(unsigned l=0; l<3; l++) for(unsigned k=0; k<3; k++) {
285 : // k over real coordinates
286 : // l over scaled coordinates
287 6507 : dcenter_sin[i][l][k]=ccos[l]*invbox2pi[k][l];
288 6507 : dcenter_cos[i][l][k]=-csin[l]*invbox2pi[k][l];
289 : }
290 : }
291 : const Vector c(
292 241 : std::atan2(center_sin[0],center_cos[0]),
293 241 : std::atan2(center_sin[1],center_cos[1]),
294 241 : std::atan2(center_sin[2],center_cos[2])
295 241 : );
296 :
297 : // normalization is convenient for doing derivatives later
298 964 : for(unsigned l=0; l<3; l++) {
299 723 : double norm=1.0/(center_sin[l]*center_sin[l]+center_cos[l]*center_cos[l]);
300 723 : center_sin[l]*=norm;
301 723 : center_cos[l]*=norm;
302 : }
303 :
304 964 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
305 723 : Tensor dd;
306 9399 : for(unsigned l=0; l<3; l++) for(unsigned k=0; k<3; k++) {
307 : // k over real coordinates
308 : // l over scaled coordinates
309 6507 : dd[l][k]= (center_cos[l]*dcenter_sin[i][l][k] - center_sin[l]*dcenter_cos[i][l][k]);
310 : }
311 : // scaled to real
312 723 : deriv[i]=matmul(dd,box2pi);
313 : }
314 241 : setAtomsDerivatives(deriv);
315 : // scaled to real
316 241 : setPosition(matmul(c,box2pi));
317 : } else {
318 412212 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
319 379507 : double w=weights[i];
320 379507 : pos+=w*getPosition(i);
321 379507 : deriv[i]=w*Tensor::identity();
322 : }
323 32705 : setPosition(pos);
324 32705 : setAtomsDerivatives(deriv);
325 : }
326 32946 : }
327 :
328 : }
329 : }
|
