Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2019 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 :
27 : using namespace std;
28 :
29 : namespace PLMD {
30 : namespace vatom {
31 :
32 : //+PLUMEDOC VATOM CENTER
33 : /*
34 : Calculate the center for a group of atoms, with arbitrary weights.
35 :
36 : The computed
37 : center is stored as a virtual atom that can be accessed in
38 : an atom list through the label for the CENTER action that creates it.
39 : Notice that the generated virtual atom has charge equal to the sum of the
40 : charges and mass equal to the sum of the masses. If used with the MASS flag,
41 : then it provides a result identical to \ref COM.
42 :
43 : When running with periodic boundary conditions, the atoms should be
44 : in the proper periodic image. This is done automatically since PLUMED 2.2,
45 : by considering the ordered list of atoms and rebuilding PBCs with a procedure
46 : that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
47 : rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
48 : which actually modifies the coordinates stored in PLUMED.
49 :
50 : In case you want to recover the old behavior you should use the NOPBC flag.
51 : In that case you need to take care that atoms are in the correct
52 : periodic image.
53 :
54 :
55 : \par Examples
56 :
57 : \plumedfile
58 : # a point which is on the line connecting atoms 1 and 10, so that its distance
59 : # from 10 is twice its distance from 1:
60 : c1: CENTER ATOMS=1,1,10
61 : # this is another way of stating the same:
62 : c1bis: CENTER ATOMS=1,10 WEIGHTS=2,1
63 :
64 : # center of mass among these atoms:
65 : c2: CENTER ATOMS=2,3,4,5 MASS
66 :
67 : d1: DISTANCE ATOMS=c1,c2
68 :
69 : PRINT ARG=d1
70 : \endplumedfile
71 :
72 : */
73 : //+ENDPLUMEDOC
74 :
75 :
76 219 : class Center:
77 : public ActionWithVirtualAtom
78 : {
79 : std::vector<double> weights;
80 : bool weight_mass;
81 : bool nopbc;
82 : public:
83 : explicit Center(const ActionOptions&ao);
84 : void calculate();
85 : static void registerKeywords( Keywords& keys );
86 : };
87 :
88 6527 : PLUMED_REGISTER_ACTION(Center,"CENTER")
89 :
90 76 : void Center::registerKeywords(Keywords& keys) {
91 76 : ActionWithVirtualAtom::registerKeywords(keys);
92 304 : keys.add("optional","WEIGHTS","Center is computed as a weighted average.");
93 228 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
94 228 : keys.addFlag("MASS",false,"If set center is mass weighted");
95 76 : }
96 :
97 75 : Center::Center(const ActionOptions&ao):
98 : Action(ao),
99 : ActionWithVirtualAtom(ao),
100 : weight_mass(false),
101 152 : nopbc(false)
102 : {
103 : vector<AtomNumber> atoms;
104 150 : parseAtomList("ATOMS",atoms);
105 75 : if(atoms.size()==0) error("at least one atom should be specified");
106 150 : parseVector("WEIGHTS",weights);
107 150 : parseFlag("MASS",weight_mass);
108 150 : parseFlag("NOPBC",nopbc);
109 75 : checkRead();
110 75 : log.printf(" of atoms:");
111 2862 : for(unsigned i=0; i<atoms.size(); ++i) {
112 904 : if(i%25==0) log<<"\n";
113 1808 : log.printf(" %d",atoms[i].serial());
114 : }
115 75 : log<<"\n";
116 75 : if(weight_mass) {
117 3 : log<<" mass weighted\n";
118 4 : if(weights.size()!=0) error("WEIGHTS and MASS keywords should not be used simultaneously");
119 : } else {
120 72 : if( weights.size()==0) {
121 68 : log<<" using the geometric center\n";
122 68 : weights.resize( atoms.size() );
123 2695 : for(unsigned i=0; i<atoms.size(); i++) weights[i] = 1.;
124 : } else {
125 4 : log<<" with weights:";
126 5 : if( weights.size()!=atoms.size() ) error("number of elements in weight vector does not match the number of atoms");
127 75 : for(unsigned i=0; i<weights.size(); ++i) {
128 23 : if(i%25==0) log<<"\n";
129 46 : log.printf(" %f",weights[i]);
130 : }
131 3 : log.printf("\n");
132 : }
133 : }
134 73 : if(nopbc) {
135 17 : log<<" PBC will be ignored\n";
136 : } else {
137 56 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
138 : }
139 73 : requestAtoms(atoms);
140 73 : }
141 :
142 940 : void Center::calculate() {
143 940 : Vector pos;
144 : double mass(0.0);
145 940 : if(!nopbc) makeWhole();
146 940 : vector<Tensor> deriv(getNumberOfAtoms());
147 16874 : for(unsigned i=0; i<getNumberOfAtoms(); i++) mass+=getMass(i);
148 1880 : if( plumed.getAtoms().chargesWereSet() ) {
149 : double charge(0.0);
150 8907 : for(unsigned i=0; i<getNumberOfAtoms(); i++) charge+=getCharge(i);
151 : setCharge(charge);
152 : } else {
153 : setCharge(0.0);
154 : }
155 : double wtot=0.0;
156 25481 : for(unsigned i=0; i<weights.size(); i++) wtot+=weights[i];
157 16874 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
158 : double w=0;
159 8067 : if(weight_mass) w=getMass(i)/mass;
160 15734 : else w=weights[i]/wtot;
161 15934 : pos+=w*getPosition(i);
162 15934 : deriv[i]=w*Tensor::identity();
163 : }
164 : setPosition(pos);
165 : setMass(mass);
166 : setAtomsDerivatives(deriv);
167 940 : }
168 :
169 : }
170 4839 : }
|
