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