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 : {
122 : std::vector<double> weights;
123 : std::vector<Tensor> dcenter_sin;
124 : std::vector<Tensor> dcenter_cos;
125 : bool weight_mass;
126 : bool nopbc;
127 : bool first;
128 : bool phases;
129 : public:
130 : explicit Center(const ActionOptions&ao);
131 : void calculate() override;
132 : static void registerKeywords( Keywords& keys );
133 : };
134 :
135 24591 : PLUMED_REGISTER_ACTION(Center,"CENTER")
136 10567 : PLUMED_REGISTER_ACTION(Center,"COM")
137 :
138 7163 : void Center::registerKeywords(Keywords& keys) {
139 7163 : ActionWithVirtualAtom::registerKeywords(keys);
140 14326 : keys.add("optional","WEIGHTS","Center is computed as a weighted average.");
141 14326 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
142 14326 : keys.addFlag("MASS",false,"If set center is mass weighted");
143 14326 : keys.addFlag("PHASES",false,"Compute center using trigonometric phases");
144 7163 : }
145 :
146 7161 : Center::Center(const ActionOptions&ao):
147 : Action(ao),
148 : ActionWithVirtualAtom(ao),
149 7161 : weight_mass(false),
150 7161 : nopbc(false),
151 7161 : first(true),
152 7161 : phases(false)
153 : {
154 : std::vector<AtomNumber> atoms;
155 14322 : parseAtomList("ATOMS",atoms);
156 7161 : if(atoms.size()==0) error("at least one atom should be specified");
157 7161 : parseVector("WEIGHTS",weights);
158 7161 : parseFlag("MASS",weight_mass);
159 7161 : parseFlag("NOPBC",nopbc);
160 14322 : parseFlag("PHASES",phases);
161 7161 : if( getName()=="COM") weight_mass=true;
162 7161 : checkRead();
163 7161 : log.printf(" of atoms:");
164 37413 : for(unsigned i=0; i<atoms.size(); ++i) {
165 30252 : if(i%25==0) log<<"\n";
166 30252 : log.printf(" %d",atoms[i].serial());
167 : }
168 7161 : log<<"\n";
169 7161 : if(weight_mass) {
170 77 : log<<" mass weighted\n";
171 77 : if(weights.size()!=0) error("WEIGHTS and MASS keywords should not be used simultaneously");
172 : } else {
173 7084 : if( weights.size()==0) {
174 107 : log<<" using the geometric center\n";
175 107 : weights.resize( atoms.size() );
176 1318 : for(unsigned i=0; i<atoms.size(); i++) weights[i] = 1.;
177 : } else {
178 6977 : log<<" with weights:";
179 6979 : if( weights.size()!=atoms.size() ) error("number of elements in weight vector does not match the number of atoms");
180 35074 : for(unsigned i=0; i<weights.size(); ++i) {
181 28098 : if(i%25==0) log<<"\n";
182 28098 : log.printf(" %f",weights[i]);
183 : }
184 6976 : log.printf("\n");
185 : }
186 : }
187 7159 : if(phases) {
188 3 : log<<" Phases will be used to take into account PBC\n";
189 7156 : } else if(nopbc) {
190 45 : log<<" PBC will be ignored\n";
191 : } else {
192 7111 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
193 : }
194 7159 : requestAtoms(atoms);
195 7163 : }
196 :
197 13798 : void Center::calculate() {
198 13798 : Vector pos;
199 : double mass(0.0);
200 13798 : const bool dophases=(getPbc().isSet() ? phases : false);
201 :
202 13798 : if(!nopbc && !dophases) makeWhole();
203 :
204 13798 : if( first && weight_mass) {
205 725 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
206 659 : if(std::isnan(getMass(i))) {
207 0 : error(
208 : "You are trying to compute a CENTER or COM but masses are not known.\n"
209 : " If you are using plumed driver, please use the --mc option"
210 : );
211 : }
212 : }
213 66 : first=false;
214 : }
215 :
216 13798 : std::vector<Tensor> deriv(getNumberOfAtoms());
217 93241 : for(unsigned i=0; i<getNumberOfAtoms(); i++) mass+=getMass(i);
218 13798 : if( plumed.getAtoms().chargesWereSet() ) {
219 : double charge(0.0);
220 60621 : for(unsigned i=0; i<getNumberOfAtoms(); i++) charge+=getCharge(i);
221 : setCharge(charge);
222 : } else {
223 : setCharge(0.0);
224 : }
225 : double wtot=0.0;
226 57907 : for(unsigned i=0; i<weights.size(); i++) wtot+=weights[i];
227 :
228 13798 : if(dophases) {
229 240 : dcenter_sin.resize(getNumberOfAtoms());
230 240 : dcenter_cos.resize(getNumberOfAtoms());
231 240 : Vector center_sin;
232 240 : Vector center_cos;
233 240 : Tensor invbox2pi=2*pi*getPbc().getInvBox();
234 240 : Tensor box2pi=getPbc().getBox() / (2*pi);
235 960 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
236 : double w=0;
237 720 : if(weight_mass) w=getMass(i)/mass;
238 720 : else w=weights[i]/wtot;
239 :
240 : // real to scaled
241 720 : const Vector scaled=matmul(getPosition(i),invbox2pi);
242 : const Vector ccos(
243 720 : w*std::cos(scaled[0]),
244 720 : w*std::cos(scaled[1]),
245 720 : w*std::cos(scaled[2])
246 720 : );
247 : const Vector csin(
248 720 : w*std::sin(scaled[0]),
249 720 : w*std::sin(scaled[1]),
250 720 : w*std::sin(scaled[2])
251 720 : );
252 720 : center_cos+=ccos;
253 720 : center_sin+=csin;
254 9360 : for(unsigned l=0; l<3; l++) for(unsigned k=0; k<3; k++) {
255 : // k over real coordinates
256 : // l over scaled coordinates
257 6480 : dcenter_sin[i][l][k]=ccos[l]*invbox2pi[k][l];
258 6480 : dcenter_cos[i][l][k]=-csin[l]*invbox2pi[k][l];
259 : }
260 : }
261 : const Vector c(
262 240 : std::atan2(center_sin[0],center_cos[0]),
263 240 : std::atan2(center_sin[1],center_cos[1]),
264 240 : std::atan2(center_sin[2],center_cos[2])
265 240 : );
266 :
267 : // normalization is convenient for doing derivatives later
268 960 : for(unsigned l=0; l<3; l++) {
269 720 : double norm=1.0/(center_sin[l]*center_sin[l]+center_cos[l]*center_cos[l]);
270 720 : center_sin[l]*=norm;
271 720 : center_cos[l]*=norm;
272 : }
273 :
274 960 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
275 720 : Tensor dd;
276 9360 : for(unsigned l=0; l<3; l++) for(unsigned k=0; k<3; k++) {
277 : // k over real coordinates
278 : // l over scaled coordinates
279 6480 : dd[l][k]= (center_cos[l]*dcenter_sin[i][l][k] - center_sin[l]*dcenter_cos[i][l][k]);
280 : }
281 : // scaled to real
282 720 : deriv[i]=matmul(dd,box2pi);
283 : }
284 : setMass(mass);
285 : setAtomsDerivatives(deriv);
286 : // scaled to real
287 240 : setPosition(matmul(c,box2pi));
288 : } else {
289 92281 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
290 : double w=0;
291 78723 : if(weight_mass) w=getMass(i)/mass;
292 43389 : else w=weights[i]/wtot;
293 78723 : pos+=w*getPosition(i);
294 78723 : deriv[i]=w*Tensor::identity();
295 : }
296 : setPosition(pos);
297 : setMass(mass);
298 : setAtomsDerivatives(deriv);
299 : }
300 13798 : }
301 :
302 : }
303 : }
|
