Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2017 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 "core/ActionWithArguments.h"
23 : #include "core/ActionWithValue.h"
24 : #include "core/ActionRegister.h"
25 : #include "core/ActionSet.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/PbcAction.h"
28 : #include "tools/Pbc.h"
29 :
30 : //+PLUMEDOC VATOM ARGS2VATOM
31 : /*
32 : Create a virtual atom from the input scalars
33 :
34 : This action takes five scalars that are computed by other actions in input and uses them to set the
35 : x, y and z positions and the mass and charge of a virtual atom. This action is used within the
36 : [CENTER](CENTER.md) shortcut to compute a center of mass. An example input that shows how you
37 : can use this command to calculate the center of mass of atoms 1-10 is as follows:
38 :
39 : ```plumed
40 : # Calculate the total mass of the atoms
41 : m: MASS ATOMS=1-10
42 : mass: SUM ARG=m PERIODIC=NO
43 : # Calculate the totla charge of the atoms
44 : q: CHARGE ATOMS=1-10
45 : charge: SUM ARG=q PERIODIC=NO
46 : # Now get the positions of the atoms
47 : pos: POSITION WHOLEMOLECULES ATOMS=1-10
48 : # Multiply each vector of positions by the masses
49 : xwvec: CUSTOM ARG=m,pos.x FUNC=x*y PERIODIC=NO
50 : ywvec: CUSTOM ARG=m,pos.y FUNC=x*y PERIODIC=NO
51 : zwvec: CUSTOM ARG=m,pos.z FUNC=x*y PERIODIC=NO
52 : # Sum the numerators in the expression for the center of mass
53 : xnum: SUM ARG=xwvec PERIODIC=NO
54 : ynum: SUM ARG=ywvec PERIODIC=NO
55 : znum: SUM ARG=zwvec PERIODIC=NO
56 : # And compute the x, y and z positions of the center of mass
57 : x: CUSTOM ARG=xnum,mass FUNC=x/y PERIODIC=NO
58 : y: CUSTOM ARG=ynum,mass FUNC=x/y PERIODIC=NO
59 : z: CUSTOM ARG=znum,mass FUNC=x/y PERIODIC=NO
60 : # And now create the virtual atom
61 : p: ARGS2VATOM XPOS=x YPOS=y ZPOS=z MASS=mass CHARGE=charge
62 : ```
63 :
64 : In the following input by contrast we use the PHASES method that is discussed in the documentation for [CENTER](CENTER.md) to calculate the position of the center of mass:
65 :
66 : ```plumed
67 : # Calculate the total mass of the atoms
68 : m: MASS ATOMS=1-10
69 : mass: SUM ARG=m PERIODIC=NO
70 : # Calculate the totla charge of the atoms
71 : q: CHARGE ATOMS=1-10
72 : charge: SUM ARG=q PERIODIC=NO
73 : # Now get the positions of the atoms
74 : pos: POSITION SCALED_COMPONENTS ATOMS=1-10
75 : # Multiply the sins and cosines of the scaled positions by the weights
76 : sina: CUSTOM ARG=m,pos.a FUNC=x*sin(2*pi*y) PERIODIC=NO
77 : cosa: CUSTOM ARG=m,pos.a FUNC=x*cos(2*pi*y) PERIODIC=NO
78 : sinb: CUSTOM ARG=m,pos.b FUNC=x*sin(2*pi*y) PERIODIC=NO
79 : cosb: CUSTOM ARG=m,pos.b FUNC=x*cos(2*pi*y) PERIODIC=NO
80 : sinc: CUSTOM ARG=m,pos.c FUNC=x*sin(2*pi*y) PERIODIC=NO
81 : cosc: CUSTOM ARG=m,pos.c FUNC=x*cos(2*pi*y) PERIODIC=NO
82 : # And accumulate the sums
83 : sinsuma: SUM ARG=sina PERIODIC=NO
84 : cossuma: SUM ARG=cosa PERIODIC=NO
85 : sinsumb: SUM ARG=sinb PERIODIC=NO
86 : cossumb: SUM ARG=cosb PERIODIC=NO
87 : sinsumc: SUM ARG=sinc PERIODIC=NO
88 : cossumc: SUM ARG=cosc PERIODIC=NO
89 : # And get the position of the center in fractional coordinates
90 : a: CUSTOM ARG=sinsuma,cossuma FUNC=atan2(x,y)/(2*pi) PERIODIC=NO
91 : b: CUSTOM ARG=sinsumb,cossumb FUNC=atan2(x,y)/(2*pi) PERIODIC=NO
92 : c: CUSTOM ARG=sinsumc,cossumc FUNC=atan2(x,y)/(2*pi) PERIODIC=NO
93 : # And now create the virtual atom
94 : p: ARGS2VATOM XPOS=a YPOS=b ZPOS=c MASS=mass CHARGE=charge FRACTIONAL
95 : ```
96 :
97 : These inputs provide very slow ways of computing a center of mass - PLUMED contains faster implementations that do calculations that are equivalent to both of these inputs.
98 : This type of input is nevertheless useful if you are using arbitary weights when computing the sums in the numerator
99 : and denominator of the expression for the center as is detailed in the documentation for the [CENTER](CENTER.md) command.
100 :
101 : */
102 : //+ENDPLUMEDOC
103 :
104 : namespace PLMD {
105 : namespace vatom {
106 :
107 : class ArgsToVatom :
108 : public ActionWithValue,
109 : public ActionWithArguments {
110 : private:
111 : bool fractional;
112 : PbcAction* pbc_action;
113 : public:
114 : static void registerKeywords( Keywords& keys );
115 : /// Constructor
116 : explicit ArgsToVatom(const ActionOptions&);
117 : /// Get the number of derivatives
118 116 : unsigned getNumberOfDerivatives() override {
119 116 : return getNumberOfArguments();
120 : }
121 : /// Do the calculation
122 : void calculate() override;
123 : ///
124 : void apply() override;
125 : };
126 :
127 : PLUMED_REGISTER_ACTION(ArgsToVatom,"ARGS2VATOM")
128 :
129 22 : void ArgsToVatom::registerKeywords( Keywords& keys ) {
130 22 : Action::registerKeywords( keys );
131 22 : ActionWithValue::registerKeywords( keys );
132 22 : ActionWithArguments::registerKeywords( keys );
133 22 : keys.remove("NUMERICAL_DERIVATIVES");
134 44 : keys.addInputKeyword("compulsory","XPOS","scalar","the value to use for the x position of the atom");
135 44 : keys.addInputKeyword("compulsory","YPOS","scalar","the value to use for the y position of the atom");
136 44 : keys.addInputKeyword("compulsory","ZPOS","scalar","the value to use for the z position of the atom");
137 44 : keys.addInputKeyword("compulsory","MASS","scalar","the value to use for the mass of the atom");
138 44 : keys.addInputKeyword("compulsory","CHARGE","scalar","the value to use for the charge of the atom");
139 44 : keys.addInputKeyword("hidden","XBKP","scalar","x position to use in case PBC not set when using PHASES");
140 44 : keys.addInputKeyword("hidden","YBKP","scalar","y position to use in case PBC not set when using PHASES");
141 44 : keys.addInputKeyword("hidden","ZBKP","scalar","z position to use in case PBC not set when using PHASES");
142 22 : keys.addFlag("FRACTIONAL",false,"the input arguments are calculated in fractional coordinates so you need to multiply by the cell");
143 44 : keys.addOutputComponent("x","default","scalar","the x coordinate of the virtual atom");
144 44 : keys.addOutputComponent("y","default","scalar","the y coordinate of the virtual atom");
145 44 : keys.addOutputComponent("z","default","scalar","the z coordinate of the virtual atom");
146 44 : keys.addOutputComponent("mass","default","scalar","the mass of the virtual atom");
147 44 : keys.addOutputComponent("charge","default","scalar","the charge of the virtual atom");
148 22 : }
149 :
150 9 : ArgsToVatom::ArgsToVatom(const ActionOptions& ao):
151 : Action(ao),
152 : ActionWithValue(ao),
153 9 : ActionWithArguments(ao) {
154 9 : parseFlag("FRACTIONAL",fractional);
155 9 : std::vector<Value*> xpos;
156 18 : parseArgumentList("XPOS",xpos);
157 9 : if( xpos.size()!=1 || xpos[0]->getRank()!=0 ) {
158 0 : error("invalid input argument for XPOS");
159 : }
160 9 : std::vector<Value*> ypos;
161 18 : parseArgumentList("YPOS",ypos);
162 9 : if( ypos.size()!=1 || ypos[0]->getRank()!=0 ) {
163 0 : error("invalid input argument for YPOS");
164 : }
165 9 : std::vector<Value*> zpos;
166 18 : parseArgumentList("ZPOS",zpos);
167 9 : if( zpos.size()!=1 || zpos[0]->getRank()!=0 ) {
168 0 : error("invalid input argument for ZPOS");
169 : }
170 9 : std::vector<Value*> mass;
171 18 : parseArgumentList("MASS",mass);
172 9 : if( mass.size()!=1 || mass[0]->getRank()!=0 ) {
173 0 : error("invalid input argument for MASS");
174 : }
175 9 : std::vector<Value*> charge;
176 18 : parseArgumentList("CHARGE",charge);
177 9 : if( charge.size()!=1 || charge[0]->getRank()!=0 ) {
178 0 : error("invalid input argument for CHARGE");
179 : }
180 : // Make sure we have requested everything that we need in xpos
181 9 : xpos.push_back(ypos[0]);
182 9 : xpos.push_back(zpos[0]);
183 9 : xpos.push_back(mass[0]);
184 9 : xpos.push_back(charge[0]);
185 9 : if( fractional ) {
186 7 : log.printf(" creating atom from fractional pos a=%s, b=%s and c=%s \n", xpos[0]->getName().c_str(), ypos[0]->getName().c_str(), zpos[0]->getName().c_str() );
187 7 : std::vector<Value*> xbkp;
188 14 : parseArgumentList("XBKP",xbkp);
189 7 : if( xbkp.size()>0 ) {
190 0 : if( xbkp.size()!=1 || xbkp[0]->getRank()!=0 ) {
191 0 : error("invalid input argument for XBKP");
192 : }
193 0 : std::vector<Value*> ybkp;
194 0 : parseArgumentList("YBKP",ybkp);
195 0 : if( ybkp.size()!=1 || ybkp[0]->getRank()!=0 ) {
196 0 : error("invalid input argument for YBKP");
197 : }
198 0 : std::vector<Value*> zbkp;
199 0 : parseArgumentList("ZBKP",zbkp);
200 0 : if( zbkp.size()!=1 || zbkp[0]->getRank()!=0 ) {
201 0 : error("invalid input argument for ZBKP");
202 : }
203 : // Store backup for NOPBC
204 0 : xpos.push_back(xbkp[0]);
205 0 : xpos.push_back(ybkp[0]);
206 0 : xpos.push_back(zbkp[0]);
207 0 : log.printf(" using x=%s, y=%s and z=%s if PBC not set \n", xbkp[0]->getName().c_str(), ybkp[0]->getName().c_str(), zbkp[0]->getName().c_str() );
208 : }
209 : } else {
210 2 : log.printf(" creating atom at x=%s, y=%s and z=%s \n", xpos[0]->getName().c_str(), ypos[0]->getName().c_str(), zpos[0]->getName().c_str() );
211 : }
212 9 : log.printf(" mass of atom is %s and charge is %s \n", mass[0]->getName().c_str(), charge[0]->getName().c_str() );
213 : // Request the arguments
214 9 : requestArguments(xpos);
215 : // Create the components to hold the atom
216 18 : addComponentWithDerivatives("x");
217 18 : componentIsNotPeriodic("x");
218 18 : addComponentWithDerivatives("y");
219 18 : componentIsNotPeriodic("y");
220 18 : addComponentWithDerivatives("z");
221 18 : componentIsNotPeriodic("z");
222 18 : addComponent("mass");
223 18 : componentIsNotPeriodic("mass");
224 9 : if( mass[0]->isConstant() ) {
225 9 : getPntrToComponent(3)->setConstant();
226 : }
227 18 : addComponent("charge");
228 18 : componentIsNotPeriodic("charge");
229 9 : if( charge[0]->isConstant() ) {
230 9 : getPntrToComponent(4)->setConstant();
231 : }
232 9 : pbc_action = plumed.getActionSet().selectWithLabel<PbcAction*>("Box");
233 36 : for(unsigned i=0; i<3; ++i) {
234 27 : getPntrToComponent(i)->resizeDerivatives( getNumberOfArguments() );
235 : }
236 9 : }
237 :
238 26 : void ArgsToVatom::calculate() {
239 26 : if( fractional ) {
240 24 : if( pbc_action->getPbc().isSet() ) {
241 : // Get the position in fractional coordinates
242 24 : Vector fpos;
243 96 : for(unsigned i=0; i<3; ++i) {
244 72 : fpos[i] = getPntrToArgument(i)->get();
245 : }
246 : // Convert fractioanl coordinates to cartesian coordinates
247 24 : Tensor box=pbc_action->getPbc().getBox();
248 24 : Vector cpos=matmul(fpos,box);
249 : // Set the final position and derivatives
250 96 : for(unsigned i=0; i<3; ++i) {
251 72 : Value* vv=getPntrToComponent(i);
252 72 : vv->set( cpos[i] );
253 288 : for(unsigned j=0; j<3; ++j) {
254 216 : vv->addDerivative( j, box[j][i] );
255 : }
256 : }
257 : } else {
258 0 : if( getNumberOfArguments()<8 ) {
259 0 : error("cannot use PHASES option if box is not set");
260 : }
261 : // Set the values
262 0 : for(unsigned i=0; i<3; ++i) {
263 0 : getPntrToComponent(i)->set( getPntrToArgument(5+i)->get() );
264 : }
265 : // And the derivatives
266 0 : for(unsigned i=0; i<3; ++i) {
267 0 : getPntrToComponent(i)->addDerivative( 5+i, 1.0 );
268 : }
269 : }
270 : // Set the mass and charge
271 72 : for(unsigned i=3; i<5; ++i) {
272 48 : getPntrToComponent(i)->set( getPntrToArgument(i)->get() );
273 : }
274 : } else {
275 : // Set the values
276 12 : for(unsigned i=0; i<5; ++i) {
277 10 : getPntrToComponent(i)->set( getPntrToArgument(i)->get() );
278 : }
279 : // And the derivatives
280 8 : for(unsigned i=0; i<3; ++i) {
281 6 : getPntrToComponent(i)->addDerivative( i, 1.0 );
282 : }
283 : }
284 26 : }
285 :
286 26 : void ArgsToVatom::apply() {
287 26 : if( !checkForForces() ) {
288 19 : return ;
289 : }
290 :
291 7 : unsigned start=0;
292 7 : addForcesOnArguments( 0, getForcesToApply(), start );
293 : }
294 :
295 : }
296 : }
|
