Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2020 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 : using namespace std;
29 :
30 : namespace PLMD {
31 : namespace vatom {
32 :
33 : //+PLUMEDOC VATOM CENTER
34 : /*
35 : Calculate the center for a group of atoms, with arbitrary weights.
36 :
37 : The computed
38 : center is stored as a virtual atom that can be accessed in
39 : an atom list through the label for the CENTER action that creates it.
40 : Notice that the generated virtual atom has charge equal to the sum of the
41 : charges and mass equal to the sum of the masses. If used with the MASS flag,
42 : then it provides a result identical to \ref COM.
43 :
44 : When running with periodic boundary conditions, the atoms should be
45 : in the proper periodic image. This is done automatically since PLUMED 2.2,
46 : by considering the ordered list of atoms and rebuilding the molecule using a procedure
47 : that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
48 : rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
49 : which actually modifies the coordinates stored in PLUMED.
50 :
51 : In case you want to recover the old behavior you should use the NOPBC flag.
52 : In that case you need to take care that atoms are in the correct
53 : periodic image.
54 :
55 : \note As an experimental feature, CENTER also supports a keyword PHASES.
56 : 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
57 : trigonometric functions, similarly to \ref CENTER_OF_MULTICOLVAR.
58 : Notice that by construction this center position is
59 : not invariant with respect to rotations of the atoms at fixed cell lattice.
60 : In addition, for symmetric Bravais lattices, it is not invariant with respect
61 : to special symmetries. E.g., if you have an hexagonal cell, the center will
62 : not be invariant with respect to rotations of 120 degrees.
63 : On the other hand, it might make the treatment of PBC easier in difficult cases.
64 :
65 : \par Examples
66 :
67 : \plumedfile
68 : # a point which is on the line connecting atoms 1 and 10, so that its distance
69 : # from 10 is twice its distance from 1:
70 : c1: CENTER ATOMS=1,1,10
71 : # this is another way of stating the same:
72 : c1bis: CENTER ATOMS=1,10 WEIGHTS=2,1
73 :
74 : # center of mass among these atoms:
75 : c2: CENTER ATOMS=2,3,4,5 MASS
76 :
77 : d1: DISTANCE ATOMS=c1,c2
78 :
79 : PRINT ARG=d1
80 : \endplumedfile
81 :
82 : */
83 : //+ENDPLUMEDOC
84 :
85 : //+PLUMEDOC VATOM COM
86 : /*
87 : Calculate the center of mass for a group of atoms.
88 :
89 : The computed
90 : center of mass is stored as a virtual atom that can be accessed in
91 : an atom list through the label for the COM action that creates it.
92 :
93 : For arbitrary weights (e.g. geometric center) see \ref CENTER.
94 :
95 : When running with periodic boundary conditions, the atoms should be
96 : in the proper periodic image. This is done automatically since PLUMED 2.2,
97 : by considering the ordered list of atoms and rebuilding the molecule using a procedure
98 : that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
99 : rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
100 : which actually modifies the coordinates stored in PLUMED.
101 :
102 : In case you want to recover the old behavior you should use the NOPBC flag.
103 : In that case you need to take care that atoms are in the correct
104 : periodic image.
105 :
106 : \par Examples
107 :
108 : The following input instructs plumed to print the distance between the
109 : center of mass for atoms 1,2,3,4,5,6,7 and that for atoms 15,20:
110 : \plumedfile
111 : c1: COM ATOMS=1-7
112 : c2: COM ATOMS=15,20
113 : d1: DISTANCE ATOMS=c1,c2
114 : PRINT ARG=d1
115 : \endplumedfile
116 :
117 : */
118 : //+ENDPLUMEDOC
119 :
120 :
121 528 : class Center:
122 : public ActionWithVirtualAtom
123 : {
124 : std::vector<double> weights;
125 : std::vector<Tensor> dcenter_sin;
126 : std::vector<Tensor> dcenter_cos;
127 : bool weight_mass;
128 : bool nopbc;
129 : bool first;
130 : bool phases;
131 : public:
132 : explicit Center(const ActionOptions&ao);
133 : void calculate();
134 : static void registerKeywords( Keywords& keys );
135 : };
136 :
137 7580 : PLUMED_REGISTER_ACTION(Center,"CENTER")
138 7486 : PLUMED_REGISTER_ACTION(Center,"COM")
139 :
140 180 : void Center::registerKeywords(Keywords& keys) {
141 180 : ActionWithVirtualAtom::registerKeywords(keys);
142 720 : keys.add("optional","WEIGHTS","Center is computed as a weighted average.");
143 540 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
144 540 : keys.addFlag("MASS",false,"If set center is mass weighted");
145 540 : keys.addFlag("PHASES",false,"Compute center using trigonometric phases");
146 180 : }
147 :
148 178 : Center::Center(const ActionOptions&ao):
149 : Action(ao),
150 : ActionWithVirtualAtom(ao),
151 : weight_mass(false),
152 : nopbc(false),
153 : first(true),
154 358 : phases(false)
155 : {
156 : vector<AtomNumber> atoms;
157 356 : parseAtomList("ATOMS",atoms);
158 178 : if(atoms.size()==0) error("at least one atom should be specified");
159 356 : parseVector("WEIGHTS",weights);
160 356 : parseFlag("MASS",weight_mass);
161 356 : parseFlag("NOPBC",nopbc);
162 356 : parseFlag("PHASES",phases);
163 178 : if( getName()=="COM") weight_mass=true;
164 178 : checkRead();
165 178 : log.printf(" of atoms:");
166 6635 : for(unsigned i=0; i<atoms.size(); ++i) {
167 2093 : if(i%25==0) log<<"\n";
168 4186 : log.printf(" %d",atoms[i].serial());
169 : }
170 178 : log<<"\n";
171 178 : if(weight_mass) {
172 68 : log<<" mass weighted\n";
173 69 : if(weights.size()!=0) error("WEIGHTS and MASS keywords should not be used simultaneously");
174 : } else {
175 110 : if( weights.size()==0) {
176 105 : log<<" using the geometric center\n";
177 105 : weights.resize( atoms.size() );
178 3801 : for(unsigned i=0; i<atoms.size(); i++) weights[i] = 1.;
179 : } else {
180 5 : log<<" with weights:";
181 6 : if( weights.size()!=atoms.size() ) error("number of elements in weight vector does not match the number of atoms");
182 86 : for(unsigned i=0; i<weights.size(); ++i) {
183 26 : if(i%25==0) log<<"\n";
184 52 : log.printf(" %f",weights[i]);
185 : }
186 4 : log.printf("\n");
187 : }
188 : }
189 176 : if(phases) {
190 3 : log<<" Phases will be used to take into account PBC\n";
191 173 : } else if(nopbc) {
192 45 : log<<" PBC will be ignored\n";
193 : } else {
194 128 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
195 : }
196 176 : requestAtoms(atoms);
197 176 : }
198 :
199 6599 : void Center::calculate() {
200 6599 : Vector pos;
201 : double mass(0.0);
202 6599 : const bool dophases=(getPbc().isSet() ? phases : false);
203 :
204 6599 : if(!nopbc && !dophases) makeWhole();
205 :
206 6599 : if( first && weight_mass) {
207 1236 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
208 589 : if(std::isnan(getMass(i))) {
209 0 : error(
210 : "You are trying to compute a CENTER or COM but masses are not known.\n"
211 : " If you are using plumed driver, please use the --mc option"
212 : );
213 : }
214 : }
215 58 : first=false;
216 : }
217 :
218 6599 : vector<Tensor> deriv(getNumberOfAtoms());
219 105429 : for(unsigned i=0; i<getNumberOfAtoms(); i++) mass+=getMass(i);
220 13198 : if( plumed.getAtoms().chargesWereSet() ) {
221 : double charge(0.0);
222 23399 : for(unsigned i=0; i<getNumberOfAtoms(); i++) charge+=getCharge(i);
223 : setCharge(charge);
224 : } else {
225 : setCharge(0.0);
226 : }
227 : double wtot=0.0;
228 61084 : for(unsigned i=0; i<weights.size(); i++) wtot+=weights[i];
229 :
230 6599 : if(dophases) {
231 456 : dcenter_sin.resize(getNumberOfAtoms());
232 456 : dcenter_cos.resize(getNumberOfAtoms());
233 228 : Vector center_sin;
234 228 : Vector center_cos;
235 228 : Tensor invbox2pi=2*pi*getPbc().getInvBox();
236 228 : Tensor box2pi=getPbc().getBox() / (2*pi);
237 1596 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
238 : double w=0;
239 684 : if(weight_mass) w=getMass(i)/mass;
240 1368 : else w=weights[i]/wtot;
241 :
242 : // real to scaled
243 1368 : const Vector scaled=matmul(getPosition(i),invbox2pi);
244 : const Vector ccos(
245 684 : w*std::cos(scaled[0]),
246 684 : w*std::cos(scaled[1]),
247 684 : w*std::cos(scaled[2])
248 2052 : );
249 : const Vector csin(
250 684 : w*std::sin(scaled[0]),
251 684 : w*std::sin(scaled[1]),
252 684 : w*std::sin(scaled[2])
253 2052 : );
254 684 : center_cos+=ccos;
255 684 : center_sin+=csin;
256 8892 : for(unsigned l=0; l<3; l++) for(unsigned k=0; k<3; k++) {
257 : // k over real coordinates
258 : // l over scaled coordinates
259 12312 : dcenter_sin[i][l][k]=ccos[l]*invbox2pi[k][l];
260 12312 : dcenter_cos[i][l][k]=-csin[l]*invbox2pi[k][l];
261 : }
262 : }
263 : const Vector c(
264 456 : std::atan2(center_sin[0],center_cos[0]),
265 456 : std::atan2(center_sin[1],center_cos[1]),
266 456 : std::atan2(center_sin[2],center_cos[2])
267 1368 : );
268 :
269 : // normalization is convenient for doing derivatives later
270 1596 : for(unsigned l=0; l<3; l++) {
271 684 : double norm=1.0/(center_sin[l]*center_sin[l]+center_cos[l]*center_cos[l]);
272 684 : center_sin[l]*=norm;
273 684 : center_cos[l]*=norm;
274 : }
275 :
276 1596 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
277 684 : Tensor dd;
278 8892 : for(unsigned l=0; l<3; l++) for(unsigned k=0; k<3; k++) {
279 : // k over real coordinates
280 : // l over scaled coordinates
281 18468 : dd[l][k]= (center_cos[l]*dcenter_sin[i][l][k] - center_sin[l]*dcenter_cos[i][l][k]);
282 : }
283 : // scaled to real
284 1368 : deriv[i]=matmul(dd,box2pi);
285 : }
286 : setMass(mass);
287 : setAtomsDerivatives(deriv);
288 : // scaled to real
289 456 : setPosition(matmul(c,box2pi));
290 : } else {
291 103833 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
292 : double w=0;
293 82184 : if(weight_mass) w=getMass(i)/mass;
294 30556 : else w=weights[i]/wtot;
295 97462 : pos+=w*getPosition(i);
296 97462 : deriv[i]=w*Tensor::identity();
297 : }
298 : setPosition(pos);
299 : setMass(mass);
300 : setAtomsDerivatives(deriv);
301 : }
302 6599 : }
303 :
304 : }
305 5517 : }
|
