Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 "tools/HistogramBead.h"
28 : #include "VolumeShortcut.h"
29 :
30 : //+PLUMEDOC VOLUMES CAVITY
31 : /*
32 : This quantity can be used to calculate functions of the distribution of collective variables for the atoms that lie in a box defined by the positions of four atoms.
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 : In this expression $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. The first of these unit vectors points from the first reference atom to the second.
64 : The second is then the normal to the plane containing atoms 1,2 and 3 and the the third is the unit vector orthogonal to
65 : these first two vectors. $(x_o,y_o,z_o)$, meanwhile, specifies the position of the first reference atom.
66 :
67 : In the first expression above $K$ is one of the kernel functions described in the documentation for the action [BETWEEN](BETWEEN.md)
68 : and $\sigma$ is a bandwidth parameter. Furthermore, The vector connecting atom 1 to atom 4 is used to define the extent of the box in
69 : each of the $u$, $v$ and $w$ directions. This vector connecting atom 1 to atom 4 is projected onto the three unit vectors
70 : described above and the resulting projections determine the $u'$, $v'$ and $w'$ parameters in the above expression.
71 :
72 : 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
73 : in an ion channel in a protein. The extent of the channel is calculated from the positions of atoms 1, 4, 5 and 11.
74 :
75 : ```plumed
76 : cav: CAVITY ...
77 : ATOMS=20-500 BOX=1,4,5,11
78 : SIGMA=0.1 KERNEL=gaussian
79 : ...
80 : s: SUM ARG=cav PERIODIC=NO
81 : PRINT ARG=s FILE=colvar
82 : ```
83 :
84 : If you want to calculate the number of atoms that are not inthe protein chanel you can use the `OUTSIDE` flag as shown below:
85 :
86 : ```plumed
87 : cav: CAVITY ...
88 : ATOMS=20-500 BOX=1,4,5,11
89 : SIGMA=0.1 KERNEL=gaussian
90 : OUTSIDE
91 : ...
92 : s: SUM ARG=cav PERIODIC=NO
93 : PRINT ARG=s FILE=colvar
94 : ```
95 :
96 : By contrast the following command tells plumed to calculate the coordination numbers (with other water molecules) for the water
97 : molecules in the protein channel described above. The average coordination number and the number of coordination
98 : numbers more than 4 is then calculated for those molecules that are in the region of interest.
99 :
100 : ```plumed
101 : # Calculate the atoms that are in the cavity
102 : cav: CAVITY ATOMS=20-500 BOX=1,4,5,11 SIGMA=0.1
103 : # Calculate the coordination numbers of all the atoms
104 : d1: COORDINATIONNUMBER SPECIES=20-500 SWITCH={RATIONAL R_0=0.1}
105 : # Do atoms have a coordination number that is greater than 4
106 : fd1: MORE_THAN ARG=d1 SWITCH={RATIONAL R_0=4}
107 : # Calculate the mean coordination number in the channel
108 : nnn: CUSTOM ARG=cav,d1 FUNC=x*y PERIODIC=NO
109 : numer: SUM ARG=nnn PERIODIC=NO
110 : denom: SUM ARG=cav PERIODIC=NO
111 : mean: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
112 : # Calculate the number of atoms that are in the channel and that
113 : # have a coordination number that is greater than 4
114 : sss: CUSTOM ARG=fd1,cav FUNC=x*y PERIODIC=NO
115 : mt: SUM ARG=sss PERIODIC=NO
116 : # And print these two quantities to a file
117 : PRINT ARG=mean,mt FILE=colvar
118 : ```
119 :
120 : !!! note ""
121 :
122 : As with [AROUND](AROUND.md) earlier version of PLUMED used a different syntax for doing these types of calculations, which can
123 : still be used with this new version of the command. We strongly recommend using the newer syntax but if you are interested in the
124 : old syntax you can find more information in the old syntax section of the documentation for [AROUND](AROUND.md).
125 :
126 : */
127 : //+ENDPLUMEDOC
128 :
129 : namespace PLMD {
130 : namespace volumes {
131 :
132 : class VolumeCavity {
133 : public:
134 : double jacob_det;
135 : double len_bi, len_cross, len_perp, sigma;
136 : Vector bi, cross, perp;
137 : HistogramBead::KernelType kerneltype{HistogramBead::KernelType::gaussian};
138 : std::vector<Vector> dlbi, dlcross, dlperp;
139 : std::vector<Tensor> dbi, dcross, dperp;
140 : static void registerKeywords( Keywords& keys );
141 2 : VolumeCavity() : jacob_det(0), len_bi(0), len_cross(0), len_perp(0), sigma(0), dlbi(4), dlcross(4), dlperp(4), dbi(3), dcross(3), dperp(3) {}
142 : void setupRegions( ActionVolume<VolumeCavity>* action, const Pbc& pbc, const std::vector<Vector>& positions );
143 : void parseInput( ActionVolume<VolumeCavity>* action );
144 : static void parseAtoms( ActionVolume<VolumeCavity>* action, std::vector<AtomNumber>& atoms );
145 1 : VolumeCavity& operator=( const VolumeCavity& m ) {
146 1 : jacob_det=m.jacob_det;
147 1 : len_bi=m.len_bi;
148 1 : len_cross=m.len_cross;
149 1 : len_perp=m.len_perp;
150 1 : sigma=m.sigma;
151 1 : dlbi.resize(4);
152 1 : dlcross.resize(4);
153 1 : dlperp.resize(4);
154 1 : dbi.resize(3);
155 1 : dcross.resize(3);
156 1 : dperp.resize(3);
157 1 : kerneltype=m.kerneltype;
158 1 : return *this;
159 : }
160 : static void calculateNumberInside( const VolumeInput& input, const VolumeCavity& actioninput, VolumeOutput& output );
161 : };
162 :
163 : typedef ActionVolume<VolumeCavity> VolCav;
164 : PLUMED_REGISTER_ACTION(VolCav,"CAVITY_CALC")
165 : char glob_cavity[] = "CAVITY";
166 : typedef VolumeShortcut<glob_cavity> VolumeCavityShortcut;
167 : PLUMED_REGISTER_ACTION(VolumeCavityShortcut,"CAVITY")
168 :
169 7 : void VolumeCavity::registerKeywords( Keywords& keys ) {
170 14 : keys.setDisplayName("CAVITY");
171 7 : keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
172 7 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
173 7 : keys.add("atoms","BOX","the positions of four atoms that define spatial extent of the cavity");
174 7 : }
175 :
176 1 : void VolumeCavity::parseInput( ActionVolume<VolumeCavity>* action ) {
177 2 : action->parse("SIGMA",sigma);
178 : std::string mykerneltype;
179 1 : action->parse("KERNEL",mykerneltype);
180 1 : kerneltype=HistogramBead::getKernelType(mykerneltype);
181 1 : action->log.printf(" using %s kernels with a bandwidth of %d \n", mykerneltype.c_str(), sigma );
182 1 : }
183 :
184 1 : void VolumeCavity::parseAtoms( ActionVolume<VolumeCavity>* action, std::vector<AtomNumber>& atoms ) {
185 2 : action->parseAtomList("BOX",atoms);
186 1 : if( atoms.size()!=4 ) {
187 0 : action->error("number of atoms in box should be equal to four");
188 : }
189 :
190 1 : action->log.printf(" boundaries for region are calculated based on positions of atoms : ");
191 5 : for(unsigned i=0; i<atoms.size(); ++i) {
192 4 : action->log.printf("%d ",atoms[i].serial() );
193 : }
194 1 : action->log.printf("\n");
195 1 : }
196 :
197 1560 : void VolumeCavity::setupRegions( ActionVolume<VolumeCavity>* action, const Pbc& pbc, const std::vector<Vector>& positions ) {
198 : // Make some space for things
199 : Vector d1, d2, d3;
200 :
201 : // Set the position of the origin
202 1560 : Vector origin=positions[0];
203 :
204 : // Get two vectors
205 1560 : d1 = pbc.distance(origin,positions[1]);
206 1560 : double d1l=d1.modulo();
207 1560 : d2 = pbc.distance(origin,positions[2]);
208 :
209 : // Find the vector connecting the origin to the top corner of
210 : // the subregion
211 1560 : d3 = pbc.distance(origin,positions[3]);
212 :
213 : // Create a set of unit vectors
214 1560 : bi = d1 / d1l;
215 1560 : len_bi=dotProduct( d3, bi );
216 1560 : cross = crossProduct( d1, d2 );
217 1560 : double crossmod=cross.modulo();
218 1560 : cross = cross / crossmod;
219 1560 : len_cross=dotProduct( d3, cross );
220 1560 : perp = crossProduct( cross, bi );
221 1560 : len_perp=dotProduct( d3, perp );
222 :
223 : // Calculate derivatives of box shape with respect to atoms
224 1560 : double d1l3=d1l*d1l*d1l;
225 1560 : dbi[0](0,0) = ( -(d1[1]*d1[1]+d1[2]*d1[2])/d1l3 ); // dx/dx
226 1560 : dbi[0](0,1) = ( d1[0]*d1[1]/d1l3 ); // dx/dy
227 1560 : dbi[0](0,2) = ( d1[0]*d1[2]/d1l3 ); // dx/dz
228 1560 : dbi[0](1,0) = ( d1[1]*d1[0]/d1l3 ); // dy/dx
229 1560 : dbi[0](1,1) = ( -(d1[0]*d1[0]+d1[2]*d1[2])/d1l3 ); // dy/dy
230 1560 : dbi[0](1,2) = ( d1[1]*d1[2]/d1l3 );
231 1560 : dbi[0](2,0) = ( d1[2]*d1[0]/d1l3 );
232 1560 : dbi[0](2,1) = ( d1[2]*d1[1]/d1l3 );
233 1560 : dbi[0](2,2) = ( -(d1[1]*d1[1]+d1[0]*d1[0])/d1l3 );
234 :
235 1560 : dbi[1](0,0) = ( (d1[1]*d1[1]+d1[2]*d1[2])/d1l3 );
236 1560 : dbi[1](0,1) = ( -d1[0]*d1[1]/d1l3 );
237 1560 : dbi[1](0,2) = ( -d1[0]*d1[2]/d1l3 );
238 1560 : dbi[1](1,0) = ( -d1[1]*d1[0]/d1l3 );
239 1560 : dbi[1](1,1) = ( (d1[0]*d1[0]+d1[2]*d1[2])/d1l3 );
240 1560 : dbi[1](1,2) = ( -d1[1]*d1[2]/d1l3 );
241 1560 : dbi[1](2,0) = ( -d1[2]*d1[0]/d1l3 );
242 1560 : dbi[1](2,1) = ( -d1[2]*d1[1]/d1l3 );
243 1560 : dbi[1](2,2) = ( (d1[1]*d1[1]+d1[0]*d1[0])/d1l3 );
244 : dbi[2].zero();
245 :
246 : Tensor tcderiv;
247 1560 : double cmod3=crossmod*crossmod*crossmod;
248 : Vector ucross=crossmod*cross;
249 3120 : tcderiv.setCol( 0, crossProduct( d1, Vector(-1.0,0.0,0.0) ) + crossProduct( Vector(-1.0,0.0,0.0), d2 ) );
250 3120 : tcderiv.setCol( 1, crossProduct( d1, Vector(0.0,-1.0,0.0) ) + crossProduct( Vector(0.0,-1.0,0.0), d2 ) );
251 1560 : tcderiv.setCol( 2, crossProduct( d1, Vector(0.0,0.0,-1.0) ) + crossProduct( Vector(0.0,0.0,-1.0), d2 ) );
252 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
253 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
254 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
255 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
256 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
257 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
258 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
259 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
260 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
261 :
262 3120 : tcderiv.setCol( 0, crossProduct( Vector(1.0,0.0,0.0), d2 ) );
263 3120 : tcderiv.setCol( 1, crossProduct( Vector(0.0,1.0,0.0), d2 ) );
264 1560 : tcderiv.setCol( 2, crossProduct( Vector(0.0,0.0,1.0), d2 ) );
265 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
266 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
267 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
268 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
269 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
270 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
271 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
272 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
273 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
274 :
275 3120 : tcderiv.setCol( 0, crossProduct( d1, Vector(1.0,0.0,0.0) ) );
276 3120 : tcderiv.setCol( 1, crossProduct( d1, Vector(0.0,1.0,0.0) ) );
277 1560 : tcderiv.setCol( 2, crossProduct( d1, Vector(0.0,0.0,1.0) ) );
278 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
279 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
280 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
281 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
282 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
283 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
284 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
285 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
286 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
287 :
288 4680 : dperp[0].setCol( 0, ( crossProduct( dcross[0].getCol(0), bi ) + crossProduct( cross, dbi[0].getCol(0) ) ) );
289 4680 : dperp[0].setCol( 1, ( crossProduct( dcross[0].getCol(1), bi ) + crossProduct( cross, dbi[0].getCol(1) ) ) );
290 4680 : dperp[0].setCol( 2, ( crossProduct( dcross[0].getCol(2), bi ) + crossProduct( cross, dbi[0].getCol(2) ) ) );
291 :
292 4680 : dperp[1].setCol( 0, ( crossProduct( dcross[1].getCol(0), bi ) + crossProduct( cross, dbi[1].getCol(0) ) ) );
293 4680 : dperp[1].setCol( 1, ( crossProduct( dcross[1].getCol(1), bi ) + crossProduct( cross, dbi[1].getCol(1) ) ) );
294 4680 : dperp[1].setCol( 2, ( crossProduct( dcross[1].getCol(2), bi ) + crossProduct( cross, dbi[1].getCol(2) ) ) );
295 :
296 3120 : dperp[2].setCol( 0, ( crossProduct( dcross[2].getCol(0), bi ) ) );
297 3120 : dperp[2].setCol( 1, ( crossProduct( dcross[2].getCol(1), bi ) ) );
298 1560 : dperp[2].setCol( 2, ( crossProduct( dcross[2].getCol(2), bi ) ) );
299 :
300 : // Ensure that all lengths are positive
301 1560 : if( len_bi<0 ) {
302 0 : bi=-bi;
303 0 : len_bi=-len_bi;
304 0 : for(unsigned i=0; i<3; ++i) {
305 0 : dbi[i]*=-1.0;
306 : }
307 : }
308 1560 : if( len_cross<0 ) {
309 0 : cross=-cross;
310 0 : len_cross=-len_cross;
311 0 : for(unsigned i=0; i<3; ++i) {
312 0 : dcross[i]*=-1.0;
313 : }
314 : }
315 1560 : if( len_perp<0 ) {
316 0 : perp=-perp;
317 0 : len_perp=-len_perp;
318 0 : for(unsigned i=0; i<3; ++i) {
319 0 : dperp[i]*=-1.0;
320 : }
321 : }
322 1560 : if( len_bi<=0 || len_cross<=0 || len_bi<=0 ) {
323 0 : plumed_merror("Invalid box coordinates");
324 : }
325 :
326 : // Now derivatives of lengths
327 : Tensor dd3( Tensor::identity() );
328 1560 : dlbi[0] = matmul(d3,dbi[0]) - matmul(bi,dd3);
329 1560 : dlbi[1] = matmul(d3,dbi[1]);
330 1560 : dlbi[2] = matmul(d3,dbi[2]);
331 1560 : dlbi[3] = matmul(bi,dd3);
332 :
333 1560 : dlcross[0] = matmul(d3,dcross[0]) - matmul(cross,dd3);
334 1560 : dlcross[1] = matmul(d3,dcross[1]);
335 1560 : dlcross[2] = matmul(d3,dcross[2]);
336 1560 : dlcross[3] = matmul(cross,dd3);
337 :
338 1560 : dlperp[0] = matmul(d3,dperp[0]) - matmul(perp,dd3);
339 1560 : dlperp[1] = matmul(d3,dperp[1]);
340 1560 : dlperp[2] = matmul(d3,dperp[2]);
341 1560 : dlperp[3] = matmul(perp,dd3);
342 :
343 : // Need to calculate the jacobian
344 : Tensor jacob;
345 1560 : jacob(0,0)=bi[0];
346 1560 : jacob(1,0)=bi[1];
347 1560 : jacob(2,0)=bi[2];
348 1560 : jacob(0,1)=cross[0];
349 1560 : jacob(1,1)=cross[1];
350 1560 : jacob(2,1)=cross[2];
351 1560 : jacob(0,2)=perp[0];
352 1560 : jacob(1,2)=perp[1];
353 1560 : jacob(2,2)=perp[2];
354 1560 : jacob_det = fabs( jacob.determinant() );
355 1560 : }
356 :
357 1560 : void VolumeCavity::calculateNumberInside( const VolumeInput& input, const VolumeCavity& actioninput, VolumeOutput& output ) {
358 : // Calculate distance of atom from origin of new coordinate frame
359 1560 : Vector datom=input.pbc.distance( Vector(input.refpos[0][0],input.refpos[0][1],input.refpos[0][2]),
360 1560 : Vector(input.cpos[0],input.cpos[1],input.cpos[2]) );
361 : double ucontr, uder, vcontr, vder, wcontr, wder;
362 :
363 : // Setup the histogram bead
364 1560 : HistogramBead bead( actioninput.kerneltype, 0, actioninput.len_bi, actioninput.sigma);
365 : // Calculate contribution from integral along bi
366 : double upos=dotProduct( datom, actioninput.bi );
367 1560 : ucontr=bead.calculate( upos, uder );
368 1560 : double udlen=bead.uboundDerivative( upos );
369 1560 : double uder2 = bead.lboundDerivative( upos ) - udlen;
370 :
371 : // Calculate contribution from integral along cross
372 1560 : bead.set( 0, actioninput.len_cross, actioninput.sigma );
373 : double vpos=dotProduct( datom, actioninput.cross );
374 1560 : vcontr=bead.calculate( vpos, vder );
375 1560 : double vdlen=bead.uboundDerivative( vpos );
376 1560 : double vder2 = bead.lboundDerivative( vpos ) - vdlen;
377 :
378 : // Calculate contribution from integral along perp
379 1560 : bead.set( 0, actioninput.len_perp, actioninput.sigma );
380 : double wpos=dotProduct( datom, actioninput.perp );
381 1560 : wcontr=bead.calculate( wpos, wder );
382 1560 : double wdlen=bead.uboundDerivative( wpos );
383 1560 : double wder2 = bead.lboundDerivative( wpos ) - wdlen;
384 :
385 : Vector dfd;
386 1560 : dfd[0]=uder*vcontr*wcontr;
387 1560 : dfd[1]=ucontr*vder*wcontr;
388 1560 : dfd[2]=ucontr*vcontr*wder;
389 1560 : output.derivatives[0] = (dfd[0]*actioninput.bi[0]+dfd[1]*actioninput.cross[0]+dfd[2]*actioninput.perp[0]);
390 1560 : output.derivatives[1] = (dfd[0]*actioninput.bi[1]+dfd[1]*actioninput.cross[1]+dfd[2]*actioninput.perp[1]);
391 1560 : output.derivatives[2] = (dfd[0]*actioninput.bi[2]+dfd[1]*actioninput.cross[2]+dfd[2]*actioninput.perp[2]);
392 1560 : output.values[0] = ucontr*vcontr*wcontr*actioninput.jacob_det;
393 :
394 : // Add reference atom derivatives
395 1560 : dfd[0]=uder2*vcontr*wcontr;
396 1560 : dfd[1]=ucontr*vder2*wcontr;
397 1560 : dfd[2]=ucontr*vcontr*wder2;
398 : Vector dfld;
399 1560 : dfld[0]=udlen*vcontr*wcontr;
400 1560 : dfld[1]=ucontr*vdlen*wcontr;
401 1560 : dfld[2]=ucontr*vcontr*wdlen;
402 1560 : output.refders[0] = dfd[0]*matmul(datom,actioninput.dbi[0]) + dfd[1]*matmul(datom,actioninput.dcross[0]) + dfd[2]*matmul(datom,actioninput.dperp[0]) +
403 1560 : dfld[0]*actioninput.dlbi[0] + dfld[1]*actioninput.dlcross[0] + dfld[2]*actioninput.dlperp[0] - Vector(output.derivatives[0],output.derivatives[1],output.derivatives[2]);
404 1560 : output.refders[1] = dfd[0]*matmul(datom,actioninput.dbi[1]) + dfd[1]*matmul(datom,actioninput.dcross[1]) + dfd[2]*matmul(datom,actioninput.dperp[1]) +
405 : dfld[0]*actioninput.dlbi[1] + dfld[1]*actioninput.dlcross[1] + dfld[2]*actioninput.dlperp[1];
406 1560 : output.refders[2] = dfd[0]*matmul(datom,actioninput.dbi[2]) + dfd[1]*matmul(datom,actioninput.dcross[2]) + dfd[2]*matmul(datom,actioninput.dperp[2]) +
407 : dfld[0]*actioninput.dlbi[2] + dfld[1]*actioninput.dlcross[2] + dfld[2]*actioninput.dlperp[2];
408 : output.refders[3] = dfld[0]*actioninput.dlbi[3] + dfld[1]*actioninput.dlcross[3] + dfld[2]*actioninput.dlperp[3];
409 :
410 : Tensor vir;
411 1560 : vir=-Tensor( Vector(input.cpos[0],input.cpos[1],input.cpos[2]), Vector(output.derivatives[0],output.derivatives[1],output.derivatives[2]) );
412 7800 : for(unsigned i=0; i<4; ++i) {
413 12480 : vir -= Tensor( Vector(input.refpos[i][0],input.refpos[i][1],input.refpos[i][2]), Vector(output.refders[i][0],output.refders[i][1],output.refders[i][2]) );
414 : }
415 1560 : output.virial.set( 0, vir );
416 1560 : }
417 :
418 : }
419 : }
|
