src/gmxlib/ewald_util.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdio.h>
  40 #include <math.h>
  41 #include "maths.h"
  42 #include "typedefs.h"
  43 #include "vec.h"
  44 #include "coulomb.h"
  45 #include "smalloc.h"
  46 #include "physics.h"
  47 #include "txtdump.h"
  48 #include "futil.h"
  49 #include "names.h"
  50 #include "writeps.h"
  51 #include "macros.h"
  52
  53 real calc_ewaldcoeff(real rc,real dtol)
  54 {
  55   real x=5,low,high;
  56   int n,i=0;
  57
  58
  59   do {
  60     i++;
  61     x*=2;
  62   } while (gmx_erfc(x*rc) > dtol);
  63
  64   n=i+60; /* search tolerance is 2^-60 */
  65   low=0;
  66   high=x;
  67   for(i=0;i<n;i++) {
  68     x=(low+high)/2;
  69     if (gmx_erfc(x*rc) > dtol)
  70       low=x;
  71     else
  72       high=x;
  73   }
  74   return x;
  75 }
  76
  77
  78
  79 real ewald_LRcorrection(FILE *fplog,
  80                         int start,int end,
  81                         t_commrec *cr,t_forcerec *fr,
  82                         real *chargeA,real *chargeB,
  83                         t_blocka *excl,rvec x[],
  84                         matrix box,rvec mu_tot[],
  85                         int ewald_geometry,real epsilon_surface,
  86                         real lambda,real *dvdlambda,
  87                         real *vdip,real *vcharge)
  88 {
  89   int     i,i1,i2,j,k,m,iv,jv,q;
  90   atom_id *AA;
  91   double  q2sumA,q2sumB,Vexcl,dvdl_excl; /* Necessary for precision */
  92   real    one_4pi_eps;
  93   real    v,vc,qiA,qiB,dr,dr2,rinv,fscal,enercorr;
  94   real    VselfA,VselfB=0,Vcharge[2],Vdipole[2],rinv2,ewc=fr->ewaldcoeff,ewcdr;
  95   rvec    df,dx,mutot[2],dipcorrA,dipcorrB;
  96   rvec    *f=fr->f_novirsum;
  97   tensor  dxdf;
  98   real    vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
  99   real    L1,dipole_coeff,qqA,qqB,qqL,vr0;
 100   /*#define TABLES*/
 101 #ifdef TABLES
 102   real    tabscale=fr->tabscale;
 103   real    eps,eps2,VV,FF,F,Y,Geps,Heps2,Fp,fijC,r1t;
 104   real    *VFtab=fr->coulvdwtab;
 105   int     n0,n1,nnn;
 106 #else
 107   double  isp=0.564189583547756;
 108 #endif
 109   int     niat;
 110   gmx_bool    bFreeEnergy = (chargeB != NULL);
 111   gmx_bool    bMolPBC = fr->bMolPBC;
 112
 113   one_4pi_eps = ONE_4PI_EPS0/fr->epsilon_r;
 114   vr0 = ewc*2/sqrt(M_PI);
 115
 116   AA         = excl->a;
 117   Vexcl      = 0;
 118   dvdl_excl  = 0;
 119   q2sumA     = 0;
 120   q2sumB     = 0;
 121   Vdipole[0] = 0;
 122   Vdipole[1] = 0;
 123   Vcharge[0] = 0;
 124   Vcharge[1] = 0;
 125   L1         = 1.0-lambda;
 126
 127   /* Note that we have to transform back to gromacs units, since
 128    * mu_tot contains the dipole in debye units (for output).
 129    */
 130   for(i=0; (i<DIM); i++) {
 131     mutot[0][i] = mu_tot[0][i]*DEBYE2ENM;
 132     mutot[1][i] = mu_tot[1][i]*DEBYE2ENM;
 133     dipcorrA[i] = 0;
 134     dipcorrB[i] = 0;
 135   }
 136   dipole_coeff=0;
 137   switch (ewald_geometry) {
 138   case eewg3D:
 139     if (epsilon_surface != 0) {
 140       dipole_coeff =
 141         2*M_PI*ONE_4PI_EPS0/((2*epsilon_surface + fr->epsilon_r)*vol);
 142       for(i=0; (i<DIM); i++) {
 143         dipcorrA[i] = 2*dipole_coeff*mutot[0][i];
 144         dipcorrB[i] = 2*dipole_coeff*mutot[1][i];
 145       }
 146     }
 147     break;
 148   case eewg3DC:
 149     dipole_coeff = 2*M_PI*one_4pi_eps/vol;
 150     dipcorrA[ZZ] = 2*dipole_coeff*mutot[0][ZZ];
 151     dipcorrB[ZZ] = 2*dipole_coeff*mutot[1][ZZ];
 152     break;
 153   default:
 154     gmx_incons("Unsupported Ewald geometry");
 155     break;
 156   }
 157   if (debug) {
 158     fprintf(debug,"dipcorr = %8.3f  %8.3f  %8.3f\n",
 159             dipcorrA[XX],dipcorrA[YY],dipcorrA[ZZ]);
 160     fprintf(debug,"mutot   = %8.3f  %8.3f  %8.3f\n",
 161             mutot[0][XX],mutot[0][YY],mutot[0][ZZ]);
 162   }
 163
 164   if (DOMAINDECOMP(cr))
 165     niat = excl->nr;
 166   else
 167     niat = end;
 168
 169   clear_mat(dxdf);
 170   if (!bFreeEnergy) {
 171     for(i=start; (i<niat); i++) {
 172       /* Initiate local variables (for this i-particle) to 0 */
 173       qiA = chargeA[i]*one_4pi_eps;
 174       i1  = excl->index[i];
 175       i2  = excl->index[i+1];
 176       if (i < end)
 177         q2sumA += chargeA[i]*chargeA[i];
 178
 179       /* Loop over excluded neighbours */
 180       for(j=i1; (j<i2); j++) {
 181         k = AA[j];
 182         /*
 183          * First we must test whether k <> i, and then, because the
 184          * exclusions are all listed twice i->k and k->i we must select
 185          * just one of the two.
 186          * As a minor optimization we only compute forces when the charges
 187          * are non-zero.
 188          */
 189         if (k > i) {
 190           qqA = qiA*chargeA[k];
 191           if (qqA != 0.0) {
 192             rvec_sub(x[i],x[k],dx);
 193             if (bMolPBC) {
 194               /* Cheap pbc_dx, assume excluded pairs are at short distance. */
 195               for(m=DIM-1; (m>=0); m--) {
 196                 if (dx[m] > 0.5*box[m][m])
 197                   rvec_dec(dx,box[m]);
 198                 else if (dx[m] < -0.5*box[m][m])
 199                   rvec_inc(dx,box[m]);
 200               }
 201             }
 202             dr2 = norm2(dx);
 203             /* Distance between two excluded particles may be zero in the
 204              * case of shells
 205              */
 206             if (dr2 != 0) {
 207               rinv              = gmx_invsqrt(dr2);
 208               rinv2             = rinv*rinv;
 209               dr                = 1.0/rinv;
 210 #ifdef TABLES
 211               r1t               = tabscale*dr;
 212               n0                = r1t;
 213               assert(n0 >= 3);
 214               n1                = 12*n0;
 215               eps               = r1t-n0;
 216               eps2              = eps*eps;
 217               nnn               = n1;
 218               Y                 = VFtab[nnn];
 219               F                 = VFtab[nnn+1];
 220               Geps              = eps*VFtab[nnn+2];
 221               Heps2             = eps2*VFtab[nnn+3];
 222               Fp                = F+Geps+Heps2;
 223               VV                = Y+eps*Fp;
 224               FF                = Fp+Geps+2.0*Heps2;
 225               vc                = qqA*(rinv-VV);
 226               fijC              = qqA*FF;
 227               Vexcl            += vc;
 228
 229               fscal             = vc*rinv2+fijC*tabscale*rinv;
 230             /* End of tabulated interaction part */
 231 #else
 232
 233             /* This is the code you would want instead if not using
 234              * tables:
 235              */
 236               ewcdr   = ewc*dr;
 237               vc      = qqA*gmx_erf(ewcdr)*rinv;
 238               Vexcl  += vc;
 239 #ifdef GMX_DOUBLE
 240               /* Relative accuracy at R_ERF_R_INACC of 3e-10 */
 241 #define       R_ERF_R_INACC 0.006
 242 #else
 243               /* Relative accuracy at R_ERF_R_INACC of 2e-5 */
 244 #define       R_ERF_R_INACC 0.1
 245 #endif
 246               if (ewcdr > R_ERF_R_INACC) {
 247                 fscal = rinv2*(vc - 2.0*qqA*ewc*isp*exp(-ewcdr*ewcdr));
 248               } else {
 249                 /* Use a fourth order series expansion for small ewcdr */
 250                 fscal = ewc*ewc*qqA*vr0*(2.0/3.0 - 0.4*ewcdr*ewcdr);
 251               }
 252 #endif
 253               /* The force vector is obtained by multiplication with the
 254                * distance vector
 255                */
 256               svmul(fscal,dx,df);
 257               rvec_inc(f[k],df);
 258               rvec_dec(f[i],df);
 259               for(iv=0; (iv<DIM); iv++)
 260                 for(jv=0; (jv<DIM); jv++)
 261                   dxdf[iv][jv] += dx[iv]*df[jv];
 262             } else {
 263               Vexcl += qqA*vr0;
 264             }
 265           }
 266         }
 267       }
 268       /* Dipole correction on force */
 269       if (dipole_coeff != 0) {
 270         for(j=0; (j<DIM); j++)
 271           f[i][j] -= dipcorrA[j]*chargeA[i];
 272       }
 273     }
 274   } else {
 275     for(i=start; (i<niat); i++) {
 276       /* Initiate local variables (for this i-particle) to 0 */
 277       qiA = chargeA[i]*one_4pi_eps;
 278       qiB = chargeB[i]*one_4pi_eps;
 279       i1  = excl->index[i];
 280       i2  = excl->index[i+1];
 281       if (i < end) {
 282         q2sumA += chargeA[i]*chargeA[i];
 283         q2sumB += chargeB[i]*chargeB[i];
 284       }
 285
 286       /* Loop over excluded neighbours */
 287       for(j=i1; (j<i2); j++) {
 288         k = AA[j];
 289         if (k > i) {
 290           qqA = qiA*chargeA[k];
 291           qqB = qiB*chargeB[k];
 292           if (qqA != 0.0 || qqB != 0.0) {
 293             qqL = L1*qqA + lambda*qqB;
 294             rvec_sub(x[i],x[k],dx);
 295             if (bMolPBC) {
 296               /* Cheap pbc_dx, assume excluded pairs are at short distance. */
 297               for(m=DIM-1; (m>=0); m--) {
 298                 if (dx[m] > 0.5*box[m][m])
 299                   rvec_dec(dx,box[m]);
 300                 else if (dx[m] < -0.5*box[m][m])
 301                   rvec_inc(dx,box[m]);
 302               }
 303             }
 304             dr2 = norm2(dx);
 305             if (dr2 != 0) {
 306               rinv   = gmx_invsqrt(dr2);
 307               rinv2  = rinv*rinv;
 308               dr     = 1.0/rinv;
 309               v      = gmx_erf(ewc*dr)*rinv;
 310               vc     = qqL*v;
 311               Vexcl += vc;
 312               fscal  = rinv2*(vc-2.0*qqL*ewc*isp*exp(-ewc*ewc*dr2));
 313               svmul(fscal,dx,df);
 314               rvec_inc(f[k],df);
 315               rvec_dec(f[i],df);
 316               for(iv=0; (iv<DIM); iv++)
 317                 for(jv=0; (jv<DIM); jv++)
 318                   dxdf[iv][jv] += dx[iv]*df[jv];
 319               dvdl_excl += (qqB - qqA)*v;
 320             } else {
 321               Vexcl     +=         qqL*vr0;
 322               dvdl_excl += (qqB - qqA)*vr0;
 323             }
 324           }
 325         }
 326       }
 327       /* Dipole correction on force */
 328       if (dipole_coeff != 0) {
 329         for(j=0; (j<DIM); j++)
 330           f[i][j] -= L1*dipcorrA[j]*chargeA[i]
 331             + lambda*dipcorrB[j]*chargeB[i];
 332       }
 333     }
 334   }
 335   for(iv=0; (iv<DIM); iv++)
 336     for(jv=0; (jv<DIM); jv++)
 337       fr->vir_el_recip[iv][jv] += 0.5*dxdf[iv][jv];
 338
 339   /* Global corrections only on master process */
 340   if (MASTER(cr)) {
 341     for(q=0; q<(bFreeEnergy ? 2 : 1); q++) {
 342       /* Apply charge correction */
 343       /* use vc as a dummy variable */
 344       vc = fr->qsum[q]*fr->qsum[q]*M_PI*one_4pi_eps/(2.0*vol*vol*ewc*ewc);
 345       for(iv=0; (iv<DIM); iv++)
 346         fr->vir_el_recip[iv][iv] +=
 347           (bFreeEnergy ? (q==0 ? L1*vc : lambda*vc) : vc);
 348       Vcharge[q] = -vol*vc;
 349
 350       /* Apply surface dipole correction:
 351        * correction = dipole_coeff * (dipole)^2
 352        */
 353       if (dipole_coeff != 0) {
 354         if (ewald_geometry == eewg3D)
 355           Vdipole[q] = dipole_coeff*iprod(mutot[q],mutot[q]);
 356         else if (ewald_geometry == eewg3DC)
 357           Vdipole[q] = dipole_coeff*mutot[q][ZZ]*mutot[q][ZZ];
 358       }
 359     }
 360   }
 361
 362   VselfA = ewc*one_4pi_eps*q2sumA/sqrt(M_PI);
 363
 364   if (!bFreeEnergy) {
 365     *vcharge = Vcharge[0];
 366     *vdip    = Vdipole[0];
 367     enercorr = *vcharge + *vdip - VselfA - Vexcl;
 368    } else {
 369     VselfB = ewc*one_4pi_eps*q2sumB/sqrt(M_PI);
 370     *vcharge = L1*Vcharge[0] + lambda*Vcharge[1];
 371     *vdip    = L1*Vdipole[0] + lambda*Vdipole[1];
 372     enercorr = *vcharge + *vdip - (L1*VselfA + lambda*VselfB) - Vexcl;
 373     *dvdlambda += Vdipole[1] + Vcharge[1] - VselfB
 374       - (Vdipole[0] + Vcharge[0] - VselfA) - dvdl_excl;
 375   }
 376
 377   if (debug) {
 378     fprintf(debug,"Long Range corrections for Ewald interactions:\n");
 379     fprintf(debug,"start=%d,natoms=%d\n",start,end-start);
 380     fprintf(debug,"q2sum = %g, Vself=%g\n",
 381             L1*q2sumA+lambda*q2sumB,L1*VselfA+lambda*VselfB);
 382     fprintf(debug,"Long Range correction: Vexcl=%g\n",Vexcl);
 383     if (MASTER(cr)) {
 384       fprintf(debug,"Total charge correction: Vcharge=%g\n",
 385               L1*Vcharge[0]+lambda*Vcharge[1]);
 386       if (epsilon_surface > 0 || ewald_geometry == eewg3DC) {
 387         fprintf(debug,"Total dipole correction: Vdipole=%g\n",
 388                 L1*Vdipole[0]+lambda*Vdipole[1]);
 389       }
 390     }
 391   }
 392
 393   /* Return the correction to the energy */
 394   return enercorr;
 395 }