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 : double jacob_det{0.0};
141 : double len_bi{0.0}, len_cross{0.0}, len_perp{0.0}, sigma{0.0};
142 : Vector bi, cross, perp;
143 : HistogramBead::KernelType kerneltype{HistogramBead::KernelType::gaussian};
144 : std::vector<Vector> dlbi, dlcross, dlperp;
145 : std::vector<Tensor> dbi, dcross, dperp;
146 : static void registerKeywords( Keywords& keys );
147 2 : VolumeTetrapore() : dlbi(4), dlcross(4), dlperp(4), dbi(3), dcross(3), dperp(3) {}
148 : void setupRegions( ActionVolume<VolumeTetrapore>* action,
149 : const Pbc& pbc,
150 : const std::vector<Vector>& positions );
151 : void parseInput( ActionVolume<VolumeTetrapore>* action );
152 : static void parseAtoms( ActionVolume<VolumeTetrapore>* action,
153 : std::vector<AtomNumber>& atoms );
154 1 : VolumeTetrapore& operator=( const VolumeTetrapore& m ) {
155 1 : jacob_det=m.jacob_det;
156 1 : len_bi=m.len_bi;
157 1 : len_cross=m.len_cross;
158 1 : len_perp=m.len_perp;
159 1 : sigma=m.sigma;
160 1 : dlbi.resize(4);
161 1 : dlcross.resize(4);
162 1 : dlperp.resize(4);
163 1 : dbi.resize(3);
164 1 : dcross.resize(3);
165 1 : dperp.resize(3);
166 1 : kerneltype=m.kerneltype;
167 1 : return *this;
168 : }
169 : static void calculateNumberInside( const VolumeInput& input, const VolumeTetrapore& actioninput, VolumeOutput& output );
170 : };
171 :
172 7 : void VolumeTetrapore::registerKeywords( Keywords& keys ) {
173 14 : keys.setDisplayName("TETRAHEDRALPORE");
174 7 : keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
175 7 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
176 7 : keys.add("atoms","BOX","the positions of four atoms that define spatial extent of the cavity");
177 7 : }
178 :
179 1 : void VolumeTetrapore::parseInput( ActionVolume<VolumeTetrapore>* action ) {
180 2 : action->parse("SIGMA",sigma);
181 : std::string mykerneltype;
182 1 : action->parse("KERNEL",mykerneltype);
183 1 : kerneltype=HistogramBead::getKernelType(mykerneltype);
184 1 : }
185 :
186 1 : void VolumeTetrapore::parseAtoms( ActionVolume<VolumeTetrapore>* action, std::vector<AtomNumber>& atoms ) {
187 2 : action->parseAtomList("BOX",atoms);
188 1 : if( atoms.size()!=4 ) {
189 0 : action->error("number of atoms should be equal to four");
190 : }
191 :
192 1 : action->log.printf(" boundaries for region are calculated based on positions of atoms : ");
193 5 : for(unsigned i=0; i<atoms.size(); ++i) {
194 4 : action->log.printf("%d ",atoms[i].serial() );
195 : }
196 1 : action->log.printf("\n");
197 1 : }
198 :
199 : typedef ActionVolume<VolumeTetrapore> VolTet;
200 : PLUMED_REGISTER_ACTION(VolTet,"TETRAHEDRALPORE_CALC")
201 : char glob_tetrapore[] = "TETRAHEDRALPORE";
202 : typedef VolumeShortcut<glob_tetrapore> VolumeTetraporeShortcut;
203 : PLUMED_REGISTER_ACTION(VolumeTetraporeShortcut,"TETRAHEDRALPORE")
204 :
205 1560 : void VolumeTetrapore::setupRegions( ActionVolume<VolumeTetrapore>* action, const Pbc& pbc, const std::vector<Vector>& positions ) {
206 : // Make some space for things
207 : Vector d1, d2, d3;
208 :
209 : // Set the position of the origin
210 1560 : Vector origin=positions[0];
211 :
212 : // Get two vectors
213 1560 : d1 = pbc.distance(origin,positions[1]);
214 1560 : d2 = pbc.distance(origin,positions[2]);
215 :
216 : // Find the vector connecting the origin to the top corner of
217 : // the subregion
218 1560 : d3 = pbc.distance(origin,positions[3]);
219 :
220 : // Create a set of unit vectors
221 : Vector bisector = d1 + d2;
222 1560 : double bmod=bisector.modulo();
223 1560 : bisector=bisector/bmod;
224 :
225 : // bi = d1 / d1l; len_bi=dotProduct( d3, bi );
226 1560 : cross = crossProduct( d1, d2 );
227 1560 : double crossmod=cross.modulo();
228 1560 : cross = cross / crossmod;
229 1560 : len_cross=dotProduct( d3, cross );
230 1560 : Vector truep = crossProduct( cross, bisector );
231 :
232 : // These are our true vectors 45 degrees from bisector
233 1560 : bi = cos(pi/4.0)*bisector + sin(pi/4.0)*truep;
234 1560 : perp = cos(pi/4.0)*bisector - sin(pi/4.0)*truep;
235 :
236 : // And the lengths of the various parts average distance to opposite corners of tetetrahedron
237 1560 : len_bi = dotProduct( d1, bi );
238 : double len_bi2 = dotProduct( d2, bi );
239 : unsigned lbi=1;
240 1560 : if( len_bi2>len_bi ) {
241 1560 : len_bi=len_bi2;
242 : lbi=2;
243 : }
244 1560 : len_perp = dotProduct( d1, perp );
245 : double len_perp2 = dotProduct( d2, perp );
246 : unsigned lpi=1;
247 1560 : if( len_perp2>len_perp ) {
248 0 : len_perp=len_perp2;
249 : lpi=2;
250 : }
251 1560 : plumed_assert( lbi!=lpi );
252 :
253 : Tensor tcderiv;
254 1560 : double cmod3=crossmod*crossmod*crossmod;
255 : Vector ucross=crossmod*cross;
256 3120 : tcderiv.setCol( 0, crossProduct( d1, Vector(-1.0,0.0,0.0) ) + crossProduct( Vector(-1.0,0.0,0.0), d2 ) );
257 3120 : tcderiv.setCol( 1, crossProduct( d1, Vector(0.0,-1.0,0.0) ) + crossProduct( Vector(0.0,-1.0,0.0), d2 ) );
258 1560 : tcderiv.setCol( 2, crossProduct( d1, Vector(0.0,0.0,-1.0) ) + crossProduct( Vector(0.0,0.0,-1.0), d2 ) );
259 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
260 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
261 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
262 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
263 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
264 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
265 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
266 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
267 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
268 :
269 3120 : tcderiv.setCol( 0, crossProduct( Vector(1.0,0.0,0.0), d2 ) );
270 3120 : tcderiv.setCol( 1, crossProduct( Vector(0.0,1.0,0.0), d2 ) );
271 1560 : tcderiv.setCol( 2, crossProduct( Vector(0.0,0.0,1.0), d2 ) );
272 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
273 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
274 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
275 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
276 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
277 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
278 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
279 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
280 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
281 :
282 3120 : tcderiv.setCol( 0, crossProduct( d1, Vector(1.0,0.0,0.0) ) );
283 3120 : tcderiv.setCol( 1, crossProduct( d1, Vector(0.0,1.0,0.0) ) );
284 1560 : tcderiv.setCol( 2, crossProduct( d1, Vector(0.0,0.0,1.0) ) );
285 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
286 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
287 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
288 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
289 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
290 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
291 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
292 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
293 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
294 :
295 1560 : std::vector<Tensor> dbisector(3);
296 1560 : double bmod3=bmod*bmod*bmod;
297 : Vector ubisector=bmod*bisector;
298 1560 : dbisector[0](0,0)= -2.0/bmod + 2*ubisector[0]*ubisector[0]/bmod3;
299 1560 : dbisector[0](0,1)= 2*ubisector[0]*ubisector[1]/bmod3;
300 1560 : dbisector[0](0,2)= 2*ubisector[0]*ubisector[2]/bmod3;
301 1560 : dbisector[0](1,0)= 2*ubisector[1]*ubisector[0]/bmod3;
302 1560 : dbisector[0](1,1)= -2.0/bmod + 2*ubisector[1]*ubisector[1]/bmod3;
303 1560 : dbisector[0](1,2)= 2*ubisector[1]*ubisector[2]/bmod3;
304 1560 : dbisector[0](2,0)= 2*ubisector[2]*ubisector[0]/bmod3;
305 1560 : dbisector[0](2,1)= 2*ubisector[2]*ubisector[1]/bmod3;
306 1560 : dbisector[0](2,2)= -2.0/bmod + 2*ubisector[2]*ubisector[2]/bmod3;
307 :
308 1560 : dbisector[1](0,0)= 1.0/bmod - ubisector[0]*ubisector[0]/bmod3;
309 1560 : dbisector[1](0,1)= -ubisector[0]*ubisector[1]/bmod3;
310 1560 : dbisector[1](0,2)= -ubisector[0]*ubisector[2]/bmod3;
311 1560 : dbisector[1](1,0)= -ubisector[1]*ubisector[0]/bmod3;
312 1560 : dbisector[1](1,1)= 1.0/bmod - ubisector[1]*ubisector[1]/bmod3;
313 1560 : dbisector[1](1,2)= -ubisector[1]*ubisector[2]/bmod3;
314 1560 : dbisector[1](2,0)= -ubisector[2]*ubisector[0]/bmod3;
315 1560 : dbisector[1](2,1)= -ubisector[2]*ubisector[1]/bmod3;
316 1560 : dbisector[1](2,2)=1.0/bmod - ubisector[2]*ubisector[2]/bmod3;
317 :
318 1560 : dbisector[2](0,0)=1.0/bmod - ubisector[0]*ubisector[0]/bmod3;
319 1560 : dbisector[2](0,1)= -ubisector[0]*ubisector[1]/bmod3;
320 1560 : dbisector[2](0,2)= -ubisector[0]*ubisector[2]/bmod3;
321 1560 : dbisector[2](1,0)= -ubisector[1]*ubisector[0]/bmod3;
322 1560 : dbisector[2](1,1)=1.0/bmod - ubisector[1]*ubisector[1]/bmod3;
323 1560 : dbisector[2](1,2)= -ubisector[1]*ubisector[2]/bmod3;
324 1560 : dbisector[2](2,0)= -ubisector[2]*ubisector[0]/bmod3;
325 1560 : dbisector[2](2,1)= -ubisector[2]*ubisector[1]/bmod3;
326 1560 : dbisector[2](2,2)=1.0/bmod - ubisector[2]*ubisector[2]/bmod3;
327 :
328 1560 : std::vector<Tensor> dtruep(3);
329 4680 : dtruep[0].setCol( 0, ( crossProduct( dcross[0].getCol(0), bisector ) + crossProduct( cross, dbisector[0].getCol(0) ) ) );
330 4680 : dtruep[0].setCol( 1, ( crossProduct( dcross[0].getCol(1), bisector ) + crossProduct( cross, dbisector[0].getCol(1) ) ) );
331 4680 : dtruep[0].setCol( 2, ( crossProduct( dcross[0].getCol(2), bisector ) + crossProduct( cross, dbisector[0].getCol(2) ) ) );
332 :
333 4680 : dtruep[1].setCol( 0, ( crossProduct( dcross[1].getCol(0), bisector ) + crossProduct( cross, dbisector[1].getCol(0) ) ) );
334 4680 : dtruep[1].setCol( 1, ( crossProduct( dcross[1].getCol(1), bisector ) + crossProduct( cross, dbisector[1].getCol(1) ) ) );
335 4680 : dtruep[1].setCol( 2, ( crossProduct( dcross[1].getCol(2), bisector ) + crossProduct( cross, dbisector[1].getCol(2) ) ) );
336 :
337 4680 : dtruep[2].setCol( 0, ( crossProduct( dcross[2].getCol(0), bisector ) + crossProduct( cross, dbisector[2].getCol(0) ) ) );
338 4680 : dtruep[2].setCol( 1, ( crossProduct( dcross[2].getCol(1), bisector ) + crossProduct( cross, dbisector[2].getCol(1) ) ) );
339 3120 : dtruep[2].setCol( 2, ( crossProduct( dcross[2].getCol(2), bisector ) + crossProduct( cross, dbisector[2].getCol(2) ) ) );
340 :
341 : // Now convert these to the derivatives of the true axis
342 6240 : for(unsigned i=0; i<3; ++i) {
343 9360 : dbi[i] = cos(pi/4.0)*dbisector[i] + sin(pi/4.0)*dtruep[i];
344 4680 : dperp[i] = cos(pi/4.0)*dbisector[i] - sin(pi/4.0)*dtruep[i];
345 : }
346 :
347 : // Ensure that all lengths are positive
348 1560 : if( len_bi<0 ) {
349 0 : bi=-bi;
350 0 : len_bi=-len_bi;
351 0 : for(unsigned i=0; i<3; ++i) {
352 0 : dbi[i]*=-1.0;
353 : }
354 : }
355 1560 : if( len_cross<0 ) {
356 0 : cross=-cross;
357 0 : len_cross=-len_cross;
358 0 : for(unsigned i=0; i<3; ++i) {
359 0 : dcross[i]*=-1.0;
360 : }
361 : }
362 1560 : if( len_perp<0 ) {
363 0 : perp=-perp;
364 0 : len_perp=-len_perp;
365 0 : for(unsigned i=0; i<3; ++i) {
366 0 : dperp[i]*=-1.0;
367 : }
368 : }
369 1560 : if( len_bi<=0 || len_cross<=0 || len_bi<=0 ) {
370 0 : plumed_merror("Invalid box coordinates");
371 : }
372 :
373 : // Now derivatives of lengths
374 : Tensor dd3( Tensor::identity() );
375 1560 : Vector ddb2=d1;
376 1560 : if( lbi==2 ) {
377 1560 : ddb2=d2;
378 : }
379 : dlbi[1].zero();
380 : dlbi[2].zero();
381 : dlbi[3].zero();
382 1560 : dlbi[0] = matmul(ddb2,dbi[0]) - matmul(bi,dd3);
383 1560 : dlbi[lbi] = matmul(ddb2,dbi[lbi]) + matmul(bi,dd3); // Derivative wrt d1
384 :
385 1560 : dlcross[0] = matmul(d3,dcross[0]) - matmul(cross,dd3);
386 1560 : dlcross[1] = matmul(d3,dcross[1]);
387 1560 : dlcross[2] = matmul(d3,dcross[2]);
388 1560 : dlcross[3] = matmul(cross,dd3);
389 :
390 1560 : ddb2=d1;
391 1560 : if( lpi==2 ) {
392 0 : ddb2=d2;
393 : }
394 : dlperp[1].zero();
395 : dlperp[2].zero();
396 : dlperp[3].zero();
397 1560 : dlperp[0] = matmul(ddb2,dperp[0]) - matmul( perp, dd3 );
398 1560 : dlperp[lpi] = matmul(ddb2,dperp[lpi]) + matmul(perp, dd3);
399 :
400 : // Need to calculate the jacobian
401 : Tensor jacob;
402 1560 : jacob(0,0)=bi[0];
403 1560 : jacob(1,0)=bi[1];
404 1560 : jacob(2,0)=bi[2];
405 1560 : jacob(0,1)=cross[0];
406 1560 : jacob(1,1)=cross[1];
407 1560 : jacob(2,1)=cross[2];
408 1560 : jacob(0,2)=perp[0];
409 1560 : jacob(1,2)=perp[1];
410 1560 : jacob(2,2)=perp[2];
411 1560 : jacob_det = fabs( jacob.determinant() );
412 1560 : }
413 :
414 1560 : void VolumeTetrapore::calculateNumberInside( const VolumeInput& input,
415 : const VolumeTetrapore& actioninput,
416 : VolumeOutput& output ) {
417 : // Calculate distance of atom from origin of new coordinate frame
418 1560 : Vector datom=input.pbc.distance(
419 1560 : Vector(input.refpos[0][0],input.refpos[0][1],input.refpos[0][2]),
420 1560 : Vector(input.cpos[0],input.cpos[1],input.cpos[2]) );
421 : double ucontr, uder, vcontr, vder, wcontr, wder;
422 :
423 : // Calculate contribution from integral along bi
424 : // Setup the histogram bead
425 1560 : HistogramBead bead( actioninput.kerneltype,
426 1560 : 0, actioninput.len_bi, actioninput.sigma );
427 : double upos=dotProduct( datom, actioninput.bi );
428 1560 : ucontr=bead.calculate( upos, uder );
429 1560 : double udlen=bead.uboundDerivative( upos );
430 1560 : double uder2 = bead.lboundDerivative( upos ) - udlen;
431 :
432 : // Calculate contribution from integral along cross
433 1560 : bead.set( 0, actioninput.len_cross, actioninput.sigma );
434 : double vpos=dotProduct( datom, actioninput.cross );
435 1560 : vcontr=bead.calculate( vpos, vder );
436 1560 : double vdlen=bead.uboundDerivative( vpos );
437 1560 : double vder2 = bead.lboundDerivative( vpos ) - vdlen;
438 :
439 : // Calculate contribution from integral along perp
440 1560 : bead.set( 0, actioninput.len_perp, actioninput.sigma );
441 : double wpos=dotProduct( datom, actioninput.perp );
442 1560 : wcontr=bead.calculate( wpos, wder );
443 1560 : double wdlen=bead.uboundDerivative( wpos );
444 1560 : double wder2 = bead.lboundDerivative( wpos ) - wdlen;
445 :
446 : Vector dfd;
447 1560 : dfd[0]=uder*vcontr*wcontr;
448 1560 : dfd[1]=ucontr*vder*wcontr;
449 1560 : dfd[2]=ucontr*vcontr*wder;
450 1560 : output.derivatives[0] = (dfd[0]*actioninput.bi[0]+dfd[1]*actioninput.cross[0]+dfd[2]*actioninput.perp[0]);
451 1560 : output.derivatives[1] = (dfd[0]*actioninput.bi[1]+dfd[1]*actioninput.cross[1]+dfd[2]*actioninput.perp[1]);
452 1560 : output.derivatives[2] = (dfd[0]*actioninput.bi[2]+dfd[1]*actioninput.cross[2]+dfd[2]*actioninput.perp[2]);
453 1560 : output.values[0] = ucontr*vcontr*wcontr*actioninput.jacob_det;
454 :
455 : // Add reference atom derivatives
456 1560 : dfd[0]=uder2*vcontr*wcontr;
457 1560 : dfd[1]=ucontr*vder2*wcontr;
458 1560 : dfd[2]=ucontr*vcontr*wder2;
459 : Vector dfld;
460 1560 : dfld[0]=udlen*vcontr*wcontr;
461 1560 : dfld[1]=ucontr*vdlen*wcontr;
462 1560 : dfld[2]=ucontr*vcontr*wdlen;
463 1560 : output.refders[0] = dfd[0]*matmul(datom,actioninput.dbi[0])
464 1560 : + dfd[1]*matmul(datom,actioninput.dcross[0])
465 1560 : + dfd[2]*matmul(datom,actioninput.dperp[0])
466 : + dfld[0]*actioninput.dlbi[0]
467 : + dfld[1]*actioninput.dlcross[0]
468 : + dfld[2]*actioninput.dlperp[0]
469 1560 : - Vector(output.derivatives[0],output.derivatives[1],output.derivatives[2]);
470 1560 : output.refders[1] = dfd[0]*matmul(datom,actioninput.dbi[1])
471 1560 : + dfd[1]*matmul(datom,actioninput.dcross[1])
472 1560 : + dfd[2]*matmul(datom,actioninput.dperp[1])
473 : + dfld[0]*actioninput.dlbi[1]
474 : + dfld[1]*actioninput.dlcross[1]
475 : + dfld[2]*actioninput.dlperp[1];
476 1560 : output.refders[2] = dfd[0]*matmul(datom,actioninput.dbi[2])
477 1560 : + dfd[1]*matmul(datom,actioninput.dcross[2])
478 1560 : + dfd[2]*matmul(datom,actioninput.dperp[2])
479 : + dfld[0]*actioninput.dlbi[2]
480 : + dfld[1]*actioninput.dlcross[2]
481 : + dfld[2]*actioninput.dlperp[2];
482 : output.refders[3] = dfld[0]*actioninput.dlbi[3]
483 : + dfld[1]*actioninput.dlcross[3]
484 : + dfld[2]*actioninput.dlperp[3];
485 :
486 : Tensor vir;
487 1560 : vir=-Tensor( Vector(input.cpos[0],input.cpos[1],input.cpos[2]),
488 1560 : Vector(output.derivatives[0],output.derivatives[1],output.derivatives[2]) );
489 7800 : for(unsigned i=0; i<4; ++i) {
490 6240 : vir -= Tensor( Vector(input.refpos[i][0],input.refpos[i][1],input.refpos[i][2]),
491 12480 : Vector(output.refders[i][0],output.refders[i][1],output.refders[i][2]) );
492 : }
493 1560 : output.virial.set( 0, vir );
494 1560 : }
495 :
496 : }
497 : }
|
