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 "VolumeShortcut.h"
28 :
29 : //+PLUMEDOC VOLUMES CAVITY
30 : /*
31 : 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.
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 : In this expression $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. The first of these unit vectors points from the first reference atom to the second.
63 : 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
64 : these first two vectors. $(x_o,y_o,z_o)$, meanwhile, specifies the position of the first reference atom.
65 :
66 : In the first expression above $K$ is one of the kernel functions described in the documentation for the action [BETWEEN](BETWEEN.md)
67 : 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
68 : each of the $u$, $v$ and $w$ directions. This vector connecting atom 1 to atom 4 is projected onto the three unit vectors
69 : described above and the resulting projections determine the $u'$, $v'$ and $w'$ parameters in the above expression.
70 :
71 : 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
72 : in an ion channel in a protein. The extent of the channel is calculated from the positions of atoms 1, 4, 5 and 11.
73 :
74 : ```plumed
75 : cav: CAVITY ATOMS=20-500 BOX=1,4,5,11 SIGMA=0.1
76 : s: SUM ARG=cav PERIODIC=NO
77 : PRINT ARG=s FILE=colvar
78 : ```
79 :
80 : By contrst the following command tells plumed to calculate the coordination numbers (with other water molecules) for the water
81 : molecules in the protein channel described above. The average coordination number and the number of coordination
82 : numbers more than 4 is then calculated for those molecules that are in the region of interest.
83 :
84 : ```plumed
85 : # Calculate the atoms that are in the cavity
86 : cav: CAVITY ATOMS=20-500 BOX=1,4,5,11 SIGMA=0.1
87 : # Calculate the coordination numbers of all the atoms
88 : d1: COORDINATIONNUMBER SPECIES=20-500 SWITCH={RATIONAL R_0=0.1}
89 : # Do atoms have a coordination number that is greater than 4
90 : fd1: MORE_THAN ARG=d1 SWITCH={RATIONAL R_0=4}
91 : # Calculate the mean coordination number in the channel
92 : nnn: CUSTOM ARG=cav,d1 FUNC=x*y PERIODIC=NO
93 : numer: SUM ARG=nnn PERIODIC=NO
94 : denom: SUM ARG=cav PERIODIC=NO
95 : mean: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
96 : # Calculate the number of atoms that are in the channel and that have a coordination number that is greater than 4
97 : sss: CUSTOM ARG=fd1,cav FUNC=x*y PERIODIC=NO
98 : mt: SUM ARG=sss PERIODIC=NO
99 : # And print these two quantities to a file
100 : PRINT ARG=mean,mt FILE=colvar
101 : ```
102 :
103 : As with [AROUND](AROUND.md) earlier version of PLUMED used a different syntax for doing these types of calculations, which can
104 : still be used with this new version of the command. However, we strongly recommend using the newer syntax.
105 :
106 : */
107 : //+ENDPLUMEDOC
108 :
109 : //+PLUMEDOC MCOLVAR CAVITY_CALC
110 : /*
111 : 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
112 :
113 : \par Examples
114 :
115 : */
116 : //+ENDPLUMEDOC
117 :
118 : namespace PLMD {
119 : namespace volumes {
120 :
121 : class VolumeCavity : public ActionVolume {
122 : private:
123 : bool boxout;
124 : OFile boxfile;
125 : double lenunit;
126 : double jacob_det;
127 : double len_bi, len_cross, len_perp, sigma;
128 : Vector origin, bi, cross, perp;
129 : std::vector<Vector> dlbi, dlcross, dlperp;
130 : std::vector<Tensor> dbi, dcross, dperp;
131 : public:
132 : static void registerKeywords( Keywords& keys );
133 : explicit VolumeCavity(const ActionOptions& ao);
134 : ~VolumeCavity();
135 : void setupRegions() override;
136 : void update() override;
137 : double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const override;
138 : };
139 :
140 : PLUMED_REGISTER_ACTION(VolumeCavity,"CAVITY_CALC")
141 : char glob_cavity[] = "CAVITY";
142 : typedef VolumeShortcut<glob_cavity> VolumeCavityShortcut;
143 : PLUMED_REGISTER_ACTION(VolumeCavityShortcut,"CAVITY")
144 :
145 7 : void VolumeCavity::registerKeywords( Keywords& keys ) {
146 7 : ActionVolume::registerKeywords( keys );
147 14 : keys.setDisplayName("CAVITY");
148 7 : keys.add("atoms","BOX","the positions of four atoms that define spatial extent of the cavity");
149 7 : keys.addFlag("PRINT_BOX",false,"write out the positions of the corners of the box to an xyz file");
150 7 : keys.add("optional","FILE","the file on which to write out the box coordinates");
151 7 : keys.add("optional","UNITS","( default=nm ) the units in which to write out the corners of the box");
152 7 : }
153 :
154 1 : VolumeCavity::VolumeCavity(const ActionOptions& ao):
155 : Action(ao),
156 : ActionVolume(ao),
157 1 : boxout(false),
158 1 : lenunit(1.0),
159 1 : dlbi(4),
160 1 : dlcross(4),
161 1 : dlperp(4),
162 1 : dbi(3),
163 1 : dcross(3),
164 2 : dperp(3) {
165 : std::vector<AtomNumber> atoms;
166 2 : parseAtomList("BOX",atoms);
167 1 : if( atoms.size()!=4 ) {
168 0 : error("number of atoms in box should be equal to four");
169 : }
170 :
171 1 : log.printf(" boundaries for region are calculated based on positions of atoms : ");
172 5 : for(unsigned i=0; i<atoms.size(); ++i) {
173 4 : log.printf("%d ",atoms[i].serial() );
174 : }
175 1 : log.printf("\n");
176 1 : requestAtoms( atoms );
177 :
178 1 : boxout=false;
179 1 : parseFlag("PRINT_BOX",boxout);
180 1 : if(boxout) {
181 : std::string boxfname;
182 0 : parse("FILE",boxfname);
183 0 : if(boxfname.length()==0) {
184 0 : error("no name for box file specified");
185 : }
186 : std::string unitname;
187 0 : parse("UNITS",unitname);
188 0 : if ( unitname.length()>0 ) {
189 0 : Units u;
190 0 : u.setLength(unitname);
191 0 : lenunit=getUnits().getLength()/u.getLength();
192 0 : } else {
193 : unitname="nm";
194 : }
195 0 : boxfile.link(*this);
196 0 : boxfile.open( boxfname );
197 0 : log.printf(" printing box coordinates on file named %s in %s \n",boxfname.c_str(), unitname.c_str() );
198 : }
199 :
200 1 : checkRead();
201 1 : }
202 :
203 2 : VolumeCavity::~VolumeCavity() {
204 2 : }
205 :
206 1560 : void VolumeCavity::setupRegions() {
207 : // Make some space for things
208 1560 : Vector d1, d2, d3;
209 :
210 : // Retrieve the sigma value
211 1560 : sigma=getSigma();
212 : // Set the position of the origin
213 1560 : origin=getPosition(0);
214 :
215 : // Get two vectors
216 1560 : d1 = pbcDistance(origin,getPosition(1));
217 1560 : double d1l=d1.modulo();
218 1560 : d2 = pbcDistance(origin,getPosition(2));
219 :
220 : // Find the vector connecting the origin to the top corner of
221 : // the subregion
222 1560 : d3 = pbcDistance(origin,getPosition(3));
223 :
224 : // Create a set of unit vectors
225 1560 : bi = d1 / d1l;
226 1560 : len_bi=dotProduct( d3, bi );
227 1560 : cross = crossProduct( d1, d2 );
228 1560 : double crossmod=cross.modulo();
229 1560 : cross = cross / crossmod;
230 1560 : len_cross=dotProduct( d3, cross );
231 1560 : perp = crossProduct( cross, bi );
232 1560 : len_perp=dotProduct( d3, perp );
233 :
234 : // Calculate derivatives of box shape with respect to atoms
235 1560 : double d1l3=d1l*d1l*d1l;
236 1560 : dbi[0](0,0) = ( -(d1[1]*d1[1]+d1[2]*d1[2])/d1l3 ); // dx/dx
237 1560 : dbi[0](0,1) = ( d1[0]*d1[1]/d1l3 ); // dx/dy
238 1560 : dbi[0](0,2) = ( d1[0]*d1[2]/d1l3 ); // dx/dz
239 1560 : dbi[0](1,0) = ( d1[1]*d1[0]/d1l3 ); // dy/dx
240 1560 : dbi[0](1,1) = ( -(d1[0]*d1[0]+d1[2]*d1[2])/d1l3 ); // dy/dy
241 1560 : dbi[0](1,2) = ( d1[1]*d1[2]/d1l3 );
242 1560 : dbi[0](2,0) = ( d1[2]*d1[0]/d1l3 );
243 1560 : dbi[0](2,1) = ( d1[2]*d1[1]/d1l3 );
244 1560 : dbi[0](2,2) = ( -(d1[1]*d1[1]+d1[0]*d1[0])/d1l3 );
245 :
246 1560 : dbi[1](0,0) = ( (d1[1]*d1[1]+d1[2]*d1[2])/d1l3 );
247 1560 : dbi[1](0,1) = ( -d1[0]*d1[1]/d1l3 );
248 1560 : dbi[1](0,2) = ( -d1[0]*d1[2]/d1l3 );
249 1560 : dbi[1](1,0) = ( -d1[1]*d1[0]/d1l3 );
250 1560 : dbi[1](1,1) = ( (d1[0]*d1[0]+d1[2]*d1[2])/d1l3 );
251 1560 : dbi[1](1,2) = ( -d1[1]*d1[2]/d1l3 );
252 1560 : dbi[1](2,0) = ( -d1[2]*d1[0]/d1l3 );
253 1560 : dbi[1](2,1) = ( -d1[2]*d1[1]/d1l3 );
254 1560 : dbi[1](2,2) = ( (d1[1]*d1[1]+d1[0]*d1[0])/d1l3 );
255 1560 : dbi[2].zero();
256 :
257 1560 : Tensor tcderiv;
258 1560 : double cmod3=crossmod*crossmod*crossmod;
259 1560 : Vector ucross=crossmod*cross;
260 1560 : tcderiv.setCol( 0, crossProduct( d1, Vector(-1.0,0.0,0.0) ) + crossProduct( Vector(-1.0,0.0,0.0), d2 ) );
261 1560 : tcderiv.setCol( 1, crossProduct( d1, Vector(0.0,-1.0,0.0) ) + crossProduct( Vector(0.0,-1.0,0.0), d2 ) );
262 1560 : tcderiv.setCol( 2, crossProduct( d1, Vector(0.0,0.0,-1.0) ) + crossProduct( Vector(0.0,0.0,-1.0), d2 ) );
263 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
264 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
265 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
266 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
267 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
268 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
269 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
270 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
271 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
272 :
273 1560 : tcderiv.setCol( 0, crossProduct( Vector(1.0,0.0,0.0), d2 ) );
274 1560 : tcderiv.setCol( 1, crossProduct( Vector(0.0,1.0,0.0), d2 ) );
275 1560 : tcderiv.setCol( 2, crossProduct( Vector(0.0,0.0,1.0), d2 ) );
276 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
277 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
278 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
279 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
280 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
281 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
282 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
283 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
284 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
285 :
286 1560 : tcderiv.setCol( 0, crossProduct( d1, Vector(1.0,0.0,0.0) ) );
287 1560 : tcderiv.setCol( 1, crossProduct( d1, Vector(0.0,1.0,0.0) ) );
288 1560 : tcderiv.setCol( 2, crossProduct( d1, Vector(0.0,0.0,1.0) ) );
289 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
290 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
291 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
292 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
293 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
294 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
295 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
296 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
297 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
298 :
299 1560 : dperp[0].setCol( 0, ( crossProduct( dcross[0].getCol(0), bi ) + crossProduct( cross, dbi[0].getCol(0) ) ) );
300 1560 : dperp[0].setCol( 1, ( crossProduct( dcross[0].getCol(1), bi ) + crossProduct( cross, dbi[0].getCol(1) ) ) );
301 1560 : dperp[0].setCol( 2, ( crossProduct( dcross[0].getCol(2), bi ) + crossProduct( cross, dbi[0].getCol(2) ) ) );
302 :
303 1560 : dperp[1].setCol( 0, ( crossProduct( dcross[1].getCol(0), bi ) + crossProduct( cross, dbi[1].getCol(0) ) ) );
304 1560 : dperp[1].setCol( 1, ( crossProduct( dcross[1].getCol(1), bi ) + crossProduct( cross, dbi[1].getCol(1) ) ) );
305 1560 : dperp[1].setCol( 2, ( crossProduct( dcross[1].getCol(2), bi ) + crossProduct( cross, dbi[1].getCol(2) ) ) );
306 :
307 1560 : dperp[2].setCol( 0, ( crossProduct( dcross[2].getCol(0), bi ) ) );
308 1560 : dperp[2].setCol( 1, ( crossProduct( dcross[2].getCol(1), bi ) ) );
309 1560 : dperp[2].setCol( 2, ( crossProduct( dcross[2].getCol(2), bi ) ) );
310 :
311 : // Ensure that all lengths are positive
312 1560 : if( len_bi<0 ) {
313 0 : bi=-bi;
314 0 : len_bi=-len_bi;
315 0 : for(unsigned i=0; i<3; ++i) {
316 0 : dbi[i]*=-1.0;
317 : }
318 : }
319 1560 : if( len_cross<0 ) {
320 0 : cross=-cross;
321 0 : len_cross=-len_cross;
322 0 : for(unsigned i=0; i<3; ++i) {
323 0 : dcross[i]*=-1.0;
324 : }
325 : }
326 1560 : if( len_perp<0 ) {
327 0 : perp=-perp;
328 0 : len_perp=-len_perp;
329 0 : for(unsigned i=0; i<3; ++i) {
330 0 : dperp[i]*=-1.0;
331 : }
332 : }
333 1560 : if( len_bi<=0 || len_cross<=0 || len_bi<=0 ) {
334 0 : plumed_merror("Invalid box coordinates");
335 : }
336 :
337 : // Now derivatives of lengths
338 1560 : Tensor dd3( Tensor::identity() );
339 1560 : dlbi[0] = matmul(d3,dbi[0]) - matmul(bi,dd3);
340 1560 : dlbi[1] = matmul(d3,dbi[1]);
341 1560 : dlbi[2] = matmul(d3,dbi[2]);
342 1560 : dlbi[3] = matmul(bi,dd3);
343 :
344 1560 : dlcross[0] = matmul(d3,dcross[0]) - matmul(cross,dd3);
345 1560 : dlcross[1] = matmul(d3,dcross[1]);
346 1560 : dlcross[2] = matmul(d3,dcross[2]);
347 1560 : dlcross[3] = matmul(cross,dd3);
348 :
349 1560 : dlperp[0] = matmul(d3,dperp[0]) - matmul(perp,dd3);
350 1560 : dlperp[1] = matmul(d3,dperp[1]);
351 1560 : dlperp[2] = matmul(d3,dperp[2]);
352 1560 : dlperp[3] = matmul(perp,dd3);
353 :
354 : // Need to calculate the jacobian
355 1560 : Tensor jacob;
356 1560 : jacob(0,0)=bi[0];
357 1560 : jacob(1,0)=bi[1];
358 1560 : jacob(2,0)=bi[2];
359 1560 : jacob(0,1)=cross[0];
360 1560 : jacob(1,1)=cross[1];
361 1560 : jacob(2,1)=cross[2];
362 1560 : jacob(0,2)=perp[0];
363 1560 : jacob(1,2)=perp[1];
364 1560 : jacob(2,2)=perp[2];
365 1560 : jacob_det = fabs( jacob.determinant() );
366 1560 : }
367 :
368 60 : void VolumeCavity::update() {
369 60 : if(boxout) {
370 0 : boxfile.printf("%d\n",8);
371 0 : const Tensor & t(getPbc().getBox());
372 0 : if(getPbc().isOrthorombic()) {
373 0 : boxfile.printf(" %f %f %f\n",lenunit*t(0,0),lenunit*t(1,1),lenunit*t(2,2));
374 : } else {
375 0 : boxfile.printf(" %f %f %f %f %f %f %f %f %f\n",
376 0 : lenunit*t(0,0),lenunit*t(0,1),lenunit*t(0,2),
377 0 : lenunit*t(1,0),lenunit*t(1,1),lenunit*t(1,2),
378 0 : lenunit*t(2,0),lenunit*t(2,1),lenunit*t(2,2)
379 : );
380 : }
381 0 : boxfile.printf("AR %f %f %f \n",lenunit*origin[0],lenunit*origin[1],lenunit*origin[2]);
382 0 : Vector ut, vt, wt;
383 0 : ut = origin + len_bi*bi;
384 0 : vt = origin + len_cross*cross;
385 0 : wt = origin + len_perp*perp;
386 0 : boxfile.printf("AR %f %f %f \n",lenunit*(ut[0]), lenunit*(ut[1]), lenunit*(ut[2]) );
387 0 : boxfile.printf("AR %f %f %f \n",lenunit*(vt[0]), lenunit*(vt[1]), lenunit*(vt[2]) );
388 0 : boxfile.printf("AR %f %f %f \n",lenunit*(wt[0]), lenunit*(wt[1]), lenunit*(wt[2]) );
389 0 : boxfile.printf("AR %f %f %f \n",lenunit*(vt[0]+len_bi*bi[0]),
390 0 : lenunit*(vt[1]+len_bi*bi[1]),
391 0 : lenunit*(vt[2]+len_bi*bi[2]) );
392 0 : boxfile.printf("AR %f %f %f \n",lenunit*(ut[0]+len_perp*perp[0]),
393 0 : lenunit*(ut[1]+len_perp*perp[1]),
394 0 : lenunit*(ut[2]+len_perp*perp[2]) );
395 0 : boxfile.printf("AR %f %f %f \n",lenunit*(vt[0]+len_perp*perp[0]),
396 0 : lenunit*(vt[1]+len_perp*perp[1]),
397 0 : lenunit*(vt[2]+len_perp*perp[2]) );
398 0 : boxfile.printf("AR %f %f %f \n",lenunit*(vt[0]+len_perp*perp[0]+len_bi*bi[0]),
399 0 : lenunit*(vt[1]+len_perp*perp[1]+len_bi*bi[1]),
400 0 : lenunit*(vt[2]+len_perp*perp[2]+len_bi*bi[2]) );
401 : }
402 60 : }
403 :
404 1560 : double VolumeCavity::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& rderiv ) const {
405 : // Setup the histogram bead
406 1560 : HistogramBead bead;
407 : bead.isNotPeriodic();
408 1560 : bead.setKernelType( getKernelType() );
409 :
410 : // Calculate distance of atom from origin of new coordinate frame
411 1560 : Vector datom=pbcDistance( origin, cpos );
412 : double ucontr, uder, vcontr, vder, wcontr, wder;
413 :
414 : // Calculate contribution from integral along bi
415 1560 : bead.set( 0, len_bi, sigma );
416 1560 : double upos=dotProduct( datom, bi );
417 1560 : ucontr=bead.calculate( upos, uder );
418 1560 : double udlen=bead.uboundDerivative( upos );
419 1560 : double uder2 = bead.lboundDerivative( upos ) - udlen;
420 :
421 : // Calculate contribution from integral along cross
422 1560 : bead.set( 0, len_cross, sigma );
423 1560 : double vpos=dotProduct( datom, cross );
424 1560 : vcontr=bead.calculate( vpos, vder );
425 1560 : double vdlen=bead.uboundDerivative( vpos );
426 1560 : double vder2 = bead.lboundDerivative( vpos ) - vdlen;
427 :
428 : // Calculate contribution from integral along perp
429 1560 : bead.set( 0, len_perp, sigma );
430 1560 : double wpos=dotProduct( datom, perp );
431 1560 : wcontr=bead.calculate( wpos, wder );
432 1560 : double wdlen=bead.uboundDerivative( wpos );
433 1560 : double wder2 = bead.lboundDerivative( wpos ) - wdlen;
434 :
435 1560 : Vector dfd;
436 1560 : dfd[0]=uder*vcontr*wcontr;
437 1560 : dfd[1]=ucontr*vder*wcontr;
438 1560 : dfd[2]=ucontr*vcontr*wder;
439 1560 : derivatives[0] = (dfd[0]*bi[0]+dfd[1]*cross[0]+dfd[2]*perp[0]);
440 1560 : derivatives[1] = (dfd[0]*bi[1]+dfd[1]*cross[1]+dfd[2]*perp[1]);
441 1560 : derivatives[2] = (dfd[0]*bi[2]+dfd[1]*cross[2]+dfd[2]*perp[2]);
442 1560 : double tot = ucontr*vcontr*wcontr*jacob_det;
443 :
444 : // Add reference atom derivatives
445 1560 : dfd[0]=uder2*vcontr*wcontr;
446 1560 : dfd[1]=ucontr*vder2*wcontr;
447 1560 : dfd[2]=ucontr*vcontr*wder2;
448 1560 : Vector dfld;
449 1560 : dfld[0]=udlen*vcontr*wcontr;
450 1560 : dfld[1]=ucontr*vdlen*wcontr;
451 1560 : dfld[2]=ucontr*vcontr*wdlen;
452 1560 : rderiv[0] = dfd[0]*matmul(datom,dbi[0]) + dfd[1]*matmul(datom,dcross[0]) + dfd[2]*matmul(datom,dperp[0]) +
453 3120 : dfld[0]*dlbi[0] + dfld[1]*dlcross[0] + dfld[2]*dlperp[0] - derivatives;
454 1560 : rderiv[1] = dfd[0]*matmul(datom,dbi[1]) + dfd[1]*matmul(datom,dcross[1]) + dfd[2]*matmul(datom,dperp[1]) +
455 3120 : dfld[0]*dlbi[1] + dfld[1]*dlcross[1] + dfld[2]*dlperp[1];
456 1560 : rderiv[2] = dfd[0]*matmul(datom,dbi[2]) + dfd[1]*matmul(datom,dcross[2]) + dfd[2]*matmul(datom,dperp[2]) +
457 3120 : dfld[0]*dlbi[2] + dfld[1]*dlcross[2] + dfld[2]*dlperp[2];
458 1560 : rderiv[3] = dfld[0]*dlbi[3] + dfld[1]*dlcross[3] + dfld[2]*dlperp[3];
459 :
460 1560 : vir.zero();
461 1560 : vir-=Tensor( cpos,derivatives );
462 7800 : for(unsigned i=0; i<4; ++i) {
463 6240 : vir -= Tensor( getPosition(i), rderiv[i] );
464 : }
465 :
466 1560 : return tot;
467 : }
468 :
469 : }
470 : }
|
