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