1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2014-2021 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 "ActionRegister.h"
24 #include "tools/Pbc.h"
25
26 #include <string>
27 #include <cmath>
28
29 using namespace std;
30
31 namespace PLMD {
32 namespace colvar {
33
34 //+PLUMEDOC COLVAR POSITION
35 /*
36 Calculate the components of the position of an atom.
37
38 Notice that single components will not have the proper periodicity!
39 If you need the values to be consistent through PBC you should use SCALED_COMPONENTS,
40 which defines values that by construction are in the -0.5,0.5 domain. This is
41 similar to the equivalent flag for \ref DISTANCE.
42 Also notice that by default the minimal image distance from the
43 origin is considered (can be changed with NOPBC).
44
45 \attention
46 This variable should be used with extreme care since it allows to easily go into troubles. See comments below.
47
48 This variable can be safely used only if
49 Hamiltonian is not invariant for translation (i.e. there are other absolute positions which are biased, e.g. by position restraints)
50 and cell size and shapes are fixed through the simulation.
51
52 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.
53 This can be done e.g. using \ref FIT_TO_TEMPLATE.
54
55 \par Examples
56
57 \plumedfile
58 # align to a template
59 FIT_TO_TEMPLATE REFERENCE=ref.pdb
60 p: POSITION ATOM=3
61 PRINT ARG=p.x,p.y,p.z
62 \endplumedfile
63
64 The reference position is specified in a pdb file like the one shown below
65
66 \auxfile{ref.pdb}
67 ATOM 3 HT3 ALA 2 -1.480 -1.560 1.212 1.00 1.00 DIA H
68 ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 1.00 DIA C
69 ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 1.00 DIA H
70 ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 1.00 DIA H
71 ATOM 14 OY ALA 2 -1.139 0.931 -0.973 1.00 1.00 DIA O
72 END
73 \endauxfile
74
75 */
76 //+ENDPLUMEDOC
77
78 class Position : public Colvar {
79 bool scaled_components;
80 bool pbc;
81
82 public:
83 static void registerKeywords( Keywords& keys );
84 explicit Position(const ActionOptions&);
85 // active methods:
86 void calculate() override;
87 };
88
89 PLUMED_REGISTER_ACTION(Position,"POSITION")
90
registerKeywords(Keywords & keys)91 void Position::registerKeywords( Keywords& keys ) {
92 Colvar::registerKeywords( keys );
93 componentsAreNotOptional(keys);
94 keys.add("atoms","ATOM","the atom number");
95 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");
96 keys.addOutputComponent("x","default","the x-component of the atom position");
97 keys.addOutputComponent("y","default","the y-component of the atom position");
98 keys.addOutputComponent("z","default","the z-component of the atom position");
99 keys.addOutputComponent("a","SCALED_COMPONENTS","the normalized projection on the first lattice vector of the atom position");
100 keys.addOutputComponent("b","SCALED_COMPONENTS","the normalized projection on the second lattice vector of the atom position");
101 keys.addOutputComponent("c","SCALED_COMPONENTS","the normalized projection on the third lattice vector of the atom position");
102 }
103
Position(const ActionOptions & ao)104 Position::Position(const ActionOptions&ao):
105 PLUMED_COLVAR_INIT(ao),
106 scaled_components(false),
107 pbc(true)
108 {
109 vector<AtomNumber> atoms;
110 parseAtomList("ATOM",atoms);
111 if(atoms.size()!=1)
112 error("Number of specified atoms should be 1");
113 parseFlag("SCALED_COMPONENTS",scaled_components);
114 bool nopbc=!pbc;
115 parseFlag("NOPBC",nopbc);
116 pbc=!nopbc;
117 checkRead();
118
119 log.printf(" for atom %d\n",atoms[0].serial());
120 if(pbc) log.printf(" using periodic boundary conditions\n");
121 else log.printf(" without periodic boundary conditions\n");
122
123 if(scaled_components) {
124 addComponentWithDerivatives("a"); componentIsPeriodic("a","-0.5","+0.5");
125 addComponentWithDerivatives("b"); componentIsPeriodic("b","-0.5","+0.5");
126 addComponentWithDerivatives("c"); componentIsPeriodic("c","-0.5","+0.5");
127 } else {
128 addComponentWithDerivatives("x"); componentIsNotPeriodic("x");
129 addComponentWithDerivatives("y"); componentIsNotPeriodic("y");
130 addComponentWithDerivatives("z"); componentIsNotPeriodic("z");
131 log<<" WARNING: components will not have the proper periodicity - see manual\n";
132 }
133
134 requestAtoms(atoms);
135 }
136
137
138 // calculator
calculate()139 void Position::calculate() {
140
141 Vector distance;
142 if(pbc) {
143 distance=pbcDistance(Vector(0.0,0.0,0.0),getPosition(0));
144 } else {
145 distance=delta(Vector(0.0,0.0,0.0),getPosition(0));
146 }
147
148 if(scaled_components) {
149 Value* valuea=getPntrToComponent("a");
150 Value* valueb=getPntrToComponent("b");
151 Value* valuec=getPntrToComponent("c");
152 Vector d=getPbc().realToScaled(distance);
153 setAtomsDerivatives (valuea,0,matmul(getPbc().getInvBox(),Vector(+1,0,0)));
154 valuea->set(Tools::pbc(d[0]));
155 setAtomsDerivatives (valueb,0,matmul(getPbc().getInvBox(),Vector(0,+1,0)));
156 valueb->set(Tools::pbc(d[1]));
157 setAtomsDerivatives (valuec,0,matmul(getPbc().getInvBox(),Vector(0,0,+1)));
158 valuec->set(Tools::pbc(d[2]));
159 } else {
160 Value* valuex=getPntrToComponent("x");
161 Value* valuey=getPntrToComponent("y");
162 Value* valuez=getPntrToComponent("z");
163
164 setAtomsDerivatives (valuex,0,Vector(+1,0,0));
165 setBoxDerivatives (valuex,Tensor(distance,Vector(-1,0,0)));
166 valuex->set(distance[0]);
167
168 setAtomsDerivatives (valuey,0,Vector(0,+1,0));
169 setBoxDerivatives (valuey,Tensor(distance,Vector(0,-1,0)));
170 valuey->set(distance[1]);
171
172 setAtomsDerivatives (valuez,0,Vector(0,0,+1));
173 setBoxDerivatives (valuez,Tensor(distance,Vector(0,0,-1)));
174 valuez->set(distance[2]);
175 }
176 }
177
178 }
179 }
180
181
182
183