src/gmxlib/ewald_util.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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,2013, 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <stdio.h>
  43 #include <math.h>
  44 #include "maths.h"
  45 #include "typedefs.h"
  46 #include "vec.h"
  47 #include "coulomb.h"
  48 #include "smalloc.h"
  49 #include "physics.h"
  50 #include "txtdump.h"
  51 #include "futil.h"
  52 #include "names.h"
  53 #include "writeps.h"
  54 #include "macros.h"
  55
  56 real calc_ewaldcoeff(real rc,real dtol)
  57 {
  58   real x=5,low,high;
  59   int n,i=0;
  60
  61
  62   do {
  63     i++;
  64     x*=2;
  65   } while (gmx_erfc(x*rc) > dtol);
  66
  67   n=i+60; /* search tolerance is 2^-60 */
  68   low=0;
  69   high=x;
  70   for(i=0;i<n;i++) {
  71     x=(low+high)/2;
  72     if (gmx_erfc(x*rc) > dtol)
  73       low=x;
  74     else
  75       high=x;
  76   }
  77   return x;
  78 }
  79
  80
  81
  82 real ewald_LRcorrection(FILE *fplog,
  83                         int start,int end,
  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,
  90                         rvec *f,tensor vir,
  91                         real lambda,real *dvdlambda)
  92 {
  93   int     i,i1,i2,j,k,m,iv,jv,q;
  94   atom_id *AA;
  95   double  q2sumA,q2sumB,Vexcl,dvdl_excl; /* Necessary for precision */
  96   real    one_4pi_eps;
  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;
 100   tensor  dxdf;
 101   real    vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
 102   real    L1,dipole_coeff,qqA,qqB,qqL,vr0;
 103   /*#define TABLES*/
 104 #ifdef TABLES
 105   real    tabscale=fr->tabscale;
 106   real    eps,eps2,VV,FF,F,Y,Geps,Heps2,Fp,fijC,r1t;
 107   real    *VFtab=fr->coulvdwtab;
 108   int     n0,n1,nnn;
 109 #endif
 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*M_2_SQRTPI;
 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   L1         = 1.0-lambda;
 124
 125   /* Note that we have to transform back to gromacs units, since
 126    * mu_tot contains the dipole in debye units (for output).
 127    */
 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;
 131     dipcorrA[i] = 0;
 132     dipcorrB[i] = 0;
 133   }
 134   dipole_coeff=0;
 135   switch (ewald_geometry) {
 136   case eewg3D:
 137     if (epsilon_surface != 0) {
 138       dipole_coeff =
 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];
 143       }
 144     }
 145     break;
 146   case eewg3DC:
 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];
 150     break;
 151   default:
 152     gmx_incons("Unsupported Ewald geometry");
 153     break;
 154   }
 155   if (debug) {
 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]);
 160   }
 161
 162   clear_mat(dxdf);
 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;
 167
 168       if (calc_excl_corr)
 169       {
 170       i1  = excl->index[i];
 171       i2  = excl->index[i+1];
 172
 173       /* Loop over excluded neighbours */
 174       for(j=i1; (j<i2); j++) {
 175         k = AA[j];
 176         /*
 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
 181          * are non-zero.
 182          */
 183         if (k > i) {
 184           qqA = qiA*chargeA[k];
 185           if (qqA != 0.0) {
 186             rvec_sub(x[i],x[k],dx);
 187             if (bMolPBC) {
 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])
 191                   rvec_dec(dx,box[m]);
 192                 else if (dx[m] < -0.5*box[m][m])
 193                   rvec_inc(dx,box[m]);
 194               }
 195             }
 196             dr2 = norm2(dx);
 197             /* Distance between two excluded particles may be zero in the
 198              * case of shells
 199              */
 200             if (dr2 != 0) {
 201               rinv              = gmx_invsqrt(dr2);
 202               rinv2             = rinv*rinv;
 203               dr                = 1.0/rinv;
 204 #ifdef TABLES
 205               r1t               = tabscale*dr;
 206               n0                = r1t;
 207               assert(n0 >= 3);
 208               n1                = 12*n0;
 209               eps               = r1t-n0;
 210               eps2              = eps*eps;
 211               nnn               = n1;
 212               Y                 = VFtab[nnn];
 213               F                 = VFtab[nnn+1];
 214               Geps              = eps*VFtab[nnn+2];
 215               Heps2             = eps2*VFtab[nnn+3];
 216               Fp                = F+Geps+Heps2;
 217               VV                = Y+eps*Fp;
 218               FF                = Fp+Geps+2.0*Heps2;
 219               vc                = qqA*(rinv-VV);
 220               fijC              = qqA*FF;
 221               Vexcl            += vc;
 222
 223               fscal             = vc*rinv2+fijC*tabscale*rinv;
 224             /* End of tabulated interaction part */
 225 #else
 226
 227             /* This is the code you would want instead if not using
 228              * tables:
 229              */
 230               ewcdr   = ewc*dr;
 231               vc      = qqA*gmx_erf(ewcdr)*rinv;
 232               Vexcl  += vc;
 233 #ifdef GMX_DOUBLE
 234               /* Relative accuracy at R_ERF_R_INACC of 3e-10 */
 235 #define       R_ERF_R_INACC 0.006
 236 #else
 237               /* Relative accuracy at R_ERF_R_INACC of 2e-5 */
 238 #define       R_ERF_R_INACC 0.1
 239 #endif
 240               if (ewcdr > R_ERF_R_INACC) {
 241                 fscal = rinv2*(vc - qqA*ewc*M_2_SQRTPI*exp(-ewcdr*ewcdr));
 242               } else {
 243                 /* Use a fourth order series expansion for small ewcdr */
 244                 fscal = ewc*ewc*qqA*vr0*(2.0/3.0 - 0.4*ewcdr*ewcdr);
 245               }
 246 #endif
 247               /* The force vector is obtained by multiplication with the
 248                * distance vector
 249                */
 250               svmul(fscal,dx,df);
 251               rvec_inc(f[k],df);
 252               rvec_dec(f[i],df);
 253               for(iv=0; (iv<DIM); iv++)
 254                 for(jv=0; (jv<DIM); jv++)
 255                   dxdf[iv][jv] += dx[iv]*df[jv];
 256             } else {
 257               Vexcl += qqA*vr0;
 258             }
 259           }
 260         }
 261       }
 262       }
 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];
 267       }
 268     }
 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;
 274
 275       if (calc_excl_corr)
 276       {
 277       i1  = excl->index[i];
 278       i2  = excl->index[i+1];
 279
 280       /* Loop over excluded neighbours */
 281       for(j=i1; (j<i2); j++) {
 282         k = AA[j];
 283         if (k > i) {
 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);
 289             if (bMolPBC) {
 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])
 293                   rvec_dec(dx,box[m]);
 294                 else if (dx[m] < -0.5*box[m][m])
 295                   rvec_inc(dx,box[m]);
 296               }
 297             }
 298             dr2 = norm2(dx);
 299             if (dr2 != 0) {
 300               rinv   = gmx_invsqrt(dr2);
 301               rinv2  = rinv*rinv;
 302               dr     = 1.0/rinv;
 303               v      = gmx_erf(ewc*dr)*rinv;
 304               vc     = qqL*v;
 305               Vexcl += vc;
 306               fscal  = rinv2*(vc-qqL*ewc*M_2_SQRTPI*exp(-ewc*ewc*dr2));
 307               svmul(fscal,dx,df);
 308               rvec_inc(f[k],df);
 309               rvec_dec(f[i],df);
 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;
 314             } else {
 315               Vexcl     +=         qqL*vr0;
 316               dvdl_excl += (qqB - qqA)*vr0;
 317             }
 318           }
 319         }
 320       }
 321       }
 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];
 327       }
 328     }
 329   }
 330   for(iv=0; (iv<DIM); iv++)
 331     for(jv=0; (jv<DIM); jv++)
 332       vir[iv][jv] += 0.5*dxdf[iv][jv];
 333
 334
 335   Vself[0] = 0;
 336   Vself[1] = 0;
 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;
 343       }
 344
 345       /* Apply surface dipole correction:
 346        * correction = dipole_coeff * (dipole)^2
 347        */
 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];
 353       }
 354     }
 355   }
 356
 357   if (!bFreeEnergy) {
 358     enercorr = Vdipole[0] - Vself[0] - Vexcl;
 359    } else {
 360     enercorr = L1*(Vdipole[0] - Vself[0])
 361       + lambda*(Vdipole[1] - Vself[1])
 362       - Vexcl;
 363     *dvdlambda += Vdipole[1] - Vself[1]
 364       - (Vdipole[0] - Vself[0]) - dvdl_excl;
 365   }
 366
 367   if (debug) {
 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]);
 377       }
 378     }
 379   }
 380
 381   /* Return the correction to the energy */
 382   return enercorr;
 383 }
 384
 385 real ewald_charge_correction(t_commrec *cr,t_forcerec *fr,real lambda,
 386                              matrix box,
 387                              real *dvdlambda,tensor vir)
 388
 389 {
 390     real vol,fac,qs2A,qs2B,vc,enercorr;
 391     int  d;
 392
 393     if (MASTER(cr))
 394     {
 395         /* Apply charge correction */
 396         vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
 397
 398         fac = M_PI*ONE_4PI_EPS0/(fr->epsilon_r*2.0*vol*vol*sqr(fr->ewaldcoeff));
 399
 400         qs2A = fr->qsum[0]*fr->qsum[0];
 401         qs2B = fr->qsum[1]*fr->qsum[1];
 402
 403         vc = (qs2A*(1 - lambda) + qs2B*lambda)*fac;
 404
 405         enercorr = -vol*vc;
 406
 407         *dvdlambda += -vol*(qs2B - qs2A)*fac;
 408
 409         for(d=0; d<DIM; d++)
 410         {
 411             vir[d][d] += vc;
 412         }
 413
 414         if (debug)
 415         {
 416             fprintf(debug,"Total charge correction: Vcharge=%g\n",enercorr);
 417         }
 418     }
 419     else
 420     {
 421         enercorr = 0;
 422     }
 423
 424     return enercorr;
 425 }