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