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 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 : using precision=double;
135 : typedef VolumeIn<precision> VolumeInput;
136 : typedef VolumeOut<precision> VolumeOutput;
137 : double jacob_det;
138 : double len_bi, len_cross, len_perp, sigma;
139 : Vector bi, cross, perp;
140 : HistogramBead::KernelType kerneltype{HistogramBead::KernelType::gaussian};
141 : std::vector<Vector> dlbi, dlcross, dlperp;
142 : std::vector<Tensor> dbi, dcross, dperp;
143 : static void registerKeywords( Keywords& keys );
144 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) {}
145 : void setupRegions( ActionVolume<VolumeCavity>* action, const Pbc& pbc, const std::vector<Vector>& positions );
146 : void parseInput( ActionVolume<VolumeCavity>* action );
147 : static void parseAtoms( ActionVolume<VolumeCavity>* action, std::vector<AtomNumber>& atoms );
148 : // the assignment is customized (see below): copy construction is the default one
149 : VolumeCavity(const VolumeCavity&) = default;
150 1 : VolumeCavity& operator=( const VolumeCavity& m ) {
151 1 : jacob_det=m.jacob_det;
152 1 : len_bi=m.len_bi;
153 1 : len_cross=m.len_cross;
154 1 : len_perp=m.len_perp;
155 1 : sigma=m.sigma;
156 1 : dlbi.resize(4);
157 1 : dlcross.resize(4);
158 1 : dlperp.resize(4);
159 1 : dbi.resize(3);
160 1 : dcross.resize(3);
161 1 : dperp.resize(3);
162 1 : kerneltype=m.kerneltype;
163 1 : return *this;
164 : }
165 : static void calculateNumberInside( const VolumeInput& input, const VolumeCavity& actioninput, VolumeOutput& output );
166 : };
167 :
168 : typedef ActionVolume<VolumeCavity> VolCav;
169 : PLUMED_REGISTER_ACTION(VolCav,"CAVITY_CALC")
170 : char glob_cavity[] = "CAVITY";
171 : typedef VolumeShortcut<glob_cavity> VolumeCavityShortcut;
172 : PLUMED_REGISTER_ACTION(VolumeCavityShortcut,"CAVITY")
173 :
174 11 : void VolumeCavity::registerKeywords( Keywords& keys ) {
175 22 : keys.setDisplayName("CAVITY");
176 11 : keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
177 11 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
178 11 : keys.add("atoms","BOX","the positions of four atoms that define spatial extent of the cavity");
179 11 : }
180 :
181 1 : void VolumeCavity::parseInput( ActionVolume<VolumeCavity>* action ) {
182 2 : action->parse("SIGMA",sigma);
183 : std::string mykerneltype;
184 1 : action->parse("KERNEL",mykerneltype);
185 1 : kerneltype=HistogramBead::getKernelType(mykerneltype);
186 1 : action->log.printf(" using %s kernels with a bandwidth of %f \n", mykerneltype.c_str(), sigma );
187 1 : }
188 :
189 1 : void VolumeCavity::parseAtoms( ActionVolume<VolumeCavity>* action, std::vector<AtomNumber>& atoms ) {
190 2 : action->parseAtomList("BOX",atoms);
191 1 : if( atoms.size()!=4 ) {
192 0 : action->error("number of atoms in box should be equal to four");
193 : }
194 :
195 1 : action->log.printf(" boundaries for region are calculated based on positions of atoms : ");
196 5 : for(unsigned i=0; i<atoms.size(); ++i) {
197 4 : action->log.printf("%d ",atoms[i].serial() );
198 : }
199 1 : action->log.printf("\n");
200 1 : }
201 :
202 1560 : void VolumeCavity::setupRegions( ActionVolume<VolumeCavity>* action, const Pbc& pbc, const std::vector<Vector>& positions ) {
203 : // Make some space for things
204 1560 : Vector d1, d2, d3;
205 :
206 : // Set the position of the origin
207 1560 : Vector origin=positions[0];
208 :
209 : // Get two vectors
210 1560 : d1 = pbc.distance(origin,positions[1]);
211 1560 : double d1l=d1.modulo();
212 1560 : d2 = pbc.distance(origin,positions[2]);
213 :
214 : // Find the vector connecting the origin to the top corner of
215 : // the subregion
216 1560 : d3 = pbc.distance(origin,positions[3]);
217 :
218 : // Create a set of unit vectors
219 1560 : bi = d1 / d1l;
220 1560 : len_bi=dotProduct( d3, bi );
221 1560 : cross = crossProduct( d1, d2 );
222 1560 : double crossmod=cross.modulo();
223 1560 : cross = cross / crossmod;
224 1560 : len_cross=dotProduct( d3, cross );
225 1560 : perp = crossProduct( cross, bi );
226 1560 : len_perp=dotProduct( d3, perp );
227 :
228 : // Calculate derivatives of box shape with respect to atoms
229 1560 : double d1l3=d1l*d1l*d1l;
230 1560 : dbi[0](0,0) = ( -(d1[1]*d1[1]+d1[2]*d1[2])/d1l3 ); // dx/dx
231 1560 : dbi[0](0,1) = ( d1[0]*d1[1]/d1l3 ); // dx/dy
232 1560 : dbi[0](0,2) = ( d1[0]*d1[2]/d1l3 ); // dx/dz
233 1560 : dbi[0](1,0) = ( d1[1]*d1[0]/d1l3 ); // dy/dx
234 1560 : dbi[0](1,1) = ( -(d1[0]*d1[0]+d1[2]*d1[2])/d1l3 ); // dy/dy
235 1560 : dbi[0](1,2) = ( d1[1]*d1[2]/d1l3 );
236 1560 : dbi[0](2,0) = ( d1[2]*d1[0]/d1l3 );
237 1560 : dbi[0](2,1) = ( d1[2]*d1[1]/d1l3 );
238 1560 : dbi[0](2,2) = ( -(d1[1]*d1[1]+d1[0]*d1[0])/d1l3 );
239 :
240 1560 : dbi[1](0,0) = ( (d1[1]*d1[1]+d1[2]*d1[2])/d1l3 );
241 1560 : dbi[1](0,1) = ( -d1[0]*d1[1]/d1l3 );
242 1560 : dbi[1](0,2) = ( -d1[0]*d1[2]/d1l3 );
243 1560 : dbi[1](1,0) = ( -d1[1]*d1[0]/d1l3 );
244 1560 : dbi[1](1,1) = ( (d1[0]*d1[0]+d1[2]*d1[2])/d1l3 );
245 1560 : dbi[1](1,2) = ( -d1[1]*d1[2]/d1l3 );
246 1560 : dbi[1](2,0) = ( -d1[2]*d1[0]/d1l3 );
247 1560 : dbi[1](2,1) = ( -d1[2]*d1[1]/d1l3 );
248 1560 : dbi[1](2,2) = ( (d1[1]*d1[1]+d1[0]*d1[0])/d1l3 );
249 : dbi[2].zero();
250 :
251 1560 : Tensor tcderiv;
252 1560 : double cmod3=crossmod*crossmod*crossmod;
253 : Vector ucross=crossmod*cross;
254 3120 : tcderiv.setCol( 0, crossProduct( d1, Vector(-1.0,0.0,0.0) ) + crossProduct( Vector(-1.0,0.0,0.0), d2 ) );
255 3120 : tcderiv.setCol( 1, crossProduct( d1, Vector(0.0,-1.0,0.0) ) + crossProduct( Vector(0.0,-1.0,0.0), d2 ) );
256 1560 : tcderiv.setCol( 2, crossProduct( d1, Vector(0.0,0.0,-1.0) ) + crossProduct( Vector(0.0,0.0,-1.0), d2 ) );
257 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
258 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
259 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
260 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
261 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
262 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
263 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
264 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
265 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
266 :
267 3120 : tcderiv.setCol( 0, crossProduct( Vector(1.0,0.0,0.0), d2 ) );
268 3120 : tcderiv.setCol( 1, crossProduct( Vector(0.0,1.0,0.0), d2 ) );
269 1560 : tcderiv.setCol( 2, crossProduct( Vector(0.0,0.0,1.0), d2 ) );
270 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
271 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
272 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
273 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
274 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
275 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
276 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
277 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
278 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
279 :
280 3120 : tcderiv.setCol( 0, crossProduct( d1, Vector(1.0,0.0,0.0) ) );
281 3120 : tcderiv.setCol( 1, crossProduct( d1, Vector(0.0,1.0,0.0) ) );
282 1560 : tcderiv.setCol( 2, crossProduct( d1, Vector(0.0,0.0,1.0) ) );
283 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
284 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
285 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
286 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
287 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
288 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
289 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
290 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
291 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
292 :
293 4680 : dperp[0].setCol( 0, ( crossProduct( dcross[0].getCol(0), bi ) + crossProduct( cross, dbi[0].getCol(0) ) ) );
294 4680 : dperp[0].setCol( 1, ( crossProduct( dcross[0].getCol(1), bi ) + crossProduct( cross, dbi[0].getCol(1) ) ) );
295 4680 : dperp[0].setCol( 2, ( crossProduct( dcross[0].getCol(2), bi ) + crossProduct( cross, dbi[0].getCol(2) ) ) );
296 :
297 4680 : dperp[1].setCol( 0, ( crossProduct( dcross[1].getCol(0), bi ) + crossProduct( cross, dbi[1].getCol(0) ) ) );
298 4680 : dperp[1].setCol( 1, ( crossProduct( dcross[1].getCol(1), bi ) + crossProduct( cross, dbi[1].getCol(1) ) ) );
299 4680 : dperp[1].setCol( 2, ( crossProduct( dcross[1].getCol(2), bi ) + crossProduct( cross, dbi[1].getCol(2) ) ) );
300 :
301 3120 : dperp[2].setCol( 0, ( crossProduct( dcross[2].getCol(0), bi ) ) );
302 3120 : dperp[2].setCol( 1, ( crossProduct( dcross[2].getCol(1), bi ) ) );
303 1560 : dperp[2].setCol( 2, ( crossProduct( dcross[2].getCol(2), bi ) ) );
304 :
305 : // Ensure that all lengths are positive
306 1560 : if( len_bi<0 ) {
307 0 : bi=-bi;
308 0 : len_bi=-len_bi;
309 0 : for(unsigned i=0; i<3; ++i) {
310 0 : dbi[i]*=-1.0;
311 : }
312 : }
313 1560 : if( len_cross<0 ) {
314 0 : cross=-cross;
315 0 : len_cross=-len_cross;
316 0 : for(unsigned i=0; i<3; ++i) {
317 0 : dcross[i]*=-1.0;
318 : }
319 : }
320 1560 : if( len_perp<0 ) {
321 0 : perp=-perp;
322 0 : len_perp=-len_perp;
323 0 : for(unsigned i=0; i<3; ++i) {
324 0 : dperp[i]*=-1.0;
325 : }
326 : }
327 1560 : if( len_bi<=0 || len_cross<=0 || len_bi<=0 ) {
328 0 : plumed_merror("Invalid box coordinates");
329 : }
330 :
331 : // Now derivatives of lengths
332 1560 : Tensor dd3( Tensor::identity() );
333 1560 : dlbi[0] = matmul(d3,dbi[0]) - matmul(bi,dd3);
334 1560 : dlbi[1] = matmul(d3,dbi[1]);
335 1560 : dlbi[2] = matmul(d3,dbi[2]);
336 1560 : dlbi[3] = matmul(bi,dd3);
337 :
338 1560 : dlcross[0] = matmul(d3,dcross[0]) - matmul(cross,dd3);
339 1560 : dlcross[1] = matmul(d3,dcross[1]);
340 1560 : dlcross[2] = matmul(d3,dcross[2]);
341 1560 : dlcross[3] = matmul(cross,dd3);
342 :
343 1560 : dlperp[0] = matmul(d3,dperp[0]) - matmul(perp,dd3);
344 1560 : dlperp[1] = matmul(d3,dperp[1]);
345 1560 : dlperp[2] = matmul(d3,dperp[2]);
346 1560 : dlperp[3] = matmul(perp,dd3);
347 :
348 : // Need to calculate the jacobian
349 : Tensor jacob;
350 1560 : jacob(0,0)=bi[0];
351 1560 : jacob(1,0)=bi[1];
352 1560 : jacob(2,0)=bi[2];
353 1560 : jacob(0,1)=cross[0];
354 1560 : jacob(1,1)=cross[1];
355 1560 : jacob(2,1)=cross[2];
356 1560 : jacob(0,2)=perp[0];
357 1560 : jacob(1,2)=perp[1];
358 1560 : jacob(2,2)=perp[2];
359 1560 : jacob_det = fabs( jacob.determinant() );
360 1560 : }
361 :
362 1560 : void VolumeCavity::calculateNumberInside( const VolumeInput& input, const VolumeCavity& actioninput, VolumeOutput& output ) {
363 : // Calculate distance of atom from origin of new coordinate frame
364 1560 : Vector datom=input.pbc.distance( Vector(input.refpos[0][0],input.refpos[0][1],input.refpos[0][2]),
365 1560 : Vector(input.cpos[0],input.cpos[1],input.cpos[2]) );
366 : double ucontr, uder, vcontr, vder, wcontr, wder;
367 :
368 : // Setup the histogram bead
369 1560 : HistogramBead bead( actioninput.kerneltype, 0, actioninput.len_bi, actioninput.sigma);
370 : // Calculate contribution from integral along bi
371 : double upos=dotProduct( datom, actioninput.bi );
372 1560 : ucontr=bead.calculate( upos, uder );
373 1560 : double udlen=bead.uboundDerivative( upos );
374 1560 : double uder2 = bead.lboundDerivative( upos ) - udlen;
375 :
376 : // Calculate contribution from integral along cross
377 1560 : bead.set( 0, actioninput.len_cross, actioninput.sigma );
378 : double vpos=dotProduct( datom, actioninput.cross );
379 1560 : vcontr=bead.calculate( vpos, vder );
380 1560 : double vdlen=bead.uboundDerivative( vpos );
381 1560 : double vder2 = bead.lboundDerivative( vpos ) - vdlen;
382 :
383 : // Calculate contribution from integral along perp
384 1560 : bead.set( 0, actioninput.len_perp, actioninput.sigma );
385 : double wpos=dotProduct( datom, actioninput.perp );
386 1560 : wcontr=bead.calculate( wpos, wder );
387 1560 : double wdlen=bead.uboundDerivative( wpos );
388 1560 : double wder2 = bead.lboundDerivative( wpos ) - wdlen;
389 :
390 : Vector dfd;
391 1560 : dfd[0]=uder*vcontr*wcontr;
392 1560 : dfd[1]=ucontr*vder*wcontr;
393 1560 : dfd[2]=ucontr*vcontr*wder;
394 1560 : output.derivatives[0] = (dfd[0]*actioninput.bi[0]+dfd[1]*actioninput.cross[0]+dfd[2]*actioninput.perp[0]);
395 1560 : output.derivatives[1] = (dfd[0]*actioninput.bi[1]+dfd[1]*actioninput.cross[1]+dfd[2]*actioninput.perp[1]);
396 1560 : output.derivatives[2] = (dfd[0]*actioninput.bi[2]+dfd[1]*actioninput.cross[2]+dfd[2]*actioninput.perp[2]);
397 1560 : output.values[0] = ucontr*vcontr*wcontr*actioninput.jacob_det;
398 :
399 : // Add reference atom derivatives
400 1560 : dfd[0]=uder2*vcontr*wcontr;
401 1560 : dfd[1]=ucontr*vder2*wcontr;
402 1560 : dfd[2]=ucontr*vcontr*wder2;
403 : Vector dfld;
404 1560 : dfld[0]=udlen*vcontr*wcontr;
405 1560 : dfld[1]=ucontr*vdlen*wcontr;
406 1560 : dfld[2]=ucontr*vcontr*wdlen;
407 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]) +
408 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]);
409 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]) +
410 : dfld[0]*actioninput.dlbi[1] + dfld[1]*actioninput.dlcross[1] + dfld[2]*actioninput.dlperp[1];
411 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]) +
412 : dfld[0]*actioninput.dlbi[2] + dfld[1]*actioninput.dlcross[2] + dfld[2]*actioninput.dlperp[2];
413 : output.refders[3] = dfld[0]*actioninput.dlbi[3] + dfld[1]*actioninput.dlcross[3] + dfld[2]*actioninput.dlperp[3];
414 :
415 : Tensor vir;
416 1560 : vir=-Tensor( Vector(input.cpos[0],input.cpos[1],input.cpos[2]), Vector(output.derivatives[0],output.derivatives[1],output.derivatives[2]) );
417 7800 : for(unsigned i=0; i<4; ++i) {
418 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]) );
419 : }
420 1560 : output.virial.set( 0, vir );
421 1560 : }
422 :
423 : }
424 : }
|
