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