Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2018 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 : \verbatim
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 : \endverbatim
71 : (See also \ref DISTANCE, \ref COM and \ref PRINT).
72 :
73 : */
74 : //+ENDPLUMEDOC
75 :
76 :
77 136 : class Center:
78 : public ActionWithVirtualAtom
79 : {
80 : std::vector<double> weights;
81 : bool weight_mass;
82 : bool nopbc;
83 : public:
84 : explicit Center(const ActionOptions&ao);
85 : void calculate();
86 : static void registerKeywords( Keywords& keys );
87 : };
88 :
89 2591 : PLUMED_REGISTER_ACTION(Center,"CENTER")
90 :
91 69 : void Center::registerKeywords(Keywords& keys) {
92 69 : ActionWithVirtualAtom::registerKeywords(keys);
93 69 : keys.add("optional","WEIGHTS","Center is computed as a weighted average.");
94 69 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
95 69 : keys.addFlag("MASS",false,"If set center is mass weighted");
96 69 : }
97 :
98 68 : Center::Center(const ActionOptions&ao):
99 : Action(ao),
100 : ActionWithVirtualAtom(ao),
101 : weight_mass(false),
102 68 : nopbc(false)
103 : {
104 68 : vector<AtomNumber> atoms;
105 68 : parseAtomList("ATOMS",atoms);
106 68 : if(atoms.size()==0) error("at least one atom should be specified");
107 68 : parseVector("WEIGHTS",weights);
108 68 : parseFlag("MASS",weight_mass);
109 68 : parseFlag("NOPBC",nopbc);
110 68 : checkRead();
111 68 : log.printf(" of atoms");
112 68 : for(unsigned i=0; i<atoms.size(); ++i) log.printf(" %d",atoms[i].serial());
113 68 : if(weight_mass) {
114 2 : log<<" mass weighted\n";
115 2 : if(weights.size()!=0) error("WEIGHTS and MASS keywords should not be used simultaneously");
116 : } else {
117 66 : if( weights.size()==0) {
118 63 : weights.resize( atoms.size() );
119 63 : for(unsigned i=0; i<atoms.size(); i++) weights[i] = 1.;
120 : }
121 66 : log<<" with weights";
122 66 : if( weights.size()!=atoms.size() ) error("number of elements in weight vector does not match the number of atoms");
123 66 : for(unsigned i=0; i<weights.size(); ++i) log.printf(" %f",weights[i]);
124 66 : log.printf("\n");
125 : }
126 68 : if(nopbc) {
127 17 : log<<" PBC will be ignored\n";
128 : } else {
129 51 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
130 : }
131 68 : requestAtoms(atoms);
132 68 : }
133 :
134 931 : void Center::calculate() {
135 931 : Vector pos;
136 931 : double mass(0.0);
137 931 : if(!nopbc) makeWhole();
138 931 : vector<Tensor> deriv(getNumberOfAtoms());
139 931 : for(unsigned i=0; i<getNumberOfAtoms(); i++) mass+=getMass(i);
140 931 : if( plumed.getAtoms().chargesWereSet() ) {
141 931 : double charge(0.0);
142 931 : for(unsigned i=0; i<getNumberOfAtoms(); i++) charge+=getCharge(i);
143 931 : setCharge(charge);
144 : } else {
145 0 : setCharge(0.0);
146 : }
147 931 : double wtot=0.0;
148 931 : for(unsigned i=0; i<weights.size(); i++) wtot+=weights[i];
149 8806 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
150 7875 : double w=0;
151 7875 : if(weight_mass) w=getMass(i)/mass;
152 7775 : else w=weights[i]/wtot;
153 7875 : pos+=w*getPosition(i);
154 7875 : deriv[i]=w*Tensor::identity();
155 : }
156 931 : setPosition(pos);
157 931 : setMass(mass);
158 931 : setAtomsDerivatives(deriv);
159 931 : }
160 :
161 : }
162 2523 : }
|
