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