src/gromacs/listed-forces/pairs.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  *
  37  * \brief This file defines functions for "pair" interactions
  38  * (i.e. listed non-bonded interactions, e.g. 1-4 interactions)
  39  *
  40  * \author Mark Abraham <mark.j.abraham@gmail.com>
  41  *
  42  * \ingroup module_listed-forces
  43  */
  44 #include "gmxpre.h"
  45
  46 #include "pairs.h"
  47
  48 #include <cmath>
  49
  50 #include "gromacs/legacyheaders/types/group.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/pbcutil/ishift.h"
  53 #include "gromacs/pbcutil/mshift.h"
  54 #include "gromacs/pbcutil/pbc.h"
  55 #include "gromacs/utility/basedefinitions.h"
  56 #include "gromacs/utility/fatalerror.h"
  57
  58 #include "listed-internal.h"
  59
  60 namespace
  61 {
  62
  63 /*! \brief Issue a warning if a listed interaction is beyond a table limit */
  64 void
  65 warning_rlimit(const rvec *x, int ai, int aj, int * global_atom_index, real r, real rlimit)
  66 {
  67     gmx_warning("Listed nonbonded interaction between particles %d and %d\n"
  68                 "at distance %.3f which is larger than the table limit %.3f nm.\n\n"
  69                 "This is likely either a 1,4 interaction, or a listed interaction inside\n"
  70                 "a smaller molecule you are decoupling during a free energy calculation.\n"
  71                 "Since interactions at distances beyond the table cannot be computed,\n"
  72                 "they are skipped until they are inside the table limit again. You will\n"
  73                 "only see this message once, even if it occurs for several interactions.\n\n"
  74                 "IMPORTANT: This should not happen in a stable simulation, so there is\n"
  75                 "probably something wrong with your system. Only change the table-extension\n"
  76                 "distance in the mdp file if you are really sure that is the reason.\n",
  77                 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r, rlimit);
  78
  79     if (debug)
  80     {
  81         fprintf(debug,
  82                 "%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. Ignored\n",
  83                 x[ai][XX], x[ai][YY], x[ai][ZZ], x[aj][XX], x[aj][YY], x[aj][ZZ],
  84                 glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r);
  85     }
  86 }
  87
  88 /*! \brief Compute the energy and force for a single pair interaction */
  89 real
  90 evaluate_single(real r2, real tabscale, real *vftab,
  91                 real qq, real c6, real c12, real *velec, real *vvdw)
  92 {
  93     real       rinv, r, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VVe, FFe, VVd, FFd, VVr, FFr, fscal;
  94     int        ntab;
  95
  96     /* Do the tabulated interactions - first table lookup */
  97     rinv             = gmx_invsqrt(r2);
  98     r                = r2*rinv;
  99     rtab             = r*tabscale;
 100     ntab             = static_cast<int>(rtab);
 101     eps              = rtab-ntab;
 102     eps2             = eps*eps;
 103     ntab             = 12*ntab;
 104     /* Electrostatics */
 105     Y                = vftab[ntab];
 106     F                = vftab[ntab+1];
 107     Geps             = eps*vftab[ntab+2];
 108     Heps2            = eps2*vftab[ntab+3];
 109     Fp               = F+Geps+Heps2;
 110     VVe              = Y+eps*Fp;
 111     FFe              = Fp+Geps+2.0*Heps2;
 112     /* Dispersion */
 113     Y                = vftab[ntab+4];
 114     F                = vftab[ntab+5];
 115     Geps             = eps*vftab[ntab+6];
 116     Heps2            = eps2*vftab[ntab+7];
 117     Fp               = F+Geps+Heps2;
 118     VVd              = Y+eps*Fp;
 119     FFd              = Fp+Geps+2.0*Heps2;
 120     /* Repulsion */
 121     Y                = vftab[ntab+8];
 122     F                = vftab[ntab+9];
 123     Geps             = eps*vftab[ntab+10];
 124     Heps2            = eps2*vftab[ntab+11];
 125     Fp               = F+Geps+Heps2;
 126     VVr              = Y+eps*Fp;
 127     FFr              = Fp+Geps+2.0*Heps2;
 128
 129     *velec           = qq*VVe;
 130     *vvdw            = c6*VVd+c12*VVr;
 131
 132     fscal            = -(qq*FFe+c6*FFd+c12*FFr)*tabscale*rinv;
 133
 134     return fscal;
 135 }
 136
 137 /*! \brief Compute the energy and force for a single pair interaction under FEP */
 138 real
 139 free_energy_evaluate_single(real r2, real sc_r_power, real alpha_coul,
 140                             real alpha_vdw, real tabscale, real *vftab,
 141                             real qqA, real c6A, real c12A, real qqB,
 142                             real c6B, real c12B, real LFC[2], real LFV[2], real DLF[2],
 143                             real lfac_coul[2], real lfac_vdw[2], real dlfac_coul[2],
 144                             real dlfac_vdw[2], real sigma6_def, real sigma6_min,
 145                             real sigma2_def, real sigma2_min,
 146                             real *velectot, real *vvdwtot, real *dvdl)
 147 {
 148     real       rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal;
 149     real       qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2];
 150     real       alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw;
 151     real       rpinv, r_coul, r_vdw, velecsum, vvdwsum;
 152     real       fscal_vdw[2], fscal_elec[2];
 153     real       velec[2], vvdw[2];
 154     int        i, ntab;
 155     const real half        = 0.5;
 156     const real minusOne    = -1.0;
 157     const real oneThird    = 1.0 / 3.0;
 158     const real one         = 1.0;
 159     const real two         = 2.0;
 160     const real six         = 6.0;
 161     const real fourtyeight = 48.0;
 162
 163     qq[0]    = qqA;
 164     qq[1]    = qqB;
 165     c6[0]    = c6A;
 166     c6[1]    = c6B;
 167     c12[0]   = c12A;
 168     c12[1]   = c12B;
 169
 170     if (sc_r_power == six)
 171     {
 172         rpm2             = r2*r2;   /* r4 */
 173         rp               = rpm2*r2; /* r6 */
 174     }
 175     else if (sc_r_power == fourtyeight)
 176     {
 177         rp               = r2*r2*r2; /* r6 */
 178         rp               = rp*rp;    /* r12 */
 179         rp               = rp*rp;    /* r24 */
 180         rp               = rp*rp;    /* r48 */
 181         rpm2             = rp/r2;    /* r46 */
 182     }
 183     else
 184     {
 185         rp             = std::pow(r2, half * sc_r_power);  /* not currently supported as input, but can handle it */
 186         rpm2           = rp/r2;
 187     }
 188
 189     /* Loop over state A(0) and B(1) */
 190     for (i = 0; i < 2; i++)
 191     {
 192         if ((c6[i] > 0) && (c12[i] > 0))
 193         {
 194             /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively.
 195              * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5.
 196              */
 197             sigma6[i]       = half*c12[i]/c6[i];
 198             sigma2[i]       = std::pow(half*c12[i]/c6[i], oneThird);
 199             /* should be able to get rid of this ^^^ internal pow call eventually.  Will require agreement on
 200                what data to store externally.  Can't be fixed without larger scale changes, so not 5.0 */
 201             if (sigma6[i] < sigma6_min)   /* for disappearing coul and vdw with soft core at the same time */
 202             {
 203                 sigma6[i] = sigma6_min;
 204                 sigma2[i] = sigma2_min;
 205             }
 206         }
 207         else
 208         {
 209             sigma6[i]       = sigma6_def;
 210             sigma2[i]       = sigma2_def;
 211         }
 212         if (sc_r_power == six)
 213         {
 214             sigma_pow[i]    = sigma6[i];
 215         }
 216         else if (sc_r_power == fourtyeight)
 217         {
 218             sigma_pow[i]    = sigma6[i]*sigma6[i];       /* sigma^12 */
 219             sigma_pow[i]    = sigma_pow[i]*sigma_pow[i]; /* sigma^24 */
 220             sigma_pow[i]    = sigma_pow[i]*sigma_pow[i]; /* sigma^48 */
 221         }
 222         else
 223         {       /* not really supported as input, but in here for testing the general case*/
 224             sigma_pow[i]    = std::pow(sigma2[i], sc_r_power/2);
 225         }
 226     }
 227
 228     /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/
 229     if ((c12[0] > 0) && (c12[1] > 0))
 230     {
 231         alpha_vdw_eff    = 0;
 232         alpha_coul_eff   = 0;
 233     }
 234     else
 235     {
 236         alpha_vdw_eff    = alpha_vdw;
 237         alpha_coul_eff   = alpha_coul;
 238     }
 239
 240     /* Loop over A and B states again */
 241     for (i = 0; i < 2; i++)
 242     {
 243         fscal_elec[i] = 0;
 244         fscal_vdw[i]  = 0;
 245         velec[i]      = 0;
 246         vvdw[i]       = 0;
 247
 248         /* Only spend time on A or B state if it is non-zero */
 249         if ( (qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0) )
 250         {
 251             /* Coulomb */
 252             rpinv            = one/(alpha_coul_eff*lfac_coul[i]*sigma_pow[i]+rp);
 253             r_coul           = std::pow(rpinv, minusOne / sc_r_power);
 254
 255             /* Electrostatics table lookup data */
 256             rtab             = r_coul*tabscale;
 257             ntab             = static_cast<int>(rtab);
 258             eps              = rtab-ntab;
 259             eps2             = eps*eps;
 260             ntab             = 12*ntab;
 261             /* Electrostatics */
 262             Y                = vftab[ntab];
 263             F                = vftab[ntab+1];
 264             Geps             = eps*vftab[ntab+2];
 265             Heps2            = eps2*vftab[ntab+3];
 266             Fp               = F+Geps+Heps2;
 267             VV               = Y+eps*Fp;
 268             FF               = Fp+Geps+two*Heps2;
 269             velec[i]         = qq[i]*VV;
 270             fscal_elec[i]    = -qq[i]*FF*r_coul*rpinv*tabscale;
 271
 272             /* Vdw */
 273             rpinv            = one/(alpha_vdw_eff*lfac_vdw[i]*sigma_pow[i]+rp);
 274             r_vdw            = std::pow(rpinv, minusOne / sc_r_power);
 275             /* Vdw table lookup data */
 276             rtab             = r_vdw*tabscale;
 277             ntab             = static_cast<int>(rtab);
 278             eps              = rtab-ntab;
 279             eps2             = eps*eps;
 280             ntab             = 12*ntab;
 281             /* Dispersion */
 282             Y                = vftab[ntab+4];
 283             F                = vftab[ntab+5];
 284             Geps             = eps*vftab[ntab+6];
 285             Heps2            = eps2*vftab[ntab+7];
 286             Fp               = F+Geps+Heps2;
 287             VV               = Y+eps*Fp;
 288             FF               = Fp+Geps+two*Heps2;
 289             vvdw[i]          = c6[i]*VV;
 290             fscal_vdw[i]     = -c6[i]*FF;
 291
 292             /* Repulsion */
 293             Y                = vftab[ntab+8];
 294             F                = vftab[ntab+9];
 295             Geps             = eps*vftab[ntab+10];
 296             Heps2            = eps2*vftab[ntab+11];
 297             Fp               = F+Geps+Heps2;
 298             VV               = Y+eps*Fp;
 299             FF               = Fp+Geps+two*Heps2;
 300             vvdw[i]         += c12[i]*VV;
 301             fscal_vdw[i]    -= c12[i]*FF;
 302             fscal_vdw[i]    *= r_vdw*rpinv*tabscale;
 303         }
 304     }
 305     /* Now we have velec[i], vvdw[i], and fscal[i] for both states */
 306     /* Assemble A and B states */
 307     velecsum  = 0;
 308     vvdwsum   = 0;
 309     dvdl_coul = 0;
 310     dvdl_vdw  = 0;
 311     fscal     = 0;
 312     for (i = 0; i < 2; i++)
 313     {
 314         velecsum      += LFC[i]*velec[i];
 315         vvdwsum       += LFV[i]*vvdw[i];
 316
 317         fscal         += (LFC[i]*fscal_elec[i]+LFV[i]*fscal_vdw[i])*rpm2;
 318
 319         dvdl_coul     += velec[i]*DLF[i] + LFC[i]*alpha_coul_eff*dlfac_coul[i]*fscal_elec[i]*sigma_pow[i];
 320         dvdl_vdw      += vvdw[i]*DLF[i] + LFV[i]*alpha_vdw_eff*dlfac_vdw[i]*fscal_vdw[i]*sigma_pow[i];
 321     }
 322
 323     dvdl[efptCOUL]     += dvdl_coul;
 324     dvdl[efptVDW]      += dvdl_vdw;
 325
 326     *velectot           = velecsum;
 327     *vvdwtot            = vvdwsum;
 328
 329     return fscal;
 330 }
 331
 332 } // namespace
 333
 334 real
 335 do_pairs(int ftype, int nbonds,
 336          const t_iatom iatoms[], const t_iparams iparams[],
 337          const rvec x[], rvec f[], rvec fshift[],
 338          const struct t_pbc *pbc, const struct t_graph *g,
 339          real *lambda, real *dvdl,
 340          const t_mdatoms *md,
 341          const t_forcerec *fr, gmx_grppairener_t *grppener,
 342          int *global_atom_index)
 343 {
 344     real             qq, c6, c12;
 345     rvec             dx;
 346     ivec             dt;
 347     int              i, itype, ai, aj, gid;
 348     int              fshift_index;
 349     real             r2;
 350     real             fscal, velec, vvdw;
 351     real *           energygrp_elec;
 352     real *           energygrp_vdw;
 353     static gmx_bool  warned_rlimit = FALSE;
 354     /* Free energy stuff */
 355     gmx_bool         bFreeEnergy;
 356     real             LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
 357     real             qqB, c6B, c12B, sigma2_def, sigma2_min;
 358     const real       oneThird = 1.0 / 3.0;
 359
 360     switch (ftype)
 361     {
 362         case F_LJ14:
 363         case F_LJC14_Q:
 364             energygrp_elec = grppener->ener[egCOUL14];
 365             energygrp_vdw  = grppener->ener[egLJ14];
 366             break;
 367         case F_LJC_PAIRS_NB:
 368             energygrp_elec = grppener->ener[egCOULSR];
 369             energygrp_vdw  = grppener->ener[egLJSR];
 370             break;
 371         default:
 372             energygrp_elec = NULL; /* Keep compiler happy */
 373             energygrp_vdw  = NULL; /* Keep compiler happy */
 374             gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", ftype);
 375             break;
 376     }
 377
 378     if (fr->efep != efepNO)
 379     {
 380         /* Lambda factor for state A=1-lambda and B=lambda */
 381         LFC[0] = 1.0 - lambda[efptCOUL];
 382         LFV[0] = 1.0 - lambda[efptVDW];
 383         LFC[1] = lambda[efptCOUL];
 384         LFV[1] = lambda[efptVDW];
 385
 386         /*derivative of the lambda factor for state A and B */
 387         DLF[0] = -1;
 388         DLF[1] = 1;
 389
 390         /* precalculate */
 391         sigma2_def = pow(fr->sc_sigma6_def, oneThird);
 392         sigma2_min = pow(fr->sc_sigma6_min, oneThird);
 393
 394         for (i = 0; i < 2; i++)
 395         {
 396             lfac_coul[i]  = (fr->sc_power == 2 ? (1-LFC[i])*(1-LFC[i]) : (1-LFC[i]));
 397             dlfac_coul[i] = DLF[i]*fr->sc_power/fr->sc_r_power*(fr->sc_power == 2 ? (1-LFC[i]) : 1);
 398             lfac_vdw[i]   = (fr->sc_power == 2 ? (1-LFV[i])*(1-LFV[i]) : (1-LFV[i]));
 399             dlfac_vdw[i]  = DLF[i]*fr->sc_power/fr->sc_r_power*(fr->sc_power == 2 ? (1-LFV[i]) : 1);
 400         }
 401     }
 402     else
 403     {
 404         sigma2_min = sigma2_def = 0;
 405     }
 406
 407     bFreeEnergy = FALSE;
 408     for (i = 0; (i < nbonds); )
 409     {
 410         itype = iatoms[i++];
 411         ai    = iatoms[i++];
 412         aj    = iatoms[i++];
 413         gid   = GID(md->cENER[ai], md->cENER[aj], md->nenergrp);
 414
 415         /* Get parameters */
 416         switch (ftype)
 417         {
 418             case F_LJ14:
 419                 bFreeEnergy =
 420                     (fr->efep != efepNO &&
 421                      ((md->nPerturbed && (md->bPerturbed[ai] || md->bPerturbed[aj])) ||
 422                       iparams[itype].lj14.c6A != iparams[itype].lj14.c6B ||
 423                       iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
 424                 qq               = md->chargeA[ai]*md->chargeA[aj]*fr->epsfac*fr->fudgeQQ;
 425                 c6               = iparams[itype].lj14.c6A;
 426                 c12              = iparams[itype].lj14.c12A;
 427                 break;
 428             case F_LJC14_Q:
 429                 qq               = iparams[itype].ljc14.qi*iparams[itype].ljc14.qj*fr->epsfac*iparams[itype].ljc14.fqq;
 430                 c6               = iparams[itype].ljc14.c6;
 431                 c12              = iparams[itype].ljc14.c12;
 432                 break;
 433             case F_LJC_PAIRS_NB:
 434                 qq               = iparams[itype].ljcnb.qi*iparams[itype].ljcnb.qj*fr->epsfac;
 435                 c6               = iparams[itype].ljcnb.c6;
 436                 c12              = iparams[itype].ljcnb.c12;
 437                 break;
 438             default:
 439                 /* Cannot happen since we called gmx_fatal() above in this case */
 440                 qq = c6 = c12 = 0; /* Keep compiler happy */
 441                 break;
 442         }
 443
 444         /* To save flops in the optimized kernels, c6/c12 have 6.0/12.0 derivative prefactors
 445          * included in the general nfbp array now. This means the tables are scaled down by the
 446          * same factor, so when we use the original c6/c12 parameters from iparams[] they must
 447          * be scaled up.
 448          */
 449         c6  *= 6.0;
 450         c12 *= 12.0;
 451
 452         /* Do we need to apply full periodic boundary conditions? */
 453         if (fr->bMolPBC == TRUE)
 454         {
 455             fshift_index = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
 456         }
 457         else
 458         {
 459             fshift_index = CENTRAL;
 460             rvec_sub(x[ai], x[aj], dx);
 461         }
 462         r2           = norm2(dx);
 463
 464         if (r2 >= fr->tab14.r*fr->tab14.r)
 465         {
 466             /* This check isn't race free. But it doesn't matter because if a race occurs the only
 467              * disadvantage is that the warning is printed twice */
 468             if (warned_rlimit == FALSE)
 469             {
 470                 warning_rlimit(x, ai, aj, global_atom_index, sqrt(r2), fr->tab14.r);
 471                 warned_rlimit = TRUE;
 472             }
 473             continue;
 474         }
 475
 476         if (bFreeEnergy)
 477         {
 478             /* Currently free energy is only supported for F_LJ14, so no need to check for that if we got here */
 479             qqB              = md->chargeB[ai]*md->chargeB[aj]*fr->epsfac*fr->fudgeQQ;
 480             c6B              = iparams[itype].lj14.c6B*6.0;
 481             c12B             = iparams[itype].lj14.c12B*12.0;
 482
 483             fscal            = free_energy_evaluate_single(r2, fr->sc_r_power, fr->sc_alphacoul, fr->sc_alphavdw,
 484                                                            fr->tab14.scale, fr->tab14.data, qq, c6, c12, qqB, c6B, c12B,
 485                                                            LFC, LFV, DLF, lfac_coul, lfac_vdw, dlfac_coul, dlfac_vdw,
 486                                                            fr->sc_sigma6_def, fr->sc_sigma6_min, sigma2_def, sigma2_min, &velec, &vvdw, dvdl);
 487         }
 488         else
 489         {
 490             /* Evaluate tabulated interaction without free energy */
 491             fscal            = evaluate_single(r2, fr->tab14.scale, fr->tab14.data, qq, c6, c12, &velec, &vvdw);
 492         }
 493
 494         energygrp_elec[gid]  += velec;
 495         energygrp_vdw[gid]   += vvdw;
 496         svmul(fscal, dx, dx);
 497
 498         /* Add the forces */
 499         rvec_inc(f[ai], dx);
 500         rvec_dec(f[aj], dx);
 501
 502         if (g)
 503         {
 504             /* Correct the shift forces using the graph */
 505             ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
 506             fshift_index = IVEC2IS(dt);
 507         }
 508         if (fshift_index != CENTRAL)
 509         {
 510             rvec_inc(fshift[fshift_index], dx);
 511             rvec_dec(fshift[CENTRAL], dx);
 512         }
 513     }
 514     return 0.0;
 515 }