Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2020 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 138 : 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 : virtual void calculate();
87 : };
88 :
89 7495 : PLUMED_REGISTER_ACTION(Position,"POSITION")
90 :
91 71 : void Position::registerKeywords( Keywords& keys ) {
92 71 : Colvar::registerKeywords( keys );
93 71 : componentsAreNotOptional(keys);
94 284 : keys.add("atoms","ATOM","the atom number");
95 213 : 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 284 : keys.addOutputComponent("x","default","the x-component of the atom position");
97 284 : keys.addOutputComponent("y","default","the y-component of the atom position");
98 284 : keys.addOutputComponent("z","default","the z-component of the atom position");
99 284 : keys.addOutputComponent("a","SCALED_COMPONENTS","the normalized projection on the first lattice vector of the atom position");
100 284 : keys.addOutputComponent("b","SCALED_COMPONENTS","the normalized projection on the second lattice vector of the atom position");
101 284 : keys.addOutputComponent("c","SCALED_COMPONENTS","the normalized projection on the third lattice vector of the atom position");
102 71 : }
103 :
104 70 : Position::Position(const ActionOptions&ao):
105 : PLUMED_COLVAR_INIT(ao),
106 : scaled_components(false),
107 71 : pbc(true)
108 : {
109 : vector<AtomNumber> atoms;
110 140 : parseAtomList("ATOM",atoms);
111 70 : if(atoms.size()!=1)
112 2 : error("Number of specified atoms should be 1");
113 138 : parseFlag("SCALED_COMPONENTS",scaled_components);
114 69 : bool nopbc=!pbc;
115 138 : parseFlag("NOPBC",nopbc);
116 69 : pbc=!nopbc;
117 69 : checkRead();
118 :
119 138 : log.printf(" for atom %d\n",atoms[0].serial());
120 69 : if(pbc) log.printf(" using periodic boundary conditions\n");
121 3 : else log.printf(" without periodic boundary conditions\n");
122 :
123 69 : if(scaled_components) {
124 20 : addComponentWithDerivatives("a"); componentIsPeriodic("a","-0.5","+0.5");
125 20 : addComponentWithDerivatives("b"); componentIsPeriodic("b","-0.5","+0.5");
126 20 : addComponentWithDerivatives("c"); componentIsPeriodic("c","-0.5","+0.5");
127 : } else {
128 195 : addComponentWithDerivatives("x"); componentIsNotPeriodic("x");
129 195 : addComponentWithDerivatives("y"); componentIsNotPeriodic("y");
130 195 : addComponentWithDerivatives("z"); componentIsNotPeriodic("z");
131 65 : log<<" WARNING: components will not have the proper periodicity - see manual\n";
132 : }
133 :
134 69 : requestAtoms(atoms);
135 69 : }
136 :
137 :
138 : // calculator
139 4497 : void Position::calculate() {
140 :
141 4497 : Vector distance;
142 4497 : if(pbc) {
143 8964 : distance=pbcDistance(Vector(0.0,0.0,0.0),getPosition(0));
144 : } else {
145 15 : distance=delta(Vector(0.0,0.0,0.0),getPosition(0));
146 : }
147 :
148 4497 : if(scaled_components) {
149 82 : Value* valuea=getPntrToComponent("a");
150 82 : Value* valueb=getPntrToComponent("b");
151 82 : Value* valuec=getPntrToComponent("c");
152 41 : Vector d=getPbc().realToScaled(distance);
153 82 : setAtomsDerivatives (valuea,0,matmul(getPbc().getInvBox(),Vector(+1,0,0)));
154 41 : valuea->set(Tools::pbc(d[0]));
155 82 : setAtomsDerivatives (valueb,0,matmul(getPbc().getInvBox(),Vector(0,+1,0)));
156 41 : valueb->set(Tools::pbc(d[1]));
157 82 : setAtomsDerivatives (valuec,0,matmul(getPbc().getInvBox(),Vector(0,0,+1)));
158 41 : valuec->set(Tools::pbc(d[2]));
159 : } else {
160 8912 : Value* valuex=getPntrToComponent("x");
161 8912 : Value* valuey=getPntrToComponent("y");
162 8912 : Value* valuez=getPntrToComponent("z");
163 :
164 4456 : setAtomsDerivatives (valuex,0,Vector(+1,0,0));
165 4456 : setBoxDerivatives (valuex,Tensor(distance,Vector(-1,0,0)));
166 4456 : valuex->set(distance[0]);
167 :
168 4456 : setAtomsDerivatives (valuey,0,Vector(0,+1,0));
169 4456 : setBoxDerivatives (valuey,Tensor(distance,Vector(0,-1,0)));
170 4456 : valuey->set(distance[1]);
171 :
172 4456 : setAtomsDerivatives (valuez,0,Vector(0,0,+1));
173 4456 : setBoxDerivatives (valuez,Tensor(distance,Vector(0,0,-1)));
174 4456 : valuez->set(distance[2]);
175 : }
176 4497 : }
177 :
178 : }
179 5517 : }
180 :
181 :
182 :
|
