Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2018 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 "ActionRegister.h"
24 : #include "tools/Vector.h"
25 : #include "tools/Exception.h"
26 :
27 : using namespace std;
28 :
29 : namespace PLMD {
30 : namespace vatom {
31 :
32 : //+PLUMEDOC VATOM GHOST
33 : /*
34 : Calculate the absolute position of a ghost atom with fixed coordinates
35 : in the local reference frame formed by three atoms.
36 : The computed ghost atom is stored as a virtual atom that can be accessed in
37 : an atom list through the the label for the GHOST action that creates it.
38 :
39 : \par Examples
40 :
41 : The following input instructs plumed to print the distance between the
42 : ghost atom and the center of mass for atoms 15,20:
43 : \verbatim
44 : GHOST ATOMS=1,5,10 COORDINATES=10.0,10.0,10.0 LABEL=c1
45 : COM ATOMS=15,20 LABEL=c2
46 : DISTANCE ATOMS=c1,c2 LABEL=d1
47 : PRINT ARG=d1
48 : \endverbatim
49 : (See also \ref DISTANCE and \ref PRINT).
50 :
51 : */
52 : //+ENDPLUMEDOC
53 :
54 :
55 6 : class Ghost:
56 : public ActionWithVirtualAtom
57 : {
58 : vector<double> coord;
59 : public:
60 : explicit Ghost(const ActionOptions&ao);
61 : void calculate();
62 : static void registerKeywords( Keywords& keys );
63 : };
64 :
65 2526 : PLUMED_REGISTER_ACTION(Ghost,"GHOST")
66 :
67 4 : void Ghost::registerKeywords(Keywords& keys) {
68 4 : ActionWithVirtualAtom::registerKeywords(keys);
69 4 : keys.add("atoms","COORDINATES","coordinates of the ghost atom in the local reference frame");
70 4 : }
71 :
72 3 : Ghost::Ghost(const ActionOptions&ao):
73 : Action(ao),
74 3 : ActionWithVirtualAtom(ao)
75 : {
76 3 : vector<AtomNumber> atoms;
77 3 : parseAtomList("ATOMS",atoms);
78 3 : if(atoms.size()!=3) error("ATOMS should contain a list of three atoms");
79 :
80 3 : parseVector("COORDINATES",coord);
81 3 : if(coord.size()!=3) error("COORDINATES should be a list of three real numbers");
82 :
83 3 : checkRead();
84 3 : log.printf(" of atoms");
85 3 : for(unsigned i=0; i<atoms.size(); ++i) log.printf(" %d",atoms[i].serial());
86 3 : log.printf("\n");
87 3 : requestAtoms(atoms);
88 3 : }
89 :
90 7 : void Ghost::calculate() {
91 7 : Vector pos;
92 7 : vector<Tensor> deriv(getNumberOfAtoms());
93 14 : vector<Vector> n;
94 :
95 : // first versor
96 7 : Vector n01 = delta(getPosition(0), getPosition(1));
97 7 : n.push_back(n01/n01.modulo());
98 :
99 : // auxiliary vector
100 7 : Vector n02 = delta(getPosition(0), getPosition(2));
101 :
102 : // second versor
103 7 : Vector n03 = crossProduct(n[0],n02);
104 7 : double n03_norm = n03.modulo();
105 7 : n.push_back(n03/n03_norm);
106 :
107 : // third versor
108 7 : n.push_back(crossProduct(n[0],n[1]));
109 :
110 : // origin of the reference system
111 7 : pos = getPosition(0);
112 :
113 28 : for(unsigned i=0; i<3; ++i) {
114 21 : pos += coord[i] * n[i];
115 : }
116 :
117 7 : setPosition(pos);
118 7 : setMass(1.0);
119 7 : setCharge(0.0);
120 :
121 : // some useful tensors for derivatives
122 7 : Tensor dn0d0 = (-Tensor::identity()+Tensor(n[0],n[0]))/n01.modulo();
123 7 : Tensor dn0d1 = (+Tensor::identity()-Tensor(n[0],n[0]))/n01.modulo();
124 7 : Tensor dn02d0 = -Tensor::identity();
125 7 : Tensor dn02d2 = Tensor::identity();
126 :
127 : // derivative of n1 = n0 x n02
128 7 : Tensor dn1d0, dn1d1, dn1d2;
129 7 : Vector aux0, aux1, aux2;
130 :
131 28 : for(unsigned j=0; j<3; ++j) {
132 : // derivative of n0 x n02 with respect to point 0, coordinate j
133 21 : Vector tmp00 = Vector( dn0d0(j,0), dn0d0(j,1), dn0d0(j,2));
134 21 : Vector tmp020 = Vector(dn02d0(j,0), dn02d0(j,1), dn02d0(j,2));
135 21 : Vector tmp0 = crossProduct(tmp00,n02) + crossProduct(n[0],tmp020);
136 21 : aux0[j] = dotProduct(tmp0,n[1]);
137 : // derivative of n0 x n02 with respect to point 1, coordinate j
138 21 : Vector tmp01 = Vector( dn0d1(j,0), dn0d1(j,1), dn0d1(j,2));
139 21 : Vector tmp1 = crossProduct(tmp01,n02);
140 21 : aux1[j] = dotProduct(tmp1,n[1]);
141 : // derivative of n0 x n02 with respect to point 2, coordinate j
142 21 : Vector tmp022 = Vector(dn02d2(j,0), dn02d2(j,1), dn02d2(j,2));
143 21 : Vector tmp2 = crossProduct(n[0],tmp022);
144 21 : aux2[j] = dotProduct(tmp2,n[1]);
145 : // derivative of n1 = (n0 x n02) / || (n0 x n02) ||
146 84 : for(unsigned i=0; i<3; ++i) {
147 63 : dn1d0(j,i) = ( tmp0[i] - aux0[j] * n[1][i] ) / n03_norm;
148 63 : dn1d1(j,i) = ( tmp1[i] - aux1[j] * n[1][i] ) / n03_norm;
149 63 : dn1d2(j,i) = ( tmp2[i] - aux2[j] * n[1][i] ) / n03_norm;
150 : }
151 : }
152 :
153 : // Derivative of the last versor n2 = n0 x n1 = ( n0( n0 n02 ) - n02 ) / || n0 x n02 ||
154 : // Scalar product and derivatives
155 7 : double n0_n02 = dotProduct(n[0],n02);
156 7 : Vector dn0_n02d0, dn0_n02d1, dn0_n02d2;
157 :
158 28 : for(unsigned j=0; j<3; ++j) {
159 84 : for(unsigned i=0; i<3; ++i) {
160 63 : dn0_n02d0[j] += dn0d0(j,i)*n02[i] + n[0][i]*dn02d0(j,i);
161 63 : dn0_n02d1[j] += dn0d1(j,i)*n02[i];
162 63 : dn0_n02d2[j] += n[0][i]*dn02d2(j,i);
163 : }
164 : }
165 :
166 7 : Tensor dn2d0, dn2d1, dn2d2;
167 28 : for(unsigned j=0; j<3; ++j) {
168 84 : for(unsigned i=0; i<3; ++i) {
169 63 : 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;
170 63 : 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;
171 63 : 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;
172 : }
173 : }
174 :
175 : // Finally, the derivative tensor
176 7 : deriv[0] = Tensor::identity() + coord[0]*dn0d0 + coord[1]*dn1d0 + coord[2]*dn2d0;
177 7 : deriv[1] = coord[0]*dn0d1 + coord[1]*dn1d1 + coord[2]*dn2d1;
178 7 : deriv[2] = coord[1]*dn1d2 + coord[2]*dn2d2;
179 :
180 7 : setAtomsDerivatives(deriv);
181 :
182 : // Virial contribution
183 14 : setBoxDerivativesNoPbc();
184 7 : }
185 :
186 : }
187 2523 : }
|
