Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "core/ActionRegister.h"
23 : #include "core/PlumedMain.h"
24 : #include "tools/Units.h"
25 : #include "tools/Pbc.h"
26 : #include "ActionVolume.h"
27 : #include "VolumeShortcut.h"
28 :
29 : //+PLUMEDOC VOLUMES TETRAHEDRALPORE
30 : /*
31 : This quantity can be used to calculate functions of the distribution of collective variables for the atoms lie that lie in a box defined by the positions of four atoms at the corners of a tetrahedron.
32 :
33 : Each of the base quantities calculated by a multicolvar can can be assigned to a particular point in three
34 : dimensional space. For example, if we have the coordination numbers for all the atoms in the
35 : system each coordination number can be assumed to lie on the position of the central atom.
36 : Because each base quantity can be assigned to a particular point in space we can calculate functions of the
37 : distribution of base quantities in a particular part of the box by using:
38 :
39 : \f[
40 : \overline{s}_{\tau} = \frac{ \sum_i f(s_i) w(u_i,v_i,w_i) }{ \sum_i w(u_i,v_i,w_i) }
41 : \f]
42 :
43 : where the sum is over the collective variables, \f$s_i\f$, each of which can be thought to be at \f$ (u_i,v_i,z_i)\f$.
44 : The function \f$(s_i)\f$ can be any of the usual LESS_THAN, MORE_THAN, WITHIN etc that are used in all other multicolvars.
45 : Notice that here (at variance with what is done in \ref AROUND) we have transformed from the usual \f$(x_i,y_i,z_i)\f$
46 : position to a position in \f$ (u_i,v_i,z_i)\f$. This is done using a rotation matrix as follows:
47 :
48 : \f[
49 : \left(
50 : \begin{matrix}
51 : u_i \\
52 : v_i \\
53 : w_i
54 : \end{matrix}
55 : \right) = \mathbf{R}
56 : \left(
57 : \begin{matrix}
58 : x_i - x_o \\
59 : y_i - y_o \\
60 : z_i - z_o
61 : \end{matrix}
62 : \right)
63 : \f]
64 :
65 : where \f$\mathbf{R}\f$ is a rotation matrix that is calculated by constructing a set of three orthonormal vectors from the
66 : reference positions specified by the user. Initially unit vectors are found by calculating the bisector, \f$\mathbf{b}\f$, and
67 : cross product, \f$\mathbf{c}\f$, of the vectors connecting atoms 1 and 2. A third unit vector, \f$\mathbf{p}\f$ is then found by taking the cross
68 : product between the cross product calculated during the first step, \f$\mathbf{c}\f$ and the bisector, \f$\mathbf{b}\f$. From this
69 : second cross product \f$\mathbf{p}\f$ and the bisector \f$\mathbf{b}\f$ two new vectors are calculated using:
70 :
71 : \f[
72 : v_1 = \cos\left(\frac{\pi}{4}\right)\mathbf{b} + \sin\left(\frac{\pi}{4}\right)\mathbf{p} \qquad \textrm{and} \qquad
73 : v_2 = \cos\left(\frac{\pi}{4}\right)\mathbf{b} - \sin\left(\frac{\pi}{4}\right)\mathbf{p}
74 : \f]
75 :
76 : In the previous function \f$ w(u_i,v_i,w_i) \f$ measures whether or not the system is in the subregion of interest. It
77 : is equal to:
78 :
79 : \f[
80 : w(u_i,v_i,w_i) = \int_{0}^{u'} \int_{0}^{v'} \int_{0}^{w'} \textrm{d}u\textrm{d}v\textrm{d}w
81 : K\left( \frac{u - u_i}{\sigma} \right)K\left( \frac{v - v_i}{\sigma} \right)K\left( \frac{w - w_i}{\sigma} \right)
82 : \f]
83 :
84 : where \f$K\f$ is one of the kernel functions described on \ref histogrambead and \f$\sigma\f$ is a bandwidth parameter.
85 : The values of \f$u'\f$ and \f$v'\f$ are found by finding the projections of the vectors connecting atoms 1 and 2 and 1
86 : and 3 \f$v_1\f$ and \f$v_2\f$. This gives four projections: the largest two projections are used in the remainder of
87 : the calculations. \f$w'\f$ is calculated by taking the projection of the vector connecting atoms 1 and 4 on the vector
88 : \f$\mathbf{c}\f$. Notice that the manner by which this box is constructed differs from the way this is done in \ref CAVITY.
89 : This is in fact the only point of difference between these two actions.
90 :
91 : \par Examples
92 :
93 : The following commands tell plumed to calculate the number of atom inside a tetrahedral cavity. The extent of the tetrahedral
94 : cavity is calculated from the positions of atoms 1, 4, 5, and 11, The final value will be labeled cav.
95 :
96 : \plumedfile
97 : d1: DENSITY SPECIES=20-500
98 : TETRAHEDRALPORE DATA=d1 ATOMS=1,4,5,11 SIGMA=0.1 LABEL=cav
99 : \endplumedfile
100 :
101 : The following command tells plumed to calculate the coordination numbers (with other water molecules) for the water
102 : molecules in the tetrahedral cavity described above. The average coordination number and the number of coordination
103 : numbers more than 4 is then calculated. The values of these two quantities are given the labels cav.mean and cav.morethan
104 :
105 : \plumedfile
106 : d1: COORDINATIONNUMBER SPECIES=20-500 R_0=0.1
107 : CAVITY DATA=d1 ATOMS=1,4,5,11 SIGMA=0.1 MEAN MORE_THAN={RATIONAL R_0=4} LABEL=cav
108 : \endplumedfile
109 :
110 : */
111 : //+ENDPLUMEDOC
112 :
113 : //+PLUMEDOC MCOLVAR TETRAHEDRALPORE_CALC
114 : /*
115 : Calculate a vector from the input positions with elements equal to one when the positions are in a particular part of the cell and elements equal to zero otherwise
116 :
117 : \par Examples
118 :
119 : */
120 : //+ENDPLUMEDOC
121 :
122 : namespace PLMD {
123 : namespace volumes {
124 :
125 : class VolumeTetrapore : public ActionVolume {
126 : private:
127 : bool boxout;
128 : OFile boxfile;
129 : double lenunit;
130 : double jacob_det;
131 : double len_bi, len_cross, len_perp, sigma;
132 : Vector origin, bi, cross, perp;
133 : std::vector<Vector> dlbi, dlcross, dlperp;
134 : std::vector<Tensor> dbi, dcross, dperp;
135 : public:
136 : static void registerKeywords( Keywords& keys );
137 : explicit VolumeTetrapore(const ActionOptions& ao);
138 : ~VolumeTetrapore();
139 : void setupRegions() override;
140 : void update() override;
141 : double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const override;
142 : };
143 :
144 : PLUMED_REGISTER_ACTION(VolumeTetrapore,"TETRAHEDRALPORE_CALC")
145 : char glob_tetrapore[] = "TETRAHEDRALPORE";
146 : typedef VolumeShortcut<glob_tetrapore> VolumeTetraporeShortcut;
147 : PLUMED_REGISTER_ACTION(VolumeTetraporeShortcut,"TETRAHEDRALPORE")
148 :
149 4 : void VolumeTetrapore::registerKeywords( Keywords& keys ) {
150 4 : ActionVolume::registerKeywords( keys );
151 4 : keys.setDisplayName("TETRAHEDRALPORE");
152 8 : keys.add("atoms","BOX","the positions of four atoms that define spatial extent of the cavity");
153 8 : keys.addFlag("PRINT_BOX",false,"write out the positions of the corners of the box to an xyz file");
154 8 : keys.add("optional","FILE","the file on which to write out the box coordinates");
155 8 : keys.add("optional","UNITS","( default=nm ) the units in which to write out the corners of the box");
156 4 : }
157 :
158 0 : VolumeTetrapore::VolumeTetrapore(const ActionOptions& ao):
159 : Action(ao),
160 : ActionVolume(ao),
161 0 : boxout(false),
162 0 : lenunit(1.0),
163 0 : dlbi(4),
164 0 : dlcross(4),
165 0 : dlperp(4),
166 0 : dbi(3),
167 0 : dcross(3),
168 0 : dperp(3) {
169 : std::vector<AtomNumber> atoms;
170 0 : parseAtomList("BOX",atoms);
171 0 : if( atoms.size()!=4 ) {
172 0 : error("number of atoms should be equal to four");
173 : }
174 :
175 0 : log.printf(" boundaries for region are calculated based on positions of atoms : ");
176 0 : for(unsigned i=0; i<atoms.size(); ++i) {
177 0 : log.printf("%d ",atoms[i].serial() );
178 : }
179 0 : log.printf("\n");
180 :
181 0 : boxout=false;
182 0 : parseFlag("PRINT_BOX",boxout);
183 0 : if(boxout) {
184 : std::string boxfname;
185 0 : parse("FILE",boxfname);
186 0 : if(boxfname.length()==0) {
187 0 : error("no name for box file specified");
188 : }
189 : std::string unitname;
190 0 : parse("UNITS",unitname);
191 0 : if ( unitname.length()>0 ) {
192 0 : Units u;
193 0 : u.setLength(unitname);
194 0 : lenunit=getUnits().getLength()/u.getLength();
195 0 : } else {
196 : unitname="nm";
197 : }
198 0 : boxfile.link(*this);
199 0 : boxfile.open( boxfname );
200 0 : log.printf(" printing box coordinates on file named %s in %s \n",boxfname.c_str(), unitname.c_str() );
201 : }
202 :
203 0 : checkRead();
204 0 : requestAtoms(atoms);
205 0 : }
206 :
207 0 : VolumeTetrapore::~VolumeTetrapore() {
208 0 : }
209 :
210 0 : void VolumeTetrapore::setupRegions() {
211 : // Make some space for things
212 0 : Vector d1, d2, d3;
213 :
214 : // Retrieve the sigma value
215 0 : sigma=getSigma();
216 : // Set the position of the origin
217 0 : origin=getPosition(0);
218 :
219 : // Get two vectors
220 0 : d1 = pbcDistance(origin,getPosition(1));
221 0 : d2 = pbcDistance(origin,getPosition(2));
222 :
223 : // Find the vector connecting the origin to the top corner of
224 : // the subregion
225 0 : d3 = pbcDistance(origin,getPosition(3));
226 :
227 : // Create a set of unit vectors
228 0 : Vector bisector = d1 + d2;
229 0 : double bmod=bisector.modulo();
230 0 : bisector=bisector/bmod;
231 :
232 : // bi = d1 / d1l; len_bi=dotProduct( d3, bi );
233 0 : cross = crossProduct( d1, d2 );
234 0 : double crossmod=cross.modulo();
235 0 : cross = cross / crossmod;
236 0 : len_cross=dotProduct( d3, cross );
237 0 : Vector truep = crossProduct( cross, bisector );
238 :
239 : // These are our true vectors 45 degrees from bisector
240 0 : bi = cos(pi/4.0)*bisector + sin(pi/4.0)*truep;
241 0 : perp = cos(pi/4.0)*bisector - sin(pi/4.0)*truep;
242 :
243 : // And the lengths of the various parts average distance to opposite corners of tetetrahedron
244 0 : len_bi = dotProduct( d1, bi );
245 0 : double len_bi2 = dotProduct( d2, bi );
246 : unsigned lbi=1;
247 0 : if( len_bi2>len_bi ) {
248 0 : len_bi=len_bi2;
249 : lbi=2;
250 : }
251 0 : len_perp = dotProduct( d1, perp );
252 0 : double len_perp2 = dotProduct( d2, perp );
253 : unsigned lpi=1;
254 0 : if( len_perp2>len_perp ) {
255 0 : len_perp=len_perp2;
256 : lpi=2;
257 : }
258 0 : plumed_assert( lbi!=lpi );
259 :
260 0 : Tensor tcderiv;
261 0 : double cmod3=crossmod*crossmod*crossmod;
262 0 : Vector ucross=crossmod*cross;
263 0 : tcderiv.setCol( 0, crossProduct( d1, Vector(-1.0,0.0,0.0) ) + crossProduct( Vector(-1.0,0.0,0.0), d2 ) );
264 0 : tcderiv.setCol( 1, crossProduct( d1, Vector(0.0,-1.0,0.0) ) + crossProduct( Vector(0.0,-1.0,0.0), d2 ) );
265 0 : tcderiv.setCol( 2, crossProduct( d1, Vector(0.0,0.0,-1.0) ) + crossProduct( Vector(0.0,0.0,-1.0), d2 ) );
266 0 : dcross[0](0,0)=( tcderiv(0,0)/crossmod - ucross[0]*(ucross[0]*tcderiv(0,0) + ucross[1]*tcderiv(1,0) + ucross[2]*tcderiv(2,0))/cmod3 ); // dx/dx
267 0 : dcross[0](0,1)=( tcderiv(0,1)/crossmod - ucross[0]*(ucross[0]*tcderiv(0,1) + ucross[1]*tcderiv(1,1) + ucross[2]*tcderiv(2,1))/cmod3 ); // dx/dy
268 0 : dcross[0](0,2)=( tcderiv(0,2)/crossmod - ucross[0]*(ucross[0]*tcderiv(0,2) + ucross[1]*tcderiv(1,2) + ucross[2]*tcderiv(2,2))/cmod3 ); // dx/dz
269 0 : dcross[0](1,0)=( tcderiv(1,0)/crossmod - ucross[1]*(ucross[0]*tcderiv(0,0) + ucross[1]*tcderiv(1,0) + ucross[2]*tcderiv(2,0))/cmod3 ); // dy/dx
270 0 : dcross[0](1,1)=( tcderiv(1,1)/crossmod - ucross[1]*(ucross[0]*tcderiv(0,1) + ucross[1]*tcderiv(1,1) + ucross[2]*tcderiv(2,1))/cmod3 ); // dy/dy
271 0 : dcross[0](1,2)=( tcderiv(1,2)/crossmod - ucross[1]*(ucross[0]*tcderiv(0,2) + ucross[1]*tcderiv(1,2) + ucross[2]*tcderiv(2,2))/cmod3 ); // dy/dz
272 0 : dcross[0](2,0)=( tcderiv(2,0)/crossmod - ucross[2]*(ucross[0]*tcderiv(0,0) + ucross[1]*tcderiv(1,0) + ucross[2]*tcderiv(2,0))/cmod3 ); // dz/dx
273 0 : dcross[0](2,1)=( tcderiv(2,1)/crossmod - ucross[2]*(ucross[0]*tcderiv(0,1) + ucross[1]*tcderiv(1,1) + ucross[2]*tcderiv(2,1))/cmod3 ); // dz/dy
274 0 : dcross[0](2,2)=( tcderiv(2,2)/crossmod - ucross[2]*(ucross[0]*tcderiv(0,2) + ucross[1]*tcderiv(1,2) + ucross[2]*tcderiv(2,2))/cmod3 ); // dz/dz
275 :
276 0 : tcderiv.setCol( 0, crossProduct( Vector(1.0,0.0,0.0), d2 ) );
277 0 : tcderiv.setCol( 1, crossProduct( Vector(0.0,1.0,0.0), d2 ) );
278 0 : tcderiv.setCol( 2, crossProduct( Vector(0.0,0.0,1.0), d2 ) );
279 0 : dcross[1](0,0)=( tcderiv(0,0)/crossmod - ucross[0]*(ucross[0]*tcderiv(0,0) + ucross[1]*tcderiv(1,0) + ucross[2]*tcderiv(2,0))/cmod3 ); // dx/dx
280 0 : dcross[1](0,1)=( tcderiv(0,1)/crossmod - ucross[0]*(ucross[0]*tcderiv(0,1) + ucross[1]*tcderiv(1,1) + ucross[2]*tcderiv(2,1))/cmod3 ); // dx/dy
281 0 : dcross[1](0,2)=( tcderiv(0,2)/crossmod - ucross[0]*(ucross[0]*tcderiv(0,2) + ucross[1]*tcderiv(1,2) + ucross[2]*tcderiv(2,2))/cmod3 ); // dx/dz
282 0 : dcross[1](1,0)=( tcderiv(1,0)/crossmod - ucross[1]*(ucross[0]*tcderiv(0,0) + ucross[1]*tcderiv(1,0) + ucross[2]*tcderiv(2,0))/cmod3 ); // dy/dx
283 0 : dcross[1](1,1)=( tcderiv(1,1)/crossmod - ucross[1]*(ucross[0]*tcderiv(0,1) + ucross[1]*tcderiv(1,1) + ucross[2]*tcderiv(2,1))/cmod3 ); // dy/dy
284 0 : dcross[1](1,2)=( tcderiv(1,2)/crossmod - ucross[1]*(ucross[0]*tcderiv(0,2) + ucross[1]*tcderiv(1,2) + ucross[2]*tcderiv(2,2))/cmod3 ); // dy/dz
285 0 : dcross[1](2,0)=( tcderiv(2,0)/crossmod - ucross[2]*(ucross[0]*tcderiv(0,0) + ucross[1]*tcderiv(1,0) + ucross[2]*tcderiv(2,0))/cmod3 ); // dz/dx
286 0 : dcross[1](2,1)=( tcderiv(2,1)/crossmod - ucross[2]*(ucross[0]*tcderiv(0,1) + ucross[1]*tcderiv(1,1) + ucross[2]*tcderiv(2,1))/cmod3 ); // dz/dy
287 0 : dcross[1](2,2)=( tcderiv(2,2)/crossmod - ucross[2]*(ucross[0]*tcderiv(0,2) + ucross[1]*tcderiv(1,2) + ucross[2]*tcderiv(2,2))/cmod3 ); // dz/dz
288 :
289 0 : tcderiv.setCol( 0, crossProduct( d1, Vector(1.0,0.0,0.0) ) );
290 0 : tcderiv.setCol( 1, crossProduct( d1, Vector(0.0,1.0,0.0) ) );
291 0 : tcderiv.setCol( 2, crossProduct( d1, Vector(0.0,0.0,1.0) ) );
292 0 : dcross[2](0,0)=( tcderiv(0,0)/crossmod - ucross[0]*(ucross[0]*tcderiv(0,0) + ucross[1]*tcderiv(1,0) + ucross[2]*tcderiv(2,0))/cmod3 ); // dx/dx
293 0 : dcross[2](0,1)=( tcderiv(0,1)/crossmod - ucross[0]*(ucross[0]*tcderiv(0,1) + ucross[1]*tcderiv(1,1) + ucross[2]*tcderiv(2,1))/cmod3 ); // dx/dy
294 0 : dcross[2](0,2)=( tcderiv(0,2)/crossmod - ucross[0]*(ucross[0]*tcderiv(0,2) + ucross[1]*tcderiv(1,2) + ucross[2]*tcderiv(2,2))/cmod3 ); // dx/dz
295 0 : dcross[2](1,0)=( tcderiv(1,0)/crossmod - ucross[1]*(ucross[0]*tcderiv(0,0) + ucross[1]*tcderiv(1,0) + ucross[2]*tcderiv(2,0))/cmod3 ); // dy/dx
296 0 : dcross[2](1,1)=( tcderiv(1,1)/crossmod - ucross[1]*(ucross[0]*tcderiv(0,1) + ucross[1]*tcderiv(1,1) + ucross[2]*tcderiv(2,1))/cmod3 ); // dy/dy
297 0 : dcross[2](1,2)=( tcderiv(1,2)/crossmod - ucross[1]*(ucross[0]*tcderiv(0,2) + ucross[1]*tcderiv(1,2) + ucross[2]*tcderiv(2,2))/cmod3 ); // dy/dz
298 0 : dcross[2](2,0)=( tcderiv(2,0)/crossmod - ucross[2]*(ucross[0]*tcderiv(0,0) + ucross[1]*tcderiv(1,0) + ucross[2]*tcderiv(2,0))/cmod3 ); // dz/dx
299 0 : dcross[2](2,1)=( tcderiv(2,1)/crossmod - ucross[2]*(ucross[0]*tcderiv(0,1) + ucross[1]*tcderiv(1,1) + ucross[2]*tcderiv(2,1))/cmod3 ); // dz/dy
300 0 : dcross[2](2,2)=( tcderiv(2,2)/crossmod - ucross[2]*(ucross[0]*tcderiv(0,2) + ucross[1]*tcderiv(1,2) + ucross[2]*tcderiv(2,2))/cmod3 ); // dz/dz
301 :
302 0 : std::vector<Tensor> dbisector(3);
303 0 : double bmod3=bmod*bmod*bmod;
304 0 : Vector ubisector=bmod*bisector;
305 0 : dbisector[0](0,0)= -2.0/bmod + 2*ubisector[0]*ubisector[0]/bmod3;
306 0 : dbisector[0](0,1)= 2*ubisector[0]*ubisector[1]/bmod3;
307 0 : dbisector[0](0,2)= 2*ubisector[0]*ubisector[2]/bmod3;
308 0 : dbisector[0](1,0)= 2*ubisector[1]*ubisector[0]/bmod3;
309 0 : dbisector[0](1,1)= -2.0/bmod + 2*ubisector[1]*ubisector[1]/bmod3;
310 0 : dbisector[0](1,2)= 2*ubisector[1]*ubisector[2]/bmod3;
311 0 : dbisector[0](2,0)= 2*ubisector[2]*ubisector[0]/bmod3;
312 0 : dbisector[0](2,1)= 2*ubisector[2]*ubisector[1]/bmod3;
313 0 : dbisector[0](2,2)= -2.0/bmod + 2*ubisector[2]*ubisector[2]/bmod3;
314 :
315 0 : dbisector[1](0,0)= 1.0/bmod - ubisector[0]*ubisector[0]/bmod3;
316 0 : dbisector[1](0,1)= -ubisector[0]*ubisector[1]/bmod3;
317 0 : dbisector[1](0,2)= -ubisector[0]*ubisector[2]/bmod3;
318 0 : dbisector[1](1,0)= -ubisector[1]*ubisector[0]/bmod3;
319 0 : dbisector[1](1,1)= 1.0/bmod - ubisector[1]*ubisector[1]/bmod3;
320 0 : dbisector[1](1,2)= -ubisector[1]*ubisector[2]/bmod3;
321 0 : dbisector[1](2,0)= -ubisector[2]*ubisector[0]/bmod3;
322 0 : dbisector[1](2,1)= -ubisector[2]*ubisector[1]/bmod3;
323 0 : dbisector[1](2,2)=1.0/bmod - ubisector[2]*ubisector[2]/bmod3;
324 :
325 0 : dbisector[2](0,0)=1.0/bmod - ubisector[0]*ubisector[0]/bmod3;
326 0 : dbisector[2](0,1)= -ubisector[0]*ubisector[1]/bmod3;
327 0 : dbisector[2](0,2)= -ubisector[0]*ubisector[2]/bmod3;
328 0 : dbisector[2](1,0)= -ubisector[1]*ubisector[0]/bmod3;
329 0 : dbisector[2](1,1)=1.0/bmod - ubisector[1]*ubisector[1]/bmod3;
330 0 : dbisector[2](1,2)= -ubisector[1]*ubisector[2]/bmod3;
331 0 : dbisector[2](2,0)= -ubisector[2]*ubisector[0]/bmod3;
332 0 : dbisector[2](2,1)= -ubisector[2]*ubisector[1]/bmod3;
333 0 : dbisector[2](2,2)=1.0/bmod - ubisector[2]*ubisector[2]/bmod3;
334 :
335 0 : std::vector<Tensor> dtruep(3);
336 0 : dtruep[0].setCol( 0, ( crossProduct( dcross[0].getCol(0), bisector ) + crossProduct( cross, dbisector[0].getCol(0) ) ) );
337 0 : dtruep[0].setCol( 1, ( crossProduct( dcross[0].getCol(1), bisector ) + crossProduct( cross, dbisector[0].getCol(1) ) ) );
338 0 : dtruep[0].setCol( 2, ( crossProduct( dcross[0].getCol(2), bisector ) + crossProduct( cross, dbisector[0].getCol(2) ) ) );
339 :
340 0 : dtruep[1].setCol( 0, ( crossProduct( dcross[1].getCol(0), bisector ) + crossProduct( cross, dbisector[1].getCol(0) ) ) );
341 0 : dtruep[1].setCol( 1, ( crossProduct( dcross[1].getCol(1), bisector ) + crossProduct( cross, dbisector[1].getCol(1) ) ) );
342 0 : dtruep[1].setCol( 2, ( crossProduct( dcross[1].getCol(2), bisector ) + crossProduct( cross, dbisector[1].getCol(2) ) ) );
343 :
344 0 : dtruep[2].setCol( 0, ( crossProduct( dcross[2].getCol(0), bisector ) + crossProduct( cross, dbisector[2].getCol(0) ) ) );
345 0 : dtruep[2].setCol( 1, ( crossProduct( dcross[2].getCol(1), bisector ) + crossProduct( cross, dbisector[2].getCol(1) ) ) );
346 0 : dtruep[2].setCol( 2, ( crossProduct( dcross[2].getCol(2), bisector ) + crossProduct( cross, dbisector[2].getCol(2) ) ) );
347 :
348 : // Now convert these to the derivatives of the true axis
349 0 : for(unsigned i=0; i<3; ++i) {
350 0 : dbi[i] = cos(pi/4.0)*dbisector[i] + sin(pi/4.0)*dtruep[i];
351 0 : dperp[i] = cos(pi/4.0)*dbisector[i] - sin(pi/4.0)*dtruep[i];
352 : }
353 :
354 : // Ensure that all lengths are positive
355 0 : if( len_bi<0 ) {
356 0 : bi=-bi;
357 0 : len_bi=-len_bi;
358 0 : for(unsigned i=0; i<3; ++i) {
359 0 : dbi[i]*=-1.0;
360 : }
361 : }
362 0 : if( len_cross<0 ) {
363 0 : cross=-cross;
364 0 : len_cross=-len_cross;
365 0 : for(unsigned i=0; i<3; ++i) {
366 0 : dcross[i]*=-1.0;
367 : }
368 : }
369 0 : if( len_perp<0 ) {
370 0 : perp=-perp;
371 0 : len_perp=-len_perp;
372 0 : for(unsigned i=0; i<3; ++i) {
373 0 : dperp[i]*=-1.0;
374 : }
375 : }
376 0 : if( len_bi<=0 || len_cross<=0 || len_bi<=0 ) {
377 0 : plumed_merror("Invalid box coordinates");
378 : }
379 :
380 : // Now derivatives of lengths
381 0 : Tensor dd3( Tensor::identity() );
382 0 : Vector ddb2=d1;
383 0 : if( lbi==2 ) {
384 0 : ddb2=d2;
385 : }
386 0 : dlbi[1].zero();
387 0 : dlbi[2].zero();
388 0 : dlbi[3].zero();
389 0 : dlbi[0] = matmul(ddb2,dbi[0]) - matmul(bi,dd3);
390 0 : dlbi[lbi] = matmul(ddb2,dbi[lbi]) + matmul(bi,dd3); // Derivative wrt d1
391 :
392 0 : dlcross[0] = matmul(d3,dcross[0]) - matmul(cross,dd3);
393 0 : dlcross[1] = matmul(d3,dcross[1]);
394 0 : dlcross[2] = matmul(d3,dcross[2]);
395 0 : dlcross[3] = matmul(cross,dd3);
396 :
397 0 : ddb2=d1;
398 0 : if( lpi==2 ) {
399 0 : ddb2=d2;
400 : }
401 0 : dlperp[1].zero();
402 0 : dlperp[2].zero();
403 0 : dlperp[3].zero();
404 0 : dlperp[0] = matmul(ddb2,dperp[0]) - matmul( perp, dd3 );
405 0 : dlperp[lpi] = matmul(ddb2,dperp[lpi]) + matmul(perp, dd3);
406 :
407 : // Need to calculate the jacobian
408 0 : Tensor jacob;
409 0 : jacob(0,0)=bi[0];
410 0 : jacob(1,0)=bi[1];
411 0 : jacob(2,0)=bi[2];
412 0 : jacob(0,1)=cross[0];
413 0 : jacob(1,1)=cross[1];
414 0 : jacob(2,1)=cross[2];
415 0 : jacob(0,2)=perp[0];
416 0 : jacob(1,2)=perp[1];
417 0 : jacob(2,2)=perp[2];
418 0 : jacob_det = fabs( jacob.determinant() );
419 0 : }
420 :
421 0 : void VolumeTetrapore::update() {
422 0 : if(boxout) {
423 0 : boxfile.printf("%d\n",8);
424 0 : const Tensor & t(getPbc().getBox());
425 0 : if(getPbc().isOrthorombic()) {
426 0 : boxfile.printf(" %f %f %f\n",lenunit*t(0,0),lenunit*t(1,1),lenunit*t(2,2));
427 : } else {
428 0 : boxfile.printf(" %f %f %f %f %f %f %f %f %f\n",
429 0 : lenunit*t(0,0),lenunit*t(0,1),lenunit*t(0,2),
430 0 : lenunit*t(1,0),lenunit*t(1,1),lenunit*t(1,2),
431 0 : lenunit*t(2,0),lenunit*t(2,1),lenunit*t(2,2)
432 : );
433 : }
434 0 : boxfile.printf("AR %f %f %f \n",lenunit*origin[0],lenunit*origin[1],lenunit*origin[2]);
435 0 : Vector ut, vt, wt;
436 0 : ut = origin + len_bi*bi;
437 0 : vt = origin + len_cross*cross;
438 0 : wt = origin + len_perp*perp;
439 0 : boxfile.printf("AR %f %f %f \n",lenunit*(ut[0]), lenunit*(ut[1]), lenunit*(ut[2]) );
440 0 : boxfile.printf("AR %f %f %f \n",lenunit*(vt[0]), lenunit*(vt[1]), lenunit*(vt[2]) );
441 0 : boxfile.printf("AR %f %f %f \n",lenunit*(wt[0]), lenunit*(wt[1]), lenunit*(wt[2]) );
442 0 : boxfile.printf("AR %f %f %f \n",lenunit*(vt[0]+len_bi*bi[0]),
443 0 : lenunit*(vt[1]+len_bi*bi[1]),
444 0 : lenunit*(vt[2]+len_bi*bi[2]) );
445 0 : boxfile.printf("AR %f %f %f \n",lenunit*(ut[0]+len_perp*perp[0]),
446 0 : lenunit*(ut[1]+len_perp*perp[1]),
447 0 : lenunit*(ut[2]+len_perp*perp[2]) );
448 0 : boxfile.printf("AR %f %f %f \n",lenunit*(vt[0]+len_perp*perp[0]),
449 0 : lenunit*(vt[1]+len_perp*perp[1]),
450 0 : lenunit*(vt[2]+len_perp*perp[2]) );
451 0 : boxfile.printf("AR %f %f %f \n",lenunit*(vt[0]+len_perp*perp[0]+len_bi*bi[0]),
452 0 : lenunit*(vt[1]+len_perp*perp[1]+len_bi*bi[1]),
453 0 : lenunit*(vt[2]+len_perp*perp[2]+len_bi*bi[2]) );
454 : }
455 0 : }
456 :
457 0 : double VolumeTetrapore::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& rderiv ) const {
458 : // Setup the histogram bead
459 0 : HistogramBead bead;
460 : bead.isNotPeriodic();
461 0 : bead.setKernelType( getKernelType() );
462 :
463 : // Calculate distance of atom from origin of new coordinate frame
464 0 : Vector datom=pbcDistance( origin, cpos );
465 : double ucontr, uder, vcontr, vder, wcontr, wder;
466 :
467 : // Calculate contribution from integral along bi
468 0 : bead.set( 0, len_bi, sigma );
469 0 : double upos=dotProduct( datom, bi );
470 0 : ucontr=bead.calculate( upos, uder );
471 0 : double udlen=bead.uboundDerivative( upos );
472 0 : double uder2 = bead.lboundDerivative( upos ) - udlen;
473 :
474 : // Calculate contribution from integral along cross
475 0 : bead.set( 0, len_cross, sigma );
476 0 : double vpos=dotProduct( datom, cross );
477 0 : vcontr=bead.calculate( vpos, vder );
478 0 : double vdlen=bead.uboundDerivative( vpos );
479 0 : double vder2 = bead.lboundDerivative( vpos ) - vdlen;
480 :
481 : // Calculate contribution from integral along perp
482 0 : bead.set( 0, len_perp, sigma );
483 0 : double wpos=dotProduct( datom, perp );
484 0 : wcontr=bead.calculate( wpos, wder );
485 0 : double wdlen=bead.uboundDerivative( wpos );
486 0 : double wder2 = bead.lboundDerivative( wpos ) - wdlen;
487 :
488 0 : Vector dfd;
489 0 : dfd[0]=uder*vcontr*wcontr;
490 0 : dfd[1]=ucontr*vder*wcontr;
491 0 : dfd[2]=ucontr*vcontr*wder;
492 0 : derivatives[0] = (dfd[0]*bi[0]+dfd[1]*cross[0]+dfd[2]*perp[0]);
493 0 : derivatives[1] = (dfd[0]*bi[1]+dfd[1]*cross[1]+dfd[2]*perp[1]);
494 0 : derivatives[2] = (dfd[0]*bi[2]+dfd[1]*cross[2]+dfd[2]*perp[2]);
495 0 : double tot = ucontr*vcontr*wcontr*jacob_det;
496 :
497 : // Add reference atom derivatives
498 0 : dfd[0]=uder2*vcontr*wcontr;
499 0 : dfd[1]=ucontr*vder2*wcontr;
500 0 : dfd[2]=ucontr*vcontr*wder2;
501 0 : Vector dfld;
502 0 : dfld[0]=udlen*vcontr*wcontr;
503 0 : dfld[1]=ucontr*vdlen*wcontr;
504 0 : dfld[2]=ucontr*vcontr*wdlen;
505 0 : rderiv[0] = dfd[0]*matmul(datom,dbi[0]) + dfd[1]*matmul(datom,dcross[0]) + dfd[2]*matmul(datom,dperp[0]) +
506 0 : dfld[0]*dlbi[0] + dfld[1]*dlcross[0] + dfld[2]*dlperp[0] - derivatives;
507 0 : rderiv[1] = dfd[0]*matmul(datom,dbi[1]) + dfd[1]*matmul(datom,dcross[1]) + dfd[2]*matmul(datom,dperp[1]) +
508 0 : dfld[0]*dlbi[1] + dfld[1]*dlcross[1] + dfld[2]*dlperp[1];
509 0 : rderiv[2] = dfd[0]*matmul(datom,dbi[2]) + dfd[1]*matmul(datom,dcross[2]) + dfd[2]*matmul(datom,dperp[2]) +
510 0 : dfld[0]*dlbi[2] + dfld[1]*dlcross[2] + dfld[2]*dlperp[2];
511 0 : rderiv[3] = dfld[0]*dlbi[3] + dfld[1]*dlcross[3] + dfld[2]*dlperp[3];
512 :
513 0 : vir.zero();
514 0 : vir-=Tensor( cpos,derivatives );
515 0 : for(unsigned i=0; i<4; ++i) {
516 0 : vir -= Tensor( getPosition(i), rderiv[i] );
517 : }
518 :
519 0 : return tot;
520 : }
521 :
522 : }
523 : }
|
