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