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 ); keys.setDisplayName("POSITION");
119 936 : keys.add("atoms","ATOM","the atom number");
120 936 : keys.add("atoms","ATOMS","the atom numbers that you would like to use the positions of");
121 936 : keys.addFlag("WHOLEMOLECULES",false,"if this is a vector of positions do you want to make the positions into a whole before");
122 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");
123 936 : keys.addOutputComponent("x","default","the x-component of the atom position");
124 936 : keys.addOutputComponent("y","default","the y-component of the atom position");
125 936 : keys.addOutputComponent("z","default","the z-component of the atom position");
126 936 : keys.addOutputComponent("a","SCALED_COMPONENTS","the normalized projection on the first lattice vector of the atom position");
127 936 : keys.addOutputComponent("b","SCALED_COMPONENTS","the normalized projection on the second lattice vector of the atom position");
128 936 : keys.addOutputComponent("c","SCALED_COMPONENTS","the normalized projection on the third lattice vector of the atom position");
129 936 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
130 468 : }
131 :
132 94 : Position::Position(const ActionOptions&ao):
133 : PLUMED_COLVAR_INIT(ao),
134 94 : scaled_components(false),
135 94 : pbc(true),
136 94 : value(3),
137 95 : derivs(3),
138 188 : virial(3)
139 : {
140 376 : for(unsigned i=0; i<3; ++i) derivs[i].resize(1);
141 94 : std::vector<AtomNumber> atoms; parseAtomList(-1,atoms,this);
142 93 : unsigned mode=getModeAndSetupValues(this);
143 93 : if( mode==1 ) scaled_components=true;
144 :
145 93 : bool nopbc=!pbc;
146 94 : parseFlag("NOPBC",nopbc);
147 93 : pbc=!nopbc;
148 93 : checkRead();
149 :
150 93 : if(pbc) log.printf(" using periodic boundary conditions\n");
151 5 : else log.printf(" without periodic boundary conditions\n");
152 :
153 93 : requestAtoms(atoms);
154 96 : }
155 :
156 102 : void Position::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
157 204 : aa->parseAtomList("ATOM",num,t);
158 102 : if( t.size()==1 ) aa->log.printf(" for atom %d\n",t[0].serial());
159 3 : else if( num<0 || t.size()!=0 ) aa->error("Number of specified atoms should be 1");
160 101 : }
161 :
162 134 : unsigned Position::getModeAndSetupValues( ActionWithValue* av ) {
163 134 : bool sc; av->parseFlag("SCALED_COMPONENTS",sc);
164 134 : if(sc) {
165 45 : av->addComponentWithDerivatives("a"); av->componentIsPeriodic("a","-0.5","+0.5");
166 45 : av->addComponentWithDerivatives("b"); av->componentIsPeriodic("b","-0.5","+0.5");
167 45 : av->addComponentWithDerivatives("c"); av->componentIsPeriodic("c","-0.5","+0.5");
168 15 : return 1;
169 : }
170 238 : av->addComponentWithDerivatives("x"); av->componentIsNotPeriodic("x");
171 238 : av->addComponentWithDerivatives("y"); av->componentIsNotPeriodic("y");
172 238 : av->addComponentWithDerivatives("z"); av->componentIsNotPeriodic("z");
173 119 : av->log<<" WARNING: components will not have the proper periodicity - see manual\n";
174 : return 0;
175 : }
176 :
177 : // calculator
178 8078 : void Position::calculate() {
179 :
180 8078 : std::vector<Vector> distance(1);
181 8078 : if(pbc) {
182 16044 : distance[0]=pbcDistance(Vector(0.0,0.0,0.0),getPosition(0));
183 : } else {
184 56 : distance[0]=delta(Vector(0.0,0.0,0.0),getPosition(0));
185 : }
186 :
187 8078 : if(scaled_components) {
188 56 : calculateCV( 1, masses, charges, distance, value, derivs, virial, this );
189 56 : Value* valuea=getPntrToComponent("a");
190 56 : Value* valueb=getPntrToComponent("b");
191 56 : Value* valuec=getPntrToComponent("c");
192 56 : setAtomsDerivatives (valuea,0,derivs[0][0]);
193 56 : valuea->set(value[0]);
194 56 : setAtomsDerivatives (valueb,0,derivs[1][0]);
195 56 : valueb->set(value[1]);
196 56 : setAtomsDerivatives (valuec,0,derivs[2][0]);
197 56 : valuec->set(value[2]);
198 : } else {
199 8022 : calculateCV( 0, masses, charges, distance, value, derivs, virial, this );
200 8022 : Value* valuex=getPntrToComponent("x");
201 8022 : Value* valuey=getPntrToComponent("y");
202 8022 : Value* valuez=getPntrToComponent("z");
203 :
204 8022 : setAtomsDerivatives (valuex,0,derivs[0][0]);
205 8022 : setBoxDerivatives (valuex,virial[0]);
206 8022 : valuex->set(value[0]);
207 :
208 8022 : setAtomsDerivatives (valuey,0,derivs[1][0]);
209 8022 : setBoxDerivatives (valuey,virial[1]);
210 8022 : valuey->set(value[1]);
211 :
212 8022 : setAtomsDerivatives (valuez,0,derivs[2][0]);
213 8022 : setBoxDerivatives (valuez,virial[2]);
214 8022 : valuez->set(value[2]);
215 : }
216 8078 : }
217 :
218 151422 : void Position::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
219 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
220 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
221 151422 : if( mode==1 ) {
222 10841 : Vector d=aa->getPbc().realToScaled(pos[0]);
223 10841 : vals[0]=Tools::pbc(d[0]); vals[1]=Tools::pbc(d[1]); vals[2]=Tools::pbc(d[2]);
224 10841 : derivs[0][0]=matmul(aa->getPbc().getInvBox(),Vector(+1,0,0));
225 10841 : derivs[1][0]=matmul(aa->getPbc().getInvBox(),Vector(0,+1,0));
226 10841 : derivs[2][0]=matmul(aa->getPbc().getInvBox(),Vector(0,0,+1));
227 : } else {
228 562324 : for(unsigned i=0; i<3; ++i) vals[i]=pos[0][i];
229 140581 : derivs[0][0]=Vector(+1,0,0); derivs[1][0]=Vector(0,+1,0); derivs[2][0]=Vector(0,0,+1);
230 140581 : virial[0]=Tensor(pos[0],Vector(-1,0,0)); virial[1]=Tensor(pos[0],Vector(0,-1,0)); virial[2]=Tensor(pos[0],Vector(0,0,-1));
231 : }
232 151422 : }
233 :
234 : }
235 : }
236 :
237 :
238 :
|
