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