Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "ActionWithVirtualAtom.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/Vector.h"
25 : #include "tools/Exception.h"
26 : #include <array>
27 :
28 : namespace PLMD {
29 : namespace vatom {
30 :
31 : //+PLUMEDOC VATOM GHOST
32 : /*
33 : Calculate the absolute position of a ghost atom with fixed coordinates in the local reference frame formed by three atoms.
34 :
35 : This action allows you to create a virtual atom that has a fixed set of coordinates in a local reference
36 : frame that is formed by three atoms. The way that this action is used is illustrated below:
37 :
38 : ```plumed
39 : c1: GHOST ATOMS=1,5,10 COORDINATES=10.0,10.0,10.0
40 : c2: COM ATOMS=15,20
41 : d1: DISTANCE ATOMS=c1,c2
42 : PRINT ARG=d1
43 : ```
44 :
45 : Notice that ghost atom's position is stored as [a virtual atom](specifying_atoms.md). The position of this atom can thus be
46 : used in the DISTANCE command by using the label for the GHOST action.
47 :
48 : The position of the ghost atom `c1` for the input above is:
49 :
50 : $$
51 : r_{c1} = r_1 + 10\hat{a} + 10c + 10 \hat{b} 10\hat{c}
52 : $$
53 :
54 : where unit vectors, $\hat{a}$, $\hat{b}$ and $\hat{c}$ in the expression above are obtained by dividing each
55 : of the three (orthogonal) vectors below by their magnitudes:
56 :
57 : $$
58 : a = (r_5-r_1) \quad b = (r_5-r_1) \time (r_{10}-r_1) \quad (r_5-r_1)\times b
59 : $$
60 :
61 : In all these expressions $r_i$ is used to indicate the position of atom $i$. If you run with periodic boundary conditions
62 : these are taken into account automatically when computing the differences between position vectors above. The way this is
63 : handled is akin to the way molecules are rebuilt in the [WHOLEMOLECULES](WHOLEMOLECULES.md) command. For the example above
64 : this would ensure that atom 5 is shifted to the periodic image where it is closest to atom 1 and atom 10 is shifted to the
65 : periodic image where it is closest to atom 10. If you wish to
66 : turn off this behaviour and you wish to disregard the periodic boundaries when computing these differences you should use
67 : the NOPBC flag.
68 :
69 : */
70 : //+ENDPLUMEDOC
71 :
72 :
73 : class Ghost:
74 : public ActionWithVirtualAtom {
75 : std::vector<double> coord;
76 : std::vector<Tensor> deriv;
77 : bool nopbc=false;
78 : public:
79 : explicit Ghost(const ActionOptions&ao);
80 : void calculate() override;
81 : static void registerKeywords( Keywords& keys );
82 : };
83 :
84 : PLUMED_REGISTER_ACTION(Ghost,"GHOST")
85 :
86 711 : void Ghost::registerKeywords(Keywords& keys) {
87 711 : ActionWithVirtualAtom::registerKeywords(keys);
88 711 : keys.add("atoms","COORDINATES","coordinates of the ghost atom in the local reference frame");
89 711 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
90 711 : }
91 :
92 709 : Ghost::Ghost(const ActionOptions&ao):
93 : Action(ao),
94 709 : ActionWithVirtualAtom(ao) {
95 : std::vector<AtomNumber> atoms;
96 1418 : parseAtomList("ATOMS",atoms);
97 709 : if(atoms.size()!=3) {
98 0 : error("ATOMS should contain a list of three atoms");
99 : }
100 :
101 1418 : parseVector("COORDINATES",coord);
102 709 : if(coord.size()!=3) {
103 0 : error("COORDINATES should be a list of three real numbers");
104 : }
105 :
106 709 : parseFlag("NOPBC",nopbc);
107 :
108 709 : checkRead();
109 709 : log.printf(" of atoms");
110 2836 : for(unsigned i=0; i<atoms.size(); ++i) {
111 2127 : log.printf(" %d",atoms[i].serial());
112 : }
113 709 : log.printf("\n");
114 :
115 709 : if(nopbc) {
116 4 : log<<" PBC will be ignored\n";
117 : } else {
118 705 : log<<" broken molecules will be rebuilt assuming atoms are in the proper order\n";
119 : }
120 709 : requestAtoms(atoms);
121 709 : }
122 :
123 1413 : void Ghost::calculate() {
124 :
125 1413 : if(!nopbc) {
126 1409 : makeWhole();
127 : }
128 :
129 1413 : Vector pos;
130 1413 : deriv.resize(getNumberOfAtoms());
131 1413 : std::array<Vector,3> n;
132 :
133 : // first versor
134 1413 : Vector n01 = delta(getPosition(0), getPosition(1));
135 1413 : n[0]=n01/n01.modulo();
136 :
137 : // auxiliary vector
138 1413 : Vector n02 = delta(getPosition(0), getPosition(2));
139 :
140 : // second versor
141 1413 : Vector n03 = crossProduct(n[0],n02);
142 1413 : double n03_norm = n03.modulo();
143 1413 : n[1]=n03/n03_norm;
144 :
145 : // third versor
146 1413 : n[2]=crossProduct(n[0],n[1]);
147 :
148 : // origin of the reference system
149 1413 : pos = getPosition(0);
150 :
151 5652 : for(unsigned i=0; i<3; ++i) {
152 4239 : pos += coord[i] * n[i];
153 : }
154 :
155 1413 : setPosition(pos);
156 1413 : setMass(1.0);
157 1413 : setCharge(0.0);
158 :
159 : // some useful tensors for derivatives
160 1413 : Tensor dn0d0 = (-Tensor::identity()+Tensor(n[0],n[0]))/n01.modulo();
161 1413 : Tensor dn0d1 = (+Tensor::identity()-Tensor(n[0],n[0]))/n01.modulo();
162 1413 : Tensor dn02d0 = -Tensor::identity();
163 1413 : Tensor dn02d2 = Tensor::identity();
164 :
165 : // derivative of n1 = n0 x n02
166 1413 : Tensor dn1d0, dn1d1, dn1d2;
167 1413 : Vector aux0, aux1, aux2;
168 :
169 5652 : for(unsigned j=0; j<3; ++j) {
170 : // derivative of n0 x n02 with respect to point 0, coordinate j
171 4239 : Vector tmp00 = Vector( dn0d0(j,0), dn0d0(j,1), dn0d0(j,2));
172 4239 : Vector tmp020 = Vector(dn02d0(j,0), dn02d0(j,1), dn02d0(j,2));
173 4239 : Vector tmp0 = crossProduct(tmp00,n02) + crossProduct(n[0],tmp020);
174 4239 : aux0[j] = dotProduct(tmp0,n[1]);
175 : // derivative of n0 x n02 with respect to point 1, coordinate j
176 4239 : Vector tmp01 = Vector( dn0d1(j,0), dn0d1(j,1), dn0d1(j,2));
177 4239 : Vector tmp1 = crossProduct(tmp01,n02);
178 4239 : aux1[j] = dotProduct(tmp1,n[1]);
179 : // derivative of n0 x n02 with respect to point 2, coordinate j
180 4239 : Vector tmp022 = Vector(dn02d2(j,0), dn02d2(j,1), dn02d2(j,2));
181 4239 : Vector tmp2 = crossProduct(n[0],tmp022);
182 4239 : aux2[j] = dotProduct(tmp2,n[1]);
183 : // derivative of n1 = (n0 x n02) / || (n0 x n02) ||
184 16956 : for(unsigned i=0; i<3; ++i) {
185 12717 : dn1d0(j,i) = ( tmp0[i] - aux0[j] * n[1][i] ) / n03_norm;
186 12717 : dn1d1(j,i) = ( tmp1[i] - aux1[j] * n[1][i] ) / n03_norm;
187 12717 : dn1d2(j,i) = ( tmp2[i] - aux2[j] * n[1][i] ) / n03_norm;
188 : }
189 : }
190 :
191 : // Derivative of the last versor n2 = n0 x n1 = ( n0( n0 n02 ) - n02 ) / || n0 x n02 ||
192 : // Scalar product and derivatives
193 1413 : double n0_n02 = dotProduct(n[0],n02);
194 1413 : Vector dn0_n02d0, dn0_n02d1, dn0_n02d2;
195 :
196 5652 : for(unsigned j=0; j<3; ++j) {
197 16956 : for(unsigned i=0; i<3; ++i) {
198 12717 : dn0_n02d0[j] += dn0d0(j,i)*n02[i] + n[0][i]*dn02d0(j,i);
199 12717 : dn0_n02d1[j] += dn0d1(j,i)*n02[i];
200 12717 : dn0_n02d2[j] += n[0][i]*dn02d2(j,i);
201 : }
202 : }
203 :
204 1413 : Tensor dn2d0, dn2d1, dn2d2;
205 5652 : for(unsigned j=0; j<3; ++j) {
206 16956 : for(unsigned i=0; i<3; ++i) {
207 12717 : dn2d0(j,i) = ( dn0d0(j,i) * n0_n02 + n[0][i] * dn0_n02d0[j] - dn02d0(j,i) - ( n[0][i] * n0_n02 - n02[i] ) * aux0[j] / n03_norm ) / n03_norm;
208 12717 : dn2d1(j,i) = ( dn0d1(j,i) * n0_n02 + n[0][i] * dn0_n02d1[j] - ( n[0][i] * n0_n02 - n02[i] ) * aux1[j] / n03_norm ) / n03_norm;
209 12717 : dn2d2(j,i) = ( n[0][i] * dn0_n02d2[j] - dn02d2(j,i) - ( n[0][i] * n0_n02 - n02[i] ) * aux2[j] / n03_norm ) / n03_norm;
210 : }
211 : }
212 :
213 : // Finally, the derivative tensor
214 1413 : deriv[0] = Tensor::identity() + coord[0]*dn0d0 + coord[1]*dn1d0 + coord[2]*dn2d0;
215 1413 : deriv[1] = coord[0]*dn0d1 + coord[1]*dn1d1 + coord[2]*dn2d1;
216 1413 : deriv[2] = coord[1]*dn1d2 + coord[2]*dn2d2;
217 :
218 1413 : setAtomsDerivatives(deriv);
219 :
220 : // Virial contribution
221 1413 : setBoxDerivativesNoPbc();
222 1413 : }
223 :
224 : }
225 : }
|
