Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2023 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 "Colvar.h"
23 : #include "ColvarShortcut.h"
24 : #include "MultiColvarTemplate.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/Pbc.h"
27 :
28 : namespace PLMD {
29 : namespace colvar {
30 :
31 : //+PLUMEDOC COLVAR POSITION
32 : /*
33 : Calculate the components of the position of an atom or atoms.
34 :
35 : To print the position of atom one to a file you can use an input like this:
36 :
37 : ```plumed
38 : p: POSITION ATOM=1
39 : PRINT ARG=p.x,p.y,p.z FILE=colvar
40 : ```
41 :
42 : To print the position of four atoms you can use an input like this:
43 :
44 : ```plumed
45 : p: POSITION ATOMS=1-4
46 : PRINT ARG=p.x,p.y,p.z FILE=colvar
47 : ```
48 :
49 : The three output values, p.x, p.y and p.z, here are all four dimensional vectors.
50 :
51 : [!CAUTION]
52 : Notice that single components will not have the proper periodicity!
53 :
54 : If you need the values to be consistent through PBC you can use SCALED_COMPONENTS,
55 : which defines values that by construction are in the -0.5,0.5 domain. This is
56 : similar to the equivalent flag for [DISTANCE](DISTANCE.md).
57 : Also notice that by default the minimal image distance from the
58 : origin is considered (can be changed with NOPBC).
59 :
60 : [!CAUTION]
61 : This variable should be used with extreme care since it allows you to easily get in troubles.
62 : It can be only be used if the
63 : Hamiltonian is not invariant for translation (i.e. there are other absolute positions which are biased, e.g. by position restraints)
64 : and cell size and shapes are fixed through the simulation.
65 :
66 : If you are not in this situation and still want to use the absolute position of an atom you should first fix the reference frame
67 : by using [FIT_TO_TEMPLATE](FIT_TO_TEMPLATE.md) as shown in the example below
68 :
69 : ```plumed
70 : #SETTINGS INPUTFILES=regtest/basic/rt63/align.pdb
71 : # align to a template
72 : FIT_TO_TEMPLATE REFERENCE=regtest/basic/rt63/align.pdb
73 : p: POSITION ATOM=3
74 : PRINT ARG=p.x,p.y,p.z
75 : ```
76 :
77 : */
78 : //+ENDPLUMEDOC
79 :
80 : //+PLUMEDOC COLVAR POSITION_SCALAR
81 : /*
82 : Calculate the components of the position of an atom.
83 :
84 : \par Examples
85 :
86 : */
87 : //+ENDPLUMEDOC
88 :
89 : //+PLUMEDOC COLVAR POSITION_VECTOR
90 : /*
91 : Create a vector that holds the components of the position of a set of atoms.
92 :
93 : \par Examples
94 :
95 : */
96 : //+ENDPLUMEDOC
97 :
98 : class Position : public Colvar {
99 : bool scaled_components;
100 : bool pbc;
101 : std::vector<double> value, masses, charges;
102 : std::vector<std::vector<Vector> > derivs;
103 : std::vector<Tensor> virial;
104 : public:
105 : static void registerKeywords( Keywords& keys );
106 : explicit Position(const ActionOptions&);
107 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
108 : static unsigned getModeAndSetupValues( ActionWithValue* av );
109 : // active methods:
110 : void calculate() override;
111 : static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
112 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
113 : std::vector<Tensor>& virial, const ActionAtomistic* aa );
114 : };
115 :
116 : typedef ColvarShortcut<Position> PositionShortcut;
117 : PLUMED_REGISTER_ACTION(PositionShortcut,"POSITION")
118 : PLUMED_REGISTER_ACTION(Position,"POSITION_SCALAR")
119 : typedef MultiColvarTemplate<Position> PositionMulti;
120 : PLUMED_REGISTER_ACTION(PositionMulti,"POSITION_VECTOR")
121 :
122 468 : void Position::registerKeywords( Keywords& keys ) {
123 468 : Colvar::registerKeywords( keys );
124 936 : keys.setDisplayName("POSITION");
125 468 : keys.add("atoms","ATOM","the atom number");
126 468 : keys.add("atoms","ATOMS","the atom numbers that you would like to use the positions of");
127 468 : keys.addFlag("WHOLEMOLECULES",false,"if this is a vector of positions do you want to make the positions into a whole before");
128 468 : keys.addFlag("SCALED_COMPONENTS",false,"calculate the a, b and c scaled components of the position separately and store them as label.a, label.b and label.c");
129 936 : keys.addOutputComponent("x","default","scalar/vector","the x-component of the atom position");
130 936 : keys.addOutputComponent("y","default","scalar/vector","the y-component of the atom position");
131 936 : keys.addOutputComponent("z","default","scalar/vector","the z-component of the atom position");
132 936 : keys.addOutputComponent("a","SCALED_COMPONENTS","scalar/vector","the normalized projection on the first lattice vector of the atom position");
133 936 : keys.addOutputComponent("b","SCALED_COMPONENTS","scalar/vector","the normalized projection on the second lattice vector of the atom position");
134 936 : keys.addOutputComponent("c","SCALED_COMPONENTS","scalar/vector","the normalized projection on the third lattice vector of the atom position");
135 468 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
136 468 : }
137 :
138 94 : Position::Position(const ActionOptions&ao):
139 : PLUMED_COLVAR_INIT(ao),
140 94 : scaled_components(false),
141 94 : pbc(true),
142 94 : value(3),
143 95 : derivs(3),
144 188 : virial(3) {
145 376 : for(unsigned i=0; i<3; ++i) {
146 282 : derivs[i].resize(1);
147 : }
148 : std::vector<AtomNumber> atoms;
149 94 : parseAtomList(-1,atoms,this);
150 93 : unsigned mode=getModeAndSetupValues(this);
151 93 : if( mode==1 ) {
152 7 : scaled_components=true;
153 : }
154 :
155 93 : bool nopbc=!pbc;
156 94 : parseFlag("NOPBC",nopbc);
157 93 : pbc=!nopbc;
158 93 : checkRead();
159 :
160 93 : if(pbc) {
161 88 : log.printf(" using periodic boundary conditions\n");
162 : } else {
163 5 : log.printf(" without periodic boundary conditions\n");
164 : }
165 :
166 93 : requestAtoms(atoms);
167 96 : }
168 :
169 102 : void Position::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
170 204 : aa->parseAtomList("ATOM",num,t);
171 102 : if( t.size()==1 ) {
172 99 : aa->log.printf(" for atom %d\n",t[0].serial());
173 3 : } else if( num<0 || t.size()!=0 ) {
174 1 : aa->error("Number of specified atoms should be 1");
175 : }
176 101 : }
177 :
178 134 : unsigned Position::getModeAndSetupValues( ActionWithValue* av ) {
179 : bool sc;
180 134 : av->parseFlag("SCALED_COMPONENTS",sc);
181 134 : if(sc) {
182 30 : av->addComponentWithDerivatives("a");
183 30 : av->componentIsPeriodic("a","-0.5","+0.5");
184 30 : av->addComponentWithDerivatives("b");
185 30 : av->componentIsPeriodic("b","-0.5","+0.5");
186 30 : av->addComponentWithDerivatives("c");
187 30 : av->componentIsPeriodic("c","-0.5","+0.5");
188 15 : return 1;
189 : }
190 238 : av->addComponentWithDerivatives("x");
191 119 : av->componentIsNotPeriodic("x");
192 238 : av->addComponentWithDerivatives("y");
193 119 : av->componentIsNotPeriodic("y");
194 238 : av->addComponentWithDerivatives("z");
195 119 : av->componentIsNotPeriodic("z");
196 119 : av->log<<" WARNING: components will not have the proper periodicity - see manual\n";
197 : return 0;
198 : }
199 :
200 : // calculator
201 8078 : void Position::calculate() {
202 :
203 8078 : std::vector<Vector> distance(1);
204 8078 : if(pbc) {
205 16044 : distance[0]=pbcDistance(Vector(0.0,0.0,0.0),getPosition(0));
206 : } else {
207 56 : distance[0]=delta(Vector(0.0,0.0,0.0),getPosition(0));
208 : }
209 :
210 8078 : if(scaled_components) {
211 56 : calculateCV( 1, masses, charges, distance, value, derivs, virial, this );
212 56 : Value* valuea=getPntrToComponent("a");
213 56 : Value* valueb=getPntrToComponent("b");
214 56 : Value* valuec=getPntrToComponent("c");
215 56 : setAtomsDerivatives (valuea,0,derivs[0][0]);
216 56 : valuea->set(value[0]);
217 56 : setAtomsDerivatives (valueb,0,derivs[1][0]);
218 56 : valueb->set(value[1]);
219 56 : setAtomsDerivatives (valuec,0,derivs[2][0]);
220 56 : valuec->set(value[2]);
221 : } else {
222 8022 : calculateCV( 0, masses, charges, distance, value, derivs, virial, this );
223 8022 : Value* valuex=getPntrToComponent("x");
224 8022 : Value* valuey=getPntrToComponent("y");
225 8022 : Value* valuez=getPntrToComponent("z");
226 :
227 8022 : setAtomsDerivatives (valuex,0,derivs[0][0]);
228 8022 : setBoxDerivatives (valuex,virial[0]);
229 8022 : valuex->set(value[0]);
230 :
231 8022 : setAtomsDerivatives (valuey,0,derivs[1][0]);
232 8022 : setBoxDerivatives (valuey,virial[1]);
233 8022 : valuey->set(value[1]);
234 :
235 8022 : setAtomsDerivatives (valuez,0,derivs[2][0]);
236 8022 : setBoxDerivatives (valuez,virial[2]);
237 8022 : valuez->set(value[2]);
238 : }
239 8078 : }
240 :
241 151422 : void Position::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
242 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
243 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
244 151422 : if( mode==1 ) {
245 10841 : Vector d=aa->getPbc().realToScaled(pos[0]);
246 10841 : vals[0]=Tools::pbc(d[0]);
247 10841 : vals[1]=Tools::pbc(d[1]);
248 10841 : vals[2]=Tools::pbc(d[2]);
249 10841 : derivs[0][0]=matmul(aa->getPbc().getInvBox(),Vector(+1,0,0));
250 10841 : derivs[1][0]=matmul(aa->getPbc().getInvBox(),Vector(0,+1,0));
251 10841 : derivs[2][0]=matmul(aa->getPbc().getInvBox(),Vector(0,0,+1));
252 : } else {
253 562324 : for(unsigned i=0; i<3; ++i) {
254 421743 : vals[i]=pos[0][i];
255 : }
256 140581 : derivs[0][0]=Vector(+1,0,0);
257 140581 : derivs[1][0]=Vector(0,+1,0);
258 140581 : derivs[2][0]=Vector(0,0,+1);
259 140581 : virial[0]=Tensor(pos[0],Vector(-1,0,0));
260 140581 : virial[1]=Tensor(pos[0],Vector(0,-1,0));
261 140581 : virial[2]=Tensor(pos[0],Vector(0,0,-1));
262 : }
263 151422 : }
264 :
265 : }
266 : }
267 :
268 :
269 :
|
