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