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 { return getNumberOfArguments(); }
54 : /// Do the calculation
55 : void calculate() override;
56 : ///
57 : void apply() override;
58 : };
59 :
60 : PLUMED_REGISTER_ACTION(ArgsToVatom,"ARGS2VATOM")
61 :
62 20 : void ArgsToVatom::registerKeywords( Keywords& keys ) {
63 20 : Action::registerKeywords( keys ); ActionWithValue::registerKeywords( keys ); ActionWithArguments::registerKeywords( keys );
64 40 : keys.add("compulsory","XPOS","the x position of the atom");
65 40 : keys.add("compulsory","YPOS","the y position of the atom");
66 40 : keys.add("compulsory","ZPOS","the z position of the atom");
67 40 : keys.add("compulsory","MASS","the mass of the atom");
68 40 : keys.add("compulsory","CHARGE","the charge of the atom");
69 40 : keys.add("hidden","XBKP","x position to use in case PBC not set when using PHASES");
70 40 : keys.add("hidden","YBKP","y position to use in case PBC not set when using PHASES");
71 40 : keys.add("hidden","ZBKP","z position to use in case PBC not set when using PHASES");
72 40 : keys.addFlag("FRACTIONAL",false,"the input arguments are calculated in fractional coordinates so you need to multiply by the cell");
73 40 : keys.addOutputComponent("x","default","the x coordinate of the virtual atom");
74 40 : keys.addOutputComponent("y","default","the y coordinate of the virtual atom");
75 40 : keys.addOutputComponent("z","default","the z coordinate of the virtual atom");
76 40 : keys.addOutputComponent("mass","default","the mass of the virtual atom");
77 40 : keys.addOutputComponent("charge","default","the charge of the virtual atom");
78 20 : }
79 :
80 9 : ArgsToVatom::ArgsToVatom(const ActionOptions& ao):
81 : Action(ao),
82 : ActionWithValue(ao),
83 9 : ActionWithArguments(ao)
84 : {
85 18 : parseFlag("FRACTIONAL",fractional);
86 18 : std::vector<Value*> xpos; parseArgumentList("XPOS",xpos);
87 9 : if( xpos.size()!=1 && xpos[0]->getRank()!=0 ) error("invalid input argument for XPOS");
88 18 : std::vector<Value*> ypos; parseArgumentList("YPOS",ypos);
89 9 : if( ypos.size()!=1 && ypos[0]->getRank()!=0 ) error("invalid input argument for YPOS");
90 18 : std::vector<Value*> zpos; parseArgumentList("ZPOS",zpos);
91 9 : if( zpos.size()!=1 && zpos[0]->getRank()!=0 ) error("invalid input argument for ZPOS");
92 18 : std::vector<Value*> mass; parseArgumentList("MASS",mass);
93 9 : if( mass.size()!=1 && mass[0]->getRank()!=0 ) error("invalid input argument for MASS");
94 18 : std::vector<Value*> charge; parseArgumentList("CHARGE",charge);
95 9 : if( charge.size()!=1 && charge[0]->getRank()!=0 ) error("invalid input argument for CHARGE");
96 : // Make sure we have requested everything that we need in xpos
97 9 : xpos.push_back(ypos[0]); xpos.push_back(zpos[0]); xpos.push_back(mass[0]); xpos.push_back(charge[0]);
98 9 : if( fractional ) {
99 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() );
100 14 : std::vector<Value*> xbkp; parseArgumentList("XBKP",xbkp);
101 7 : if( xbkp.size()>0 ) {
102 0 : if( xbkp.size()!=1 && xbkp[0]->getRank()!=0 ) error("invalid input argument for XBKP");
103 0 : std::vector<Value*> ybkp; parseArgumentList("YBKP",ybkp);
104 0 : if( ybkp.size()!=1 && ybkp[0]->getRank()!=0 ) error("invalid input argument for YBKP");
105 0 : std::vector<Value*> zbkp; parseArgumentList("ZBKP",zbkp);
106 0 : if( zbkp.size()!=1 && zpos[0]->getRank()!=0 ) error("invalid input argument for ZBKP");
107 : // Store backup for NOPBC
108 0 : xpos.push_back(xbkp[0]); xpos.push_back(ybkp[0]); xpos.push_back(zbkp[0]);
109 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() );
110 : }
111 2 : } else 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() );
112 9 : log.printf(" mass of atom is %s and charge is %s \n", mass[0]->getName().c_str(), charge[0]->getName().c_str() );
113 : // Request the arguments
114 9 : requestArguments(xpos);
115 : // Create the components to hold the atom
116 27 : addComponentWithDerivatives("x"); componentIsNotPeriodic("x");
117 27 : addComponentWithDerivatives("y"); componentIsNotPeriodic("y");
118 27 : addComponentWithDerivatives("z"); componentIsNotPeriodic("z");
119 27 : addComponent("mass"); componentIsNotPeriodic("mass"); if( mass[0]->isConstant() ) getPntrToComponent(3)->setConstant();
120 27 : addComponent("charge"); componentIsNotPeriodic("charge"); if( charge[0]->isConstant() ) getPntrToComponent(4)->setConstant();
121 9 : pbc_action = plumed.getActionSet().selectWithLabel<PbcAction*>("Box");
122 36 : for(unsigned i=0; i<3; ++i) getPntrToComponent(i)->resizeDerivatives( getNumberOfArguments() );
123 9 : }
124 :
125 26 : void ArgsToVatom::calculate() {
126 26 : if( fractional ) {
127 24 : if( pbc_action->getPbc().isSet() ) {
128 : // Get the position in fractional coordinates
129 96 : Vector fpos; for(unsigned i=0; i<3; ++i) fpos[i] = getPntrToArgument(i)->get();
130 : // Convert fractioanl coordinates to cartesian coordinates
131 24 : Tensor box=pbc_action->getPbc().getBox(); Vector cpos=matmul(fpos,box);
132 : // Set the final position and derivatives
133 96 : for(unsigned i=0; i<3; ++i) {
134 72 : Value* vv=getPntrToComponent(i); vv->set( cpos[i] );
135 288 : for(unsigned j=0; j<3; ++j) vv->addDerivative( j, box[j][i] );
136 : }
137 : } else {
138 0 : if( getNumberOfArguments()<8 ) error("cannot use PHASES option if box is not set");
139 : // Set the values
140 0 : for(unsigned i=0; i<3; ++i) getPntrToComponent(i)->set( getPntrToArgument(5+i)->get() );
141 : // And the derivatives
142 0 : for(unsigned i=0; i<3; ++i) getPntrToComponent(i)->addDerivative( 5+i, 1.0 );
143 : }
144 : // Set the mass and charge
145 72 : for(unsigned i=3; i<5; ++i) getPntrToComponent(i)->set( getPntrToArgument(i)->get() );
146 : } else {
147 : // Set the values
148 12 : for(unsigned i=0; i<5; ++i) getPntrToComponent(i)->set( getPntrToArgument(i)->get() );
149 : // And the derivatives
150 8 : for(unsigned i=0; i<3; ++i) getPntrToComponent(i)->addDerivative( i, 1.0 );
151 : }
152 26 : }
153 :
154 26 : void ArgsToVatom::apply() {
155 26 : if( !checkForForces() ) return ;
156 :
157 7 : unsigned start=0; addForcesOnArguments( 0, getForcesToApply(), start, getLabel() );
158 : }
159 :
160 : }
161 : }
|
