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 : This input provides a very slow way of providing a center of mass - PLUMED contains a faster implementation that does all this.
65 : This type of input is nevertheless useful if you are using arbitary weights when computing the sums in the numerator
66 : and denominator of the expression for the center as is detailed in the documentation for the [CENTER](CENTER.md) command.
67 :
68 : */
69 : //+ENDPLUMEDOC
70 :
71 : namespace PLMD {
72 : namespace vatom {
73 :
74 : class ArgsToVatom :
75 : public ActionWithValue,
76 : public ActionWithArguments {
77 : private:
78 : bool fractional;
79 : PbcAction* pbc_action;
80 : public:
81 : static void registerKeywords( Keywords& keys );
82 : /// Constructor
83 : explicit ArgsToVatom(const ActionOptions&);
84 : /// Get the number of derivatives
85 121 : unsigned getNumberOfDerivatives() override {
86 121 : return getNumberOfArguments();
87 : }
88 : /// Do the calculation
89 : void calculate() override;
90 : ///
91 : void apply() override;
92 : };
93 :
94 : PLUMED_REGISTER_ACTION(ArgsToVatom,"ARGS2VATOM")
95 :
96 20 : void ArgsToVatom::registerKeywords( Keywords& keys ) {
97 20 : Action::registerKeywords( keys );
98 20 : ActionWithValue::registerKeywords( keys );
99 20 : ActionWithArguments::registerKeywords( keys );
100 40 : keys.addInputKeyword("compulsory","XPOS","scalar","the value to use for the x position of the atom");
101 40 : keys.addInputKeyword("compulsory","YPOS","scalar","the value to use for the y position of the atom");
102 40 : keys.addInputKeyword("compulsory","ZPOS","scalar","the value to use for the z position of the atom");
103 40 : keys.addInputKeyword("compulsory","MASS","scalar","the value to use for the mass of the atom");
104 40 : keys.addInputKeyword("compulsory","CHARGE","scalar","the value to use for the charge of the atom");
105 40 : keys.addInputKeyword("hidden","XBKP","scalar","x position to use in case PBC not set when using PHASES");
106 40 : keys.addInputKeyword("hidden","YBKP","scalar","y position to use in case PBC not set when using PHASES");
107 40 : keys.addInputKeyword("hidden","ZBKP","scalar","z position to use in case PBC not set when using PHASES");
108 20 : keys.addFlag("FRACTIONAL",false,"the input arguments are calculated in fractional coordinates so you need to multiply by the cell");
109 40 : keys.addOutputComponent("x","default","scalar","the x coordinate of the virtual atom");
110 40 : keys.addOutputComponent("y","default","scalar","the y coordinate of the virtual atom");
111 40 : keys.addOutputComponent("z","default","scalar","the z coordinate of the virtual atom");
112 40 : keys.addOutputComponent("mass","default","scalar","the mass of the virtual atom");
113 40 : keys.addOutputComponent("charge","default","scalar","the charge of the virtual atom");
114 20 : }
115 :
116 9 : ArgsToVatom::ArgsToVatom(const ActionOptions& ao):
117 : Action(ao),
118 : ActionWithValue(ao),
119 9 : ActionWithArguments(ao) {
120 18 : parseFlag("FRACTIONAL",fractional);
121 : std::vector<Value*> xpos;
122 18 : parseArgumentList("XPOS",xpos);
123 9 : if( xpos.size()!=1 && xpos[0]->getRank()!=0 ) {
124 0 : error("invalid input argument for XPOS");
125 : }
126 : std::vector<Value*> ypos;
127 18 : parseArgumentList("YPOS",ypos);
128 9 : if( ypos.size()!=1 && ypos[0]->getRank()!=0 ) {
129 0 : error("invalid input argument for YPOS");
130 : }
131 : std::vector<Value*> zpos;
132 18 : parseArgumentList("ZPOS",zpos);
133 9 : if( zpos.size()!=1 && zpos[0]->getRank()!=0 ) {
134 0 : error("invalid input argument for ZPOS");
135 : }
136 : std::vector<Value*> mass;
137 18 : parseArgumentList("MASS",mass);
138 9 : if( mass.size()!=1 && mass[0]->getRank()!=0 ) {
139 0 : error("invalid input argument for MASS");
140 : }
141 : std::vector<Value*> charge;
142 18 : parseArgumentList("CHARGE",charge);
143 9 : if( charge.size()!=1 && charge[0]->getRank()!=0 ) {
144 0 : error("invalid input argument for CHARGE");
145 : }
146 : // Make sure we have requested everything that we need in xpos
147 9 : xpos.push_back(ypos[0]);
148 9 : xpos.push_back(zpos[0]);
149 9 : xpos.push_back(mass[0]);
150 9 : xpos.push_back(charge[0]);
151 9 : if( fractional ) {
152 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() );
153 : std::vector<Value*> xbkp;
154 14 : parseArgumentList("XBKP",xbkp);
155 7 : if( xbkp.size()>0 ) {
156 0 : if( xbkp.size()!=1 && xbkp[0]->getRank()!=0 ) {
157 0 : error("invalid input argument for XBKP");
158 : }
159 : std::vector<Value*> ybkp;
160 0 : parseArgumentList("YBKP",ybkp);
161 0 : if( ybkp.size()!=1 && ybkp[0]->getRank()!=0 ) {
162 0 : error("invalid input argument for YBKP");
163 : }
164 : std::vector<Value*> zbkp;
165 0 : parseArgumentList("ZBKP",zbkp);
166 0 : if( zbkp.size()!=1 && zpos[0]->getRank()!=0 ) {
167 0 : error("invalid input argument for ZBKP");
168 : }
169 : // Store backup for NOPBC
170 0 : xpos.push_back(xbkp[0]);
171 0 : xpos.push_back(ybkp[0]);
172 0 : xpos.push_back(zbkp[0]);
173 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() );
174 : }
175 : } else {
176 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() );
177 : }
178 9 : log.printf(" mass of atom is %s and charge is %s \n", mass[0]->getName().c_str(), charge[0]->getName().c_str() );
179 : // Request the arguments
180 9 : requestArguments(xpos);
181 : // Create the components to hold the atom
182 18 : addComponentWithDerivatives("x");
183 18 : componentIsNotPeriodic("x");
184 18 : addComponentWithDerivatives("y");
185 18 : componentIsNotPeriodic("y");
186 18 : addComponentWithDerivatives("z");
187 18 : componentIsNotPeriodic("z");
188 18 : addComponent("mass");
189 18 : componentIsNotPeriodic("mass");
190 9 : if( mass[0]->isConstant() ) {
191 9 : getPntrToComponent(3)->setConstant();
192 : }
193 18 : addComponent("charge");
194 18 : componentIsNotPeriodic("charge");
195 9 : if( charge[0]->isConstant() ) {
196 9 : getPntrToComponent(4)->setConstant();
197 : }
198 9 : pbc_action = plumed.getActionSet().selectWithLabel<PbcAction*>("Box");
199 36 : for(unsigned i=0; i<3; ++i) {
200 27 : getPntrToComponent(i)->resizeDerivatives( getNumberOfArguments() );
201 : }
202 9 : }
203 :
204 26 : void ArgsToVatom::calculate() {
205 26 : if( fractional ) {
206 24 : if( pbc_action->getPbc().isSet() ) {
207 : // Get the position in fractional coordinates
208 24 : Vector fpos;
209 96 : for(unsigned i=0; i<3; ++i) {
210 72 : fpos[i] = getPntrToArgument(i)->get();
211 : }
212 : // Convert fractioanl coordinates to cartesian coordinates
213 24 : Tensor box=pbc_action->getPbc().getBox();
214 24 : Vector cpos=matmul(fpos,box);
215 : // Set the final position and derivatives
216 96 : for(unsigned i=0; i<3; ++i) {
217 72 : Value* vv=getPntrToComponent(i);
218 72 : vv->set( cpos[i] );
219 288 : for(unsigned j=0; j<3; ++j) {
220 216 : vv->addDerivative( j, box[j][i] );
221 : }
222 : }
223 : } else {
224 0 : if( getNumberOfArguments()<8 ) {
225 0 : error("cannot use PHASES option if box is not set");
226 : }
227 : // Set the values
228 0 : for(unsigned i=0; i<3; ++i) {
229 0 : getPntrToComponent(i)->set( getPntrToArgument(5+i)->get() );
230 : }
231 : // And the derivatives
232 0 : for(unsigned i=0; i<3; ++i) {
233 0 : getPntrToComponent(i)->addDerivative( 5+i, 1.0 );
234 : }
235 : }
236 : // Set the mass and charge
237 72 : for(unsigned i=3; i<5; ++i) {
238 48 : getPntrToComponent(i)->set( getPntrToArgument(i)->get() );
239 : }
240 : } else {
241 : // Set the values
242 12 : for(unsigned i=0; i<5; ++i) {
243 10 : getPntrToComponent(i)->set( getPntrToArgument(i)->get() );
244 : }
245 : // And the derivatives
246 8 : for(unsigned i=0; i<3; ++i) {
247 6 : getPntrToComponent(i)->addDerivative( i, 1.0 );
248 : }
249 : }
250 26 : }
251 :
252 26 : void ArgsToVatom::apply() {
253 26 : if( !checkForForces() ) {
254 19 : return ;
255 : }
256 :
257 7 : unsigned start=0;
258 7 : addForcesOnArguments( 0, getForcesToApply(), start, getLabel() );
259 : }
260 :
261 : }
262 : }
|
