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,int thread,t_forcerec *fr,
  82                         real *chargeA,real *chargeB,
  83                         gmx_bool calc_excl_corr,
  84                         t_blocka *excl,rvec x[],
  85                         matrix box,rvec mu_tot[],
  86                         int ewald_geometry,real epsilon_surface,
  87                         rvec *f,tensor vir,
  88                         real lambda,real *dvdlambda)
  89 {
  90   int     i,i1,i2,j,k,m,iv,jv,q;
  91   atom_id *AA;
  92   double  q2sumA,q2sumB,Vexcl,dvdl_excl; /* Necessary for precision */
  93   real    one_4pi_eps;
  94   real    v,vc,qiA,qiB,dr,dr2,rinv,fscal,enercorr;
  95   real    Vself[2],Vdipole[2],rinv2,ewc=fr->ewaldcoeff,ewcdr;
  96   rvec    df,dx,mutot[2],dipcorrA,dipcorrB;
  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   gmx_bool    bFreeEnergy = (chargeB != NULL);
 110   gmx_bool    bMolPBC = fr->bMolPBC;
 111
 112   one_4pi_eps = ONE_4PI_EPS0/fr->epsilon_r;
 113   vr0 = ewc*2/sqrt(M_PI);
 114
 115   AA         = excl->a;
 116   Vexcl      = 0;
 117   dvdl_excl  = 0;
 118   q2sumA     = 0;
 119   q2sumB     = 0;
 120   Vdipole[0] = 0;
 121   Vdipole[1] = 0;
 122   L1         = 1.0-lambda;
 123
 124   /* Note that we have to transform back to gromacs units, since
 125    * mu_tot contains the dipole in debye units (for output).
 126    */
 127   for(i=0; (i<DIM); i++) {
 128     mutot[0][i] = mu_tot[0][i]*DEBYE2ENM;
 129     mutot[1][i] = mu_tot[1][i]*DEBYE2ENM;
 130     dipcorrA[i] = 0;
 131     dipcorrB[i] = 0;
 132   }
 133   dipole_coeff=0;
 134   switch (ewald_geometry) {
 135   case eewg3D:
 136     if (epsilon_surface != 0) {
 137       dipole_coeff =
 138         2*M_PI*ONE_4PI_EPS0/((2*epsilon_surface + fr->epsilon_r)*vol);
 139       for(i=0; (i<DIM); i++) {
 140         dipcorrA[i] = 2*dipole_coeff*mutot[0][i];
 141         dipcorrB[i] = 2*dipole_coeff*mutot[1][i];
 142       }
 143     }
 144     break;
 145   case eewg3DC:
 146     dipole_coeff = 2*M_PI*one_4pi_eps/vol;
 147     dipcorrA[ZZ] = 2*dipole_coeff*mutot[0][ZZ];
 148     dipcorrB[ZZ] = 2*dipole_coeff*mutot[1][ZZ];
 149     break;
 150   default:
 151     gmx_incons("Unsupported Ewald geometry");
 152     break;
 153   }
 154   if (debug) {
 155     fprintf(debug,"dipcorr = %8.3f  %8.3f  %8.3f\n",
 156             dipcorrA[XX],dipcorrA[YY],dipcorrA[ZZ]);
 157     fprintf(debug,"mutot   = %8.3f  %8.3f  %8.3f\n",
 158             mutot[0][XX],mutot[0][YY],mutot[0][ZZ]);
 159   }
 160
 161   clear_mat(dxdf);
 162   if ((calc_excl_corr || dipole_coeff != 0) && !bFreeEnergy) {
 163     for(i=start; (i<end); i++) {
 164       /* Initiate local variables (for this i-particle) to 0 */
 165       qiA = chargeA[i]*one_4pi_eps;
 166
 167       if (calc_excl_corr)
 168       {
 169       i1  = excl->index[i];
 170       i2  = excl->index[i+1];
 171
 172       /* Loop over excluded neighbours */
 173       for(j=i1; (j<i2); j++) {
 174         k = AA[j];
 175         /*
 176          * First we must test whether k <> i, and then, because the
 177          * exclusions are all listed twice i->k and k->i we must select
 178          * just one of the two.
 179          * As a minor optimization we only compute forces when the charges
 180          * are non-zero.
 181          */
 182         if (k > i) {
 183           qqA = qiA*chargeA[k];
 184           if (qqA != 0.0) {
 185             rvec_sub(x[i],x[k],dx);
 186             if (bMolPBC) {
 187               /* Cheap pbc_dx, assume excluded pairs are at short distance. */
 188               for(m=DIM-1; (m>=0); m--) {
 189                 if (dx[m] > 0.5*box[m][m])
 190                   rvec_dec(dx,box[m]);
 191                 else if (dx[m] < -0.5*box[m][m])
 192                   rvec_inc(dx,box[m]);
 193               }
 194             }
 195             dr2 = norm2(dx);
 196             /* Distance between two excluded particles may be zero in the
 197              * case of shells
 198              */
 199             if (dr2 != 0) {
 200               rinv              = gmx_invsqrt(dr2);
 201               rinv2             = rinv*rinv;
 202               dr                = 1.0/rinv;
 203 #ifdef TABLES
 204               r1t               = tabscale*dr;
 205               n0                = r1t;
 206               assert(n0 >= 3);
 207               n1                = 12*n0;
 208               eps               = r1t-n0;
 209               eps2              = eps*eps;
 210               nnn               = n1;
 211               Y                 = VFtab[nnn];
 212               F                 = VFtab[nnn+1];
 213               Geps              = eps*VFtab[nnn+2];
 214               Heps2             = eps2*VFtab[nnn+3];
 215               Fp                = F+Geps+Heps2;
 216               VV                = Y+eps*Fp;
 217               FF                = Fp+Geps+2.0*Heps2;
 218               vc                = qqA*(rinv-VV);
 219               fijC              = qqA*FF;
 220               Vexcl            += vc;
 221
 222               fscal             = vc*rinv2+fijC*tabscale*rinv;
 223             /* End of tabulated interaction part */
 224 #else
 225
 226             /* This is the code you would want instead if not using
 227              * tables:
 228              */
 229               ewcdr   = ewc*dr;
 230               vc      = qqA*gmx_erf(ewcdr)*rinv;
 231               Vexcl  += vc;
 232 #ifdef GMX_DOUBLE
 233               /* Relative accuracy at R_ERF_R_INACC of 3e-10 */
 234 #define       R_ERF_R_INACC 0.006
 235 #else
 236               /* Relative accuracy at R_ERF_R_INACC of 2e-5 */
 237 #define       R_ERF_R_INACC 0.1
 238 #endif
 239               if (ewcdr > R_ERF_R_INACC) {
 240                 fscal = rinv2*(vc - 2.0*qqA*ewc*isp*exp(-ewcdr*ewcdr));
 241               } else {
 242                 /* Use a fourth order series expansion for small ewcdr */
 243                 fscal = ewc*ewc*qqA*vr0*(2.0/3.0 - 0.4*ewcdr*ewcdr);
 244               }
 245 #endif
 246               /* The force vector is obtained by multiplication with the
 247                * distance vector
 248                */
 249               svmul(fscal,dx,df);
 250               rvec_inc(f[k],df);
 251               rvec_dec(f[i],df);
 252               for(iv=0; (iv<DIM); iv++)
 253                 for(jv=0; (jv<DIM); jv++)
 254                   dxdf[iv][jv] += dx[iv]*df[jv];
 255             } else {
 256               Vexcl += qqA*vr0;
 257             }
 258           }
 259         }
 260       }
 261       }
 262       /* Dipole correction on force */
 263       if (dipole_coeff != 0) {
 264         for(j=0; (j<DIM); j++)
 265           f[i][j] -= dipcorrA[j]*chargeA[i];
 266       }
 267     }
 268   } else if (calc_excl_corr || dipole_coeff != 0) {
 269     for(i=start; (i<end); i++) {
 270       /* Initiate local variables (for this i-particle) to 0 */
 271       qiA = chargeA[i]*one_4pi_eps;
 272       qiB = chargeB[i]*one_4pi_eps;
 273
 274       if (calc_excl_corr)
 275       {
 276       i1  = excl->index[i];
 277       i2  = excl->index[i+1];
 278
 279       /* Loop over excluded neighbours */
 280       for(j=i1; (j<i2); j++) {
 281         k = AA[j];
 282         if (k > i) {
 283           qqA = qiA*chargeA[k];
 284           qqB = qiB*chargeB[k];
 285           if (qqA != 0.0 || qqB != 0.0) {
 286             qqL = L1*qqA + lambda*qqB;
 287             rvec_sub(x[i],x[k],dx);
 288             if (bMolPBC) {
 289               /* Cheap pbc_dx, assume excluded pairs are at short distance. */
 290               for(m=DIM-1; (m>=0); m--) {
 291                 if (dx[m] > 0.5*box[m][m])
 292                   rvec_dec(dx,box[m]);
 293                 else if (dx[m] < -0.5*box[m][m])
 294                   rvec_inc(dx,box[m]);
 295               }
 296             }
 297             dr2 = norm2(dx);
 298             if (dr2 != 0) {
 299               rinv   = gmx_invsqrt(dr2);
 300               rinv2  = rinv*rinv;
 301               dr     = 1.0/rinv;
 302               v      = gmx_erf(ewc*dr)*rinv;
 303               vc     = qqL*v;
 304               Vexcl += vc;
 305               fscal  = rinv2*(vc-2.0*qqL*ewc*isp*exp(-ewc*ewc*dr2));
 306               svmul(fscal,dx,df);
 307               rvec_inc(f[k],df);
 308               rvec_dec(f[i],df);
 309               for(iv=0; (iv<DIM); iv++)
 310                 for(jv=0; (jv<DIM); jv++)
 311                   dxdf[iv][jv] += dx[iv]*df[jv];
 312               dvdl_excl += (qqB - qqA)*v;
 313             } else {
 314               Vexcl     +=         qqL*vr0;
 315               dvdl_excl += (qqB - qqA)*vr0;
 316             }
 317           }
 318         }
 319       }
 320       }
 321       /* Dipole correction on force */
 322       if (dipole_coeff != 0) {
 323         for(j=0; (j<DIM); j++)
 324           f[i][j] -= L1*dipcorrA[j]*chargeA[i]
 325             + lambda*dipcorrB[j]*chargeB[i];
 326       }
 327     }
 328   }
 329   for(iv=0; (iv<DIM); iv++)
 330     for(jv=0; (jv<DIM); jv++)
 331       vir[iv][jv] += 0.5*dxdf[iv][jv];
 332
 333
 334   Vself[0] = 0;
 335   Vself[1] = 0;
 336   /* Global corrections only on master process */
 337   if (MASTER(cr) && thread == 0) {
 338     for(q=0; q<(bFreeEnergy ? 2 : 1); q++) {
 339       if (calc_excl_corr) {
 340         /* Self-energy correction */
 341         Vself[q] = ewc*one_4pi_eps*fr->q2sum[q]/sqrt(M_PI);
 342       }
 343
 344       /* Apply surface dipole correction:
 345        * correction = dipole_coeff * (dipole)^2
 346        */
 347       if (dipole_coeff != 0) {
 348         if (ewald_geometry == eewg3D)
 349           Vdipole[q] = dipole_coeff*iprod(mutot[q],mutot[q]);
 350         else if (ewald_geometry == eewg3DC)
 351           Vdipole[q] = dipole_coeff*mutot[q][ZZ]*mutot[q][ZZ];
 352       }
 353     }
 354   }
 355
 356   if (!bFreeEnergy) {
 357     enercorr = Vdipole[0] - Vself[0] - Vexcl;
 358    } else {
 359     enercorr = L1*(Vdipole[0] - Vself[0])
 360       + lambda*(Vdipole[1] - Vself[1])
 361       - Vexcl;
 362     *dvdlambda += Vdipole[1] - Vself[1]
 363       - (Vdipole[0] - Vself[0]) - dvdl_excl;
 364   }
 365
 366   if (debug) {
 367     fprintf(debug,"Long Range corrections for Ewald interactions:\n");
 368     fprintf(debug,"start=%d,natoms=%d\n",start,end-start);
 369     fprintf(debug,"q2sum = %g, Vself=%g\n",
 370             L1*q2sumA+lambda*q2sumB,L1*Vself[0]+lambda*Vself[1]);
 371     fprintf(debug,"Long Range correction: Vexcl=%g\n",Vexcl);
 372     if (MASTER(cr) && thread == 0) {
 373       if (epsilon_surface > 0 || ewald_geometry == eewg3DC) {
 374         fprintf(debug,"Total dipole correction: Vdipole=%g\n",
 375                 L1*Vdipole[0]+lambda*Vdipole[1]);
 376       }
 377     }
 378   }
 379
 380   /* Return the correction to the energy */
 381   return enercorr;
 382 }
 383
 384 real ewald_charge_correction(t_commrec *cr,t_forcerec *fr,real lambda,
 385                              matrix box,
 386                              real *dvdlambda,tensor vir)
 387
 388 {
 389     real vol,fac,qs2A,qs2B,vc,enercorr;
 390     int  d;
 391
 392     if (MASTER(cr))
 393     {
 394         /* Apply charge correction */
 395         vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
 396
 397         fac = M_PI*ONE_4PI_EPS0/(fr->epsilon_r*2.0*vol*vol*sqr(fr->ewaldcoeff));
 398
 399         qs2A = fr->qsum[0]*fr->qsum[0];
 400         qs2B = fr->qsum[1]*fr->qsum[1];
 401
 402         vc = (qs2A*(1 - lambda) + qs2B*lambda)*fac;
 403
 404         enercorr = -vol*vc;
 405
 406         *dvdlambda += -vol*(qs2B - qs2A)*fac;
 407
 408         for(d=0; d<DIM; d++)
 409         {
 410             vir[d][d] += vc;
 411         }
 412
 413         if (debug)
 414         {
 415             fprintf(debug,"Total charge correction: Vcharge=%g\n",enercorr);
 416         }
 417     }
 418     else
 419     {
 420         enercorr = 0;
 421     }
 422
 423     return enercorr;
 424 }