2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
56 real calc_ewaldcoeff(real rc,real dtol)
65 } while (gmx_erfc(x*rc) > dtol);
67 n=i+60; /* search tolerance is 2^-60 */
72 if (gmx_erfc(x*rc) > dtol)
82 real ewald_LRcorrection(FILE *fplog,
84 t_commrec *cr,int thread,t_forcerec *fr,
85 real *chargeA,real *chargeB,
86 gmx_bool calc_excl_corr,
87 t_blocka *excl,rvec x[],
88 matrix box,rvec mu_tot[],
89 int ewald_geometry,real epsilon_surface,
91 real lambda,real *dvdlambda)
93 int i,i1,i2,j,k,m,iv,jv,q;
95 double q2sumA,q2sumB,Vexcl,dvdl_excl; /* Necessary for precision */
97 real v,vc,qiA,qiB,dr,dr2,rinv,fscal,enercorr;
98 real Vself[2],Vdipole[2],rinv2,ewc=fr->ewaldcoeff,ewcdr;
99 rvec df,dx,mutot[2],dipcorrA,dipcorrB;
101 real vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
102 real L1,dipole_coeff,qqA,qqB,qqL,vr0;
105 real tabscale=fr->tabscale;
106 real eps,eps2,VV,FF,F,Y,Geps,Heps2,Fp,fijC,r1t;
107 real *VFtab=fr->coulvdwtab;
110 gmx_bool bFreeEnergy = (chargeB != NULL);
111 gmx_bool bMolPBC = fr->bMolPBC;
113 one_4pi_eps = ONE_4PI_EPS0/fr->epsilon_r;
114 vr0 = ewc*M_2_SQRTPI;
125 /* Note that we have to transform back to gromacs units, since
126 * mu_tot contains the dipole in debye units (for output).
128 for(i=0; (i<DIM); i++) {
129 mutot[0][i] = mu_tot[0][i]*DEBYE2ENM;
130 mutot[1][i] = mu_tot[1][i]*DEBYE2ENM;
135 switch (ewald_geometry) {
137 if (epsilon_surface != 0) {
139 2*M_PI*ONE_4PI_EPS0/((2*epsilon_surface + fr->epsilon_r)*vol);
140 for(i=0; (i<DIM); i++) {
141 dipcorrA[i] = 2*dipole_coeff*mutot[0][i];
142 dipcorrB[i] = 2*dipole_coeff*mutot[1][i];
147 dipole_coeff = 2*M_PI*one_4pi_eps/vol;
148 dipcorrA[ZZ] = 2*dipole_coeff*mutot[0][ZZ];
149 dipcorrB[ZZ] = 2*dipole_coeff*mutot[1][ZZ];
152 gmx_incons("Unsupported Ewald geometry");
156 fprintf(debug,"dipcorr = %8.3f %8.3f %8.3f\n",
157 dipcorrA[XX],dipcorrA[YY],dipcorrA[ZZ]);
158 fprintf(debug,"mutot = %8.3f %8.3f %8.3f\n",
159 mutot[0][XX],mutot[0][YY],mutot[0][ZZ]);
163 if ((calc_excl_corr || dipole_coeff != 0) && !bFreeEnergy) {
164 for(i=start; (i<end); i++) {
165 /* Initiate local variables (for this i-particle) to 0 */
166 qiA = chargeA[i]*one_4pi_eps;
171 i2 = excl->index[i+1];
173 /* Loop over excluded neighbours */
174 for(j=i1; (j<i2); j++) {
177 * First we must test whether k <> i, and then, because the
178 * exclusions are all listed twice i->k and k->i we must select
179 * just one of the two.
180 * As a minor optimization we only compute forces when the charges
184 qqA = qiA*chargeA[k];
186 rvec_sub(x[i],x[k],dx);
188 /* Cheap pbc_dx, assume excluded pairs are at short distance. */
189 for(m=DIM-1; (m>=0); m--) {
190 if (dx[m] > 0.5*box[m][m])
192 else if (dx[m] < -0.5*box[m][m])
197 /* Distance between two excluded particles may be zero in the
201 rinv = gmx_invsqrt(dr2);
214 Geps = eps*VFtab[nnn+2];
215 Heps2 = eps2*VFtab[nnn+3];
218 FF = Fp+Geps+2.0*Heps2;
223 fscal = vc*rinv2+fijC*tabscale*rinv;
224 /* End of tabulated interaction part */
227 /* This is the code you would want instead if not using
231 vc = qqA*gmx_erf(ewcdr)*rinv;
234 /* Relative accuracy at R_ERF_R_INACC of 3e-10 */
235 #define R_ERF_R_INACC 0.006
237 /* Relative accuracy at R_ERF_R_INACC of 2e-5 */
238 #define R_ERF_R_INACC 0.1
240 if (ewcdr > R_ERF_R_INACC) {
241 fscal = rinv2*(vc - qqA*ewc*M_2_SQRTPI*exp(-ewcdr*ewcdr));
243 /* Use a fourth order series expansion for small ewcdr */
244 fscal = ewc*ewc*qqA*vr0*(2.0/3.0 - 0.4*ewcdr*ewcdr);
247 /* The force vector is obtained by multiplication with the
253 for(iv=0; (iv<DIM); iv++)
254 for(jv=0; (jv<DIM); jv++)
255 dxdf[iv][jv] += dx[iv]*df[jv];
263 /* Dipole correction on force */
264 if (dipole_coeff != 0) {
265 for(j=0; (j<DIM); j++)
266 f[i][j] -= dipcorrA[j]*chargeA[i];
269 } else if (calc_excl_corr || dipole_coeff != 0) {
270 for(i=start; (i<end); i++) {
271 /* Initiate local variables (for this i-particle) to 0 */
272 qiA = chargeA[i]*one_4pi_eps;
273 qiB = chargeB[i]*one_4pi_eps;
278 i2 = excl->index[i+1];
280 /* Loop over excluded neighbours */
281 for(j=i1; (j<i2); j++) {
284 qqA = qiA*chargeA[k];
285 qqB = qiB*chargeB[k];
286 if (qqA != 0.0 || qqB != 0.0) {
287 qqL = L1*qqA + lambda*qqB;
288 rvec_sub(x[i],x[k],dx);
290 /* Cheap pbc_dx, assume excluded pairs are at short distance. */
291 for(m=DIM-1; (m>=0); m--) {
292 if (dx[m] > 0.5*box[m][m])
294 else if (dx[m] < -0.5*box[m][m])
300 rinv = gmx_invsqrt(dr2);
303 v = gmx_erf(ewc*dr)*rinv;
306 fscal = rinv2*(vc-qqL*ewc*M_2_SQRTPI*exp(-ewc*ewc*dr2));
310 for(iv=0; (iv<DIM); iv++)
311 for(jv=0; (jv<DIM); jv++)
312 dxdf[iv][jv] += dx[iv]*df[jv];
313 dvdl_excl += (qqB - qqA)*v;
316 dvdl_excl += (qqB - qqA)*vr0;
322 /* Dipole correction on force */
323 if (dipole_coeff != 0) {
324 for(j=0; (j<DIM); j++)
325 f[i][j] -= L1*dipcorrA[j]*chargeA[i]
326 + lambda*dipcorrB[j]*chargeB[i];
330 for(iv=0; (iv<DIM); iv++)
331 for(jv=0; (jv<DIM); jv++)
332 vir[iv][jv] += 0.5*dxdf[iv][jv];
337 /* Global corrections only on master process */
338 if (MASTER(cr) && thread == 0) {
339 for(q=0; q<(bFreeEnergy ? 2 : 1); q++) {
340 if (calc_excl_corr) {
341 /* Self-energy correction */
342 Vself[q] = ewc*one_4pi_eps*fr->q2sum[q]*M_1_SQRTPI;
345 /* Apply surface dipole correction:
346 * correction = dipole_coeff * (dipole)^2
348 if (dipole_coeff != 0) {
349 if (ewald_geometry == eewg3D)
350 Vdipole[q] = dipole_coeff*iprod(mutot[q],mutot[q]);
351 else if (ewald_geometry == eewg3DC)
352 Vdipole[q] = dipole_coeff*mutot[q][ZZ]*mutot[q][ZZ];
358 enercorr = Vdipole[0] - Vself[0] - Vexcl;
360 enercorr = L1*(Vdipole[0] - Vself[0])
361 + lambda*(Vdipole[1] - Vself[1])
363 *dvdlambda += Vdipole[1] - Vself[1]
364 - (Vdipole[0] - Vself[0]) - dvdl_excl;
368 fprintf(debug,"Long Range corrections for Ewald interactions:\n");
369 fprintf(debug,"start=%d,natoms=%d\n",start,end-start);
370 fprintf(debug,"q2sum = %g, Vself=%g\n",
371 L1*q2sumA+lambda*q2sumB,L1*Vself[0]+lambda*Vself[1]);
372 fprintf(debug,"Long Range correction: Vexcl=%g\n",Vexcl);
373 if (MASTER(cr) && thread == 0) {
374 if (epsilon_surface > 0 || ewald_geometry == eewg3DC) {
375 fprintf(debug,"Total dipole correction: Vdipole=%g\n",
376 L1*Vdipole[0]+lambda*Vdipole[1]);
381 /* Return the correction to the energy */
385 real ewald_charge_correction(t_commrec *cr,t_forcerec *fr,real lambda,
387 real *dvdlambda,tensor vir)
390 real vol,fac,qs2A,qs2B,vc,enercorr;
395 /* Apply charge correction */
396 vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
398 fac = M_PI*ONE_4PI_EPS0/(fr->epsilon_r*2.0*vol*vol*sqr(fr->ewaldcoeff));
400 qs2A = fr->qsum[0]*fr->qsum[0];
401 qs2B = fr->qsum[1]*fr->qsum[1];
403 vc = (qs2A*(1 - lambda) + qs2B*lambda)*fac;
407 *dvdlambda += -vol*(qs2B - qs2A)*fac;
416 fprintf(debug,"Total charge correction: Vcharge=%g\n",enercorr);