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