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 : Each of the base quantities calculated by a multicolvar can can be assigned to a particular point in three
34 : dimensional space. For example, if we have the coordination numbers for all the atoms in the
35 : system each coordination number can be assumed to lie on the position of the central atom.
36 : Because each base quantity can be assigned to a particular point in space we can calculate functions of the
37 : distribution of base quantities in a particular part of the box by using:
38 :
39 : \f[
40 : \overline{s}_{\tau} = \frac{ \sum_i f(s_i) w(u_i,v_i,w_i) }{ \sum_i w(u_i,v_i,w_i) }
41 : \f]
42 :
43 : where the sum is over the collective variables, \f$s_i\f$, each of which can be thought to be at \f$ (u_i,v_i,z_i)\f$.
44 : The function \f$(s_i)\f$ can be any of the usual LESS_THAN, MORE_THAN, WITHIN etc that are used in all other multicolvars.
45 : Notice that here (at variance with what is done in \ref AROUND) we have transformed from the usual \f$(x_i,y_i,z_i)\f$
46 : position to a position in \f$ (u_i,v_i,z_i)\f$. This is done using a rotation matrix as follows:
47 :
48 : \f[
49 : \left(
50 : \begin{matrix}
51 : u_i \\
52 : v_i \\
53 : w_i
54 : \end{matrix}
55 : \right) = \mathbf{R}
56 : \left(
57 : \begin{matrix}
58 : x_i - x_o \\
59 : y_i - y_o \\
60 : z_i - z_o
61 : \end{matrix}
62 : \right)
63 : \f]
64 :
65 : where \f$\mathbf{R}\f$ is a rotation matrix that is calculated by constructing a set of three orthonormal vectors from the
66 : reference positions specified by the user. The first of these unit vectors points from the first reference atom to the second.
67 : 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
68 : these first two vectors. \f$(x_o,y_o,z_o)\f$, meanwhile, specifies the position of the first reference atom.
69 :
70 : In the previous function \f$ w(u_i,v_i,w_i) \f$ measures whether or not the system is in the subregion of interest. It
71 : is equal to:
72 :
73 : \f[
74 : 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
75 : K\left( \frac{u - u_i}{\sigma} \right)K\left( \frac{v - v_i}{\sigma} \right)K\left( \frac{w - w_i}{\sigma} \right)
76 : \f]
77 :
78 : where \f$K\f$ is one of the kernel functions described on \ref histogrambead and \f$\sigma\f$ is a bandwidth parameter.
79 : The vector connecting atom 1 to atom 4 is used to define the extent of the box in each of the \f$u\f$, \f$v\f$ and \f$w\f$
80 : directions. Essentially the vector connecting atom 1 to atom 4 is projected onto the three unit vectors
81 : described above and the resulting projections determine the \f$u'\f$, \f$v'\f$ and \f$w'\f$ parameters in the above expression.
82 :
83 : \par Examples
84 :
85 : The following commands tell plumed to calculate the number of atoms in an ion channel in a protein.
86 : The extent of the channel is calculated from the positions of atoms 1, 4, 5 and 11. The final value will be labeled cav.
87 :
88 : \plumedfile
89 : d1: DENSITY SPECIES=20-500
90 : CAVITY DATA=d1 ATOMS=1,4,5,11 SIGMA=0.1 LABEL=cav
91 : \endplumedfile
92 :
93 : The following command tells plumed to calculate the coordination numbers (with other water molecules) for the water
94 : molecules in the protein channel described above. The average coordination number and the number of coordination
95 : numbers more than 4 is then calculated. The values of these two quantities are given the labels cav.mean and cav.morethan
96 :
97 : \plumedfile
98 : d1: COORDINATIONNUMBER SPECIES=20-500 R_0=0.1
99 : CAVITY DATA=d1 ATOMS=1,4,5,11 SIGMA=0.1 MEAN MORE_THAN={RATIONAL R_0=4} LABEL=cav
100 : \endplumedfile
101 :
102 : */
103 : //+ENDPLUMEDOC
104 :
105 : //+PLUMEDOC MCOLVAR CAVITY_CALC
106 : /*
107 : 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
108 :
109 : \par Examples
110 :
111 : */
112 : //+ENDPLUMEDOC
113 :
114 : namespace PLMD {
115 : namespace volumes {
116 :
117 : class VolumeCavity : public ActionVolume {
118 : private:
119 : bool boxout;
120 : OFile boxfile;
121 : double lenunit;
122 : double jacob_det;
123 : double len_bi, len_cross, len_perp, sigma;
124 : Vector origin, bi, cross, perp;
125 : std::vector<Vector> dlbi, dlcross, dlperp;
126 : std::vector<Tensor> dbi, dcross, dperp;
127 : public:
128 : static void registerKeywords( Keywords& keys );
129 : explicit VolumeCavity(const ActionOptions& ao);
130 : ~VolumeCavity();
131 : void setupRegions() override;
132 : void update() override;
133 : double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const override;
134 : };
135 :
136 : PLUMED_REGISTER_ACTION(VolumeCavity,"CAVITY_CALC")
137 : char glob_cavity[] = "CAVITY";
138 : typedef VolumeShortcut<glob_cavity> VolumeCavityShortcut;
139 : PLUMED_REGISTER_ACTION(VolumeCavityShortcut,"CAVITY")
140 :
141 7 : void VolumeCavity::registerKeywords( Keywords& keys ) {
142 7 : ActionVolume::registerKeywords( keys );
143 7 : keys.setDisplayName("CAVITY");
144 14 : keys.add("atoms","BOX","the positions of four atoms that define spatial extent of the cavity");
145 14 : keys.addFlag("PRINT_BOX",false,"write out the positions of the corners of the box to an xyz file");
146 14 : keys.add("optional","FILE","the file on which to write out the box coordinates");
147 14 : keys.add("optional","UNITS","( default=nm ) the units in which to write out the corners of the box");
148 7 : }
149 :
150 1 : VolumeCavity::VolumeCavity(const ActionOptions& ao):
151 : Action(ao),
152 : ActionVolume(ao),
153 1 : boxout(false),
154 1 : lenunit(1.0),
155 1 : dlbi(4),
156 1 : dlcross(4),
157 1 : dlperp(4),
158 1 : dbi(3),
159 1 : dcross(3),
160 2 : dperp(3) {
161 : std::vector<AtomNumber> atoms;
162 2 : parseAtomList("BOX",atoms);
163 1 : if( atoms.size()!=4 ) {
164 0 : error("number of atoms in box should be equal to four");
165 : }
166 :
167 1 : log.printf(" boundaries for region are calculated based on positions of atoms : ");
168 5 : for(unsigned i=0; i<atoms.size(); ++i) {
169 4 : log.printf("%d ",atoms[i].serial() );
170 : }
171 1 : log.printf("\n");
172 1 : requestAtoms( atoms );
173 :
174 1 : boxout=false;
175 1 : parseFlag("PRINT_BOX",boxout);
176 1 : if(boxout) {
177 : std::string boxfname;
178 0 : parse("FILE",boxfname);
179 0 : if(boxfname.length()==0) {
180 0 : error("no name for box file specified");
181 : }
182 : std::string unitname;
183 0 : parse("UNITS",unitname);
184 0 : if ( unitname.length()>0 ) {
185 0 : Units u;
186 0 : u.setLength(unitname);
187 0 : lenunit=getUnits().getLength()/u.getLength();
188 0 : } else {
189 : unitname="nm";
190 : }
191 0 : boxfile.link(*this);
192 0 : boxfile.open( boxfname );
193 0 : log.printf(" printing box coordinates on file named %s in %s \n",boxfname.c_str(), unitname.c_str() );
194 : }
195 :
196 1 : checkRead();
197 1 : }
198 :
199 2 : VolumeCavity::~VolumeCavity() {
200 2 : }
201 :
202 1560 : void VolumeCavity::setupRegions() {
203 : // Make some space for things
204 1560 : Vector d1, d2, d3;
205 :
206 : // Retrieve the sigma value
207 1560 : sigma=getSigma();
208 : // Set the position of the origin
209 1560 : origin=getPosition(0);
210 :
211 : // Get two vectors
212 1560 : d1 = pbcDistance(origin,getPosition(1));
213 1560 : double d1l=d1.modulo();
214 1560 : d2 = pbcDistance(origin,getPosition(2));
215 :
216 : // Find the vector connecting the origin to the top corner of
217 : // the subregion
218 1560 : d3 = pbcDistance(origin,getPosition(3));
219 :
220 : // Create a set of unit vectors
221 1560 : bi = d1 / d1l;
222 1560 : len_bi=dotProduct( d3, bi );
223 1560 : cross = crossProduct( d1, d2 );
224 1560 : double crossmod=cross.modulo();
225 1560 : cross = cross / crossmod;
226 1560 : len_cross=dotProduct( d3, cross );
227 1560 : perp = crossProduct( cross, bi );
228 1560 : len_perp=dotProduct( d3, perp );
229 :
230 : // Calculate derivatives of box shape with respect to atoms
231 1560 : double d1l3=d1l*d1l*d1l;
232 1560 : dbi[0](0,0) = ( -(d1[1]*d1[1]+d1[2]*d1[2])/d1l3 ); // dx/dx
233 1560 : dbi[0](0,1) = ( d1[0]*d1[1]/d1l3 ); // dx/dy
234 1560 : dbi[0](0,2) = ( d1[0]*d1[2]/d1l3 ); // dx/dz
235 1560 : dbi[0](1,0) = ( d1[1]*d1[0]/d1l3 ); // dy/dx
236 1560 : dbi[0](1,1) = ( -(d1[0]*d1[0]+d1[2]*d1[2])/d1l3 ); // dy/dy
237 1560 : dbi[0](1,2) = ( d1[1]*d1[2]/d1l3 );
238 1560 : dbi[0](2,0) = ( d1[2]*d1[0]/d1l3 );
239 1560 : dbi[0](2,1) = ( d1[2]*d1[1]/d1l3 );
240 1560 : dbi[0](2,2) = ( -(d1[1]*d1[1]+d1[0]*d1[0])/d1l3 );
241 :
242 1560 : dbi[1](0,0) = ( (d1[1]*d1[1]+d1[2]*d1[2])/d1l3 );
243 1560 : dbi[1](0,1) = ( -d1[0]*d1[1]/d1l3 );
244 1560 : dbi[1](0,2) = ( -d1[0]*d1[2]/d1l3 );
245 1560 : dbi[1](1,0) = ( -d1[1]*d1[0]/d1l3 );
246 1560 : dbi[1](1,1) = ( (d1[0]*d1[0]+d1[2]*d1[2])/d1l3 );
247 1560 : dbi[1](1,2) = ( -d1[1]*d1[2]/d1l3 );
248 1560 : dbi[1](2,0) = ( -d1[2]*d1[0]/d1l3 );
249 1560 : dbi[1](2,1) = ( -d1[2]*d1[1]/d1l3 );
250 1560 : dbi[1](2,2) = ( (d1[1]*d1[1]+d1[0]*d1[0])/d1l3 );
251 1560 : dbi[2].zero();
252 :
253 1560 : Tensor tcderiv;
254 1560 : double cmod3=crossmod*crossmod*crossmod;
255 1560 : Vector ucross=crossmod*cross;
256 1560 : tcderiv.setCol( 0, crossProduct( d1, Vector(-1.0,0.0,0.0) ) + crossProduct( Vector(-1.0,0.0,0.0), d2 ) );
257 1560 : 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 1560 : tcderiv.setCol( 0, crossProduct( Vector(1.0,0.0,0.0), d2 ) );
270 1560 : 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 1560 : tcderiv.setCol( 0, crossProduct( d1, Vector(1.0,0.0,0.0) ) );
283 1560 : 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 : dperp[0].setCol( 0, ( crossProduct( dcross[0].getCol(0), bi ) + crossProduct( cross, dbi[0].getCol(0) ) ) );
296 1560 : dperp[0].setCol( 1, ( crossProduct( dcross[0].getCol(1), bi ) + crossProduct( cross, dbi[0].getCol(1) ) ) );
297 1560 : dperp[0].setCol( 2, ( crossProduct( dcross[0].getCol(2), bi ) + crossProduct( cross, dbi[0].getCol(2) ) ) );
298 :
299 1560 : dperp[1].setCol( 0, ( crossProduct( dcross[1].getCol(0), bi ) + crossProduct( cross, dbi[1].getCol(0) ) ) );
300 1560 : dperp[1].setCol( 1, ( crossProduct( dcross[1].getCol(1), bi ) + crossProduct( cross, dbi[1].getCol(1) ) ) );
301 1560 : dperp[1].setCol( 2, ( crossProduct( dcross[1].getCol(2), bi ) + crossProduct( cross, dbi[1].getCol(2) ) ) );
302 :
303 1560 : dperp[2].setCol( 0, ( crossProduct( dcross[2].getCol(0), bi ) ) );
304 1560 : dperp[2].setCol( 1, ( crossProduct( dcross[2].getCol(1), bi ) ) );
305 1560 : dperp[2].setCol( 2, ( crossProduct( dcross[2].getCol(2), bi ) ) );
306 :
307 : // Ensure that all lengths are positive
308 1560 : if( len_bi<0 ) {
309 0 : bi=-bi;
310 0 : len_bi=-len_bi;
311 0 : for(unsigned i=0; i<3; ++i) {
312 0 : dbi[i]*=-1.0;
313 : }
314 : }
315 1560 : if( len_cross<0 ) {
316 0 : cross=-cross;
317 0 : len_cross=-len_cross;
318 0 : for(unsigned i=0; i<3; ++i) {
319 0 : dcross[i]*=-1.0;
320 : }
321 : }
322 1560 : if( len_perp<0 ) {
323 0 : perp=-perp;
324 0 : len_perp=-len_perp;
325 0 : for(unsigned i=0; i<3; ++i) {
326 0 : dperp[i]*=-1.0;
327 : }
328 : }
329 1560 : if( len_bi<=0 || len_cross<=0 || len_bi<=0 ) {
330 0 : plumed_merror("Invalid box coordinates");
331 : }
332 :
333 : // Now derivatives of lengths
334 1560 : Tensor dd3( Tensor::identity() );
335 1560 : dlbi[0] = matmul(d3,dbi[0]) - matmul(bi,dd3);
336 1560 : dlbi[1] = matmul(d3,dbi[1]);
337 1560 : dlbi[2] = matmul(d3,dbi[2]);
338 1560 : dlbi[3] = matmul(bi,dd3);
339 :
340 1560 : dlcross[0] = matmul(d3,dcross[0]) - matmul(cross,dd3);
341 1560 : dlcross[1] = matmul(d3,dcross[1]);
342 1560 : dlcross[2] = matmul(d3,dcross[2]);
343 1560 : dlcross[3] = matmul(cross,dd3);
344 :
345 1560 : dlperp[0] = matmul(d3,dperp[0]) - matmul(perp,dd3);
346 1560 : dlperp[1] = matmul(d3,dperp[1]);
347 1560 : dlperp[2] = matmul(d3,dperp[2]);
348 1560 : dlperp[3] = matmul(perp,dd3);
349 :
350 : // Need to calculate the jacobian
351 1560 : Tensor jacob;
352 1560 : jacob(0,0)=bi[0];
353 1560 : jacob(1,0)=bi[1];
354 1560 : jacob(2,0)=bi[2];
355 1560 : jacob(0,1)=cross[0];
356 1560 : jacob(1,1)=cross[1];
357 1560 : jacob(2,1)=cross[2];
358 1560 : jacob(0,2)=perp[0];
359 1560 : jacob(1,2)=perp[1];
360 1560 : jacob(2,2)=perp[2];
361 1560 : jacob_det = fabs( jacob.determinant() );
362 1560 : }
363 :
364 60 : void VolumeCavity::update() {
365 60 : if(boxout) {
366 0 : boxfile.printf("%d\n",8);
367 0 : const Tensor & t(getPbc().getBox());
368 0 : if(getPbc().isOrthorombic()) {
369 0 : boxfile.printf(" %f %f %f\n",lenunit*t(0,0),lenunit*t(1,1),lenunit*t(2,2));
370 : } else {
371 0 : boxfile.printf(" %f %f %f %f %f %f %f %f %f\n",
372 0 : lenunit*t(0,0),lenunit*t(0,1),lenunit*t(0,2),
373 0 : lenunit*t(1,0),lenunit*t(1,1),lenunit*t(1,2),
374 0 : lenunit*t(2,0),lenunit*t(2,1),lenunit*t(2,2)
375 : );
376 : }
377 0 : boxfile.printf("AR %f %f %f \n",lenunit*origin[0],lenunit*origin[1],lenunit*origin[2]);
378 0 : Vector ut, vt, wt;
379 0 : ut = origin + len_bi*bi;
380 0 : vt = origin + len_cross*cross;
381 0 : wt = origin + len_perp*perp;
382 0 : boxfile.printf("AR %f %f %f \n",lenunit*(ut[0]), lenunit*(ut[1]), lenunit*(ut[2]) );
383 0 : boxfile.printf("AR %f %f %f \n",lenunit*(vt[0]), lenunit*(vt[1]), lenunit*(vt[2]) );
384 0 : boxfile.printf("AR %f %f %f \n",lenunit*(wt[0]), lenunit*(wt[1]), lenunit*(wt[2]) );
385 0 : boxfile.printf("AR %f %f %f \n",lenunit*(vt[0]+len_bi*bi[0]),
386 0 : lenunit*(vt[1]+len_bi*bi[1]),
387 0 : lenunit*(vt[2]+len_bi*bi[2]) );
388 0 : boxfile.printf("AR %f %f %f \n",lenunit*(ut[0]+len_perp*perp[0]),
389 0 : lenunit*(ut[1]+len_perp*perp[1]),
390 0 : lenunit*(ut[2]+len_perp*perp[2]) );
391 0 : boxfile.printf("AR %f %f %f \n",lenunit*(vt[0]+len_perp*perp[0]),
392 0 : lenunit*(vt[1]+len_perp*perp[1]),
393 0 : lenunit*(vt[2]+len_perp*perp[2]) );
394 0 : boxfile.printf("AR %f %f %f \n",lenunit*(vt[0]+len_perp*perp[0]+len_bi*bi[0]),
395 0 : lenunit*(vt[1]+len_perp*perp[1]+len_bi*bi[1]),
396 0 : lenunit*(vt[2]+len_perp*perp[2]+len_bi*bi[2]) );
397 : }
398 60 : }
399 :
400 1560 : double VolumeCavity::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& rderiv ) const {
401 : // Setup the histogram bead
402 1560 : HistogramBead bead;
403 : bead.isNotPeriodic();
404 1560 : bead.setKernelType( getKernelType() );
405 :
406 : // Calculate distance of atom from origin of new coordinate frame
407 1560 : Vector datom=pbcDistance( origin, cpos );
408 : double ucontr, uder, vcontr, vder, wcontr, wder;
409 :
410 : // Calculate contribution from integral along bi
411 1560 : bead.set( 0, len_bi, sigma );
412 1560 : double upos=dotProduct( datom, bi );
413 1560 : ucontr=bead.calculate( upos, uder );
414 1560 : double udlen=bead.uboundDerivative( upos );
415 1560 : double uder2 = bead.lboundDerivative( upos ) - udlen;
416 :
417 : // Calculate contribution from integral along cross
418 1560 : bead.set( 0, len_cross, sigma );
419 1560 : double vpos=dotProduct( datom, cross );
420 1560 : vcontr=bead.calculate( vpos, vder );
421 1560 : double vdlen=bead.uboundDerivative( vpos );
422 1560 : double vder2 = bead.lboundDerivative( vpos ) - vdlen;
423 :
424 : // Calculate contribution from integral along perp
425 1560 : bead.set( 0, len_perp, sigma );
426 1560 : double wpos=dotProduct( datom, perp );
427 1560 : wcontr=bead.calculate( wpos, wder );
428 1560 : double wdlen=bead.uboundDerivative( wpos );
429 1560 : double wder2 = bead.lboundDerivative( wpos ) - wdlen;
430 :
431 1560 : Vector dfd;
432 1560 : dfd[0]=uder*vcontr*wcontr;
433 1560 : dfd[1]=ucontr*vder*wcontr;
434 1560 : dfd[2]=ucontr*vcontr*wder;
435 1560 : derivatives[0] = (dfd[0]*bi[0]+dfd[1]*cross[0]+dfd[2]*perp[0]);
436 1560 : derivatives[1] = (dfd[0]*bi[1]+dfd[1]*cross[1]+dfd[2]*perp[1]);
437 1560 : derivatives[2] = (dfd[0]*bi[2]+dfd[1]*cross[2]+dfd[2]*perp[2]);
438 1560 : double tot = ucontr*vcontr*wcontr*jacob_det;
439 :
440 : // Add reference atom derivatives
441 1560 : dfd[0]=uder2*vcontr*wcontr;
442 1560 : dfd[1]=ucontr*vder2*wcontr;
443 1560 : dfd[2]=ucontr*vcontr*wder2;
444 1560 : Vector dfld;
445 1560 : dfld[0]=udlen*vcontr*wcontr;
446 1560 : dfld[1]=ucontr*vdlen*wcontr;
447 1560 : dfld[2]=ucontr*vcontr*wdlen;
448 1560 : rderiv[0] = dfd[0]*matmul(datom,dbi[0]) + dfd[1]*matmul(datom,dcross[0]) + dfd[2]*matmul(datom,dperp[0]) +
449 3120 : dfld[0]*dlbi[0] + dfld[1]*dlcross[0] + dfld[2]*dlperp[0] - derivatives;
450 1560 : rderiv[1] = dfd[0]*matmul(datom,dbi[1]) + dfd[1]*matmul(datom,dcross[1]) + dfd[2]*matmul(datom,dperp[1]) +
451 3120 : dfld[0]*dlbi[1] + dfld[1]*dlcross[1] + dfld[2]*dlperp[1];
452 1560 : rderiv[2] = dfd[0]*matmul(datom,dbi[2]) + dfd[1]*matmul(datom,dcross[2]) + dfd[2]*matmul(datom,dperp[2]) +
453 3120 : dfld[0]*dlbi[2] + dfld[1]*dlcross[2] + dfld[2]*dlperp[2];
454 1560 : rderiv[3] = dfld[0]*dlbi[3] + dfld[1]*dlcross[3] + dfld[2]*dlperp[3];
455 :
456 1560 : vir.zero();
457 1560 : vir-=Tensor( cpos,derivatives );
458 7800 : for(unsigned i=0; i<4; ++i) {
459 6240 : vir -= Tensor( getPosition(i), rderiv[i] );
460 : }
461 :
462 1560 : return tot;
463 : }
464 :
465 : }
466 : }
|
