src/gromacs/mdlib/dispersioncorrection.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020,2021, 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 #include "gmxpre.h"
  36
  37 #include "dispersioncorrection.h"
  38
  39 #include <cstdio>
  40
  41 #include "gromacs/math/units.h"
  42 #include "gromacs/math/utilities.h"
  43 #include "gromacs/math/vec.h"
  44 #include "gromacs/mdtypes/forcerec.h"
  45 #include "gromacs/mdtypes/inputrec.h"
  46 #include "gromacs/mdtypes/interaction_const.h"
  47 #include "gromacs/mdtypes/md_enums.h"
  48 #include "gromacs/tables/forcetable.h"
  49 #include "gromacs/topology/mtop_util.h"
  50 #include "gromacs/topology/topology.h"
  51 #include "gromacs/utility/arrayref.h"
  52 #include "gromacs/utility/fatalerror.h"
  53 #include "gromacs/utility/logger.h"
  54
  55 /* Implementation here to avoid other files needing to include the file that defines t_nblists */
  56 DispersionCorrection::InteractionParams::~InteractionParams() = default;
  57
  58 /* Returns a matrix, as flat list, of combination rule combined LJ parameters */
  59 static std::vector<real> mk_nbfp_combination_rule(const gmx_ffparams_t& ffparams,
  60                                                   const CombinationRule comb_rule)
  61 {
  62     const int atnr = ffparams.atnr;
  63
  64     std::vector<real> nbfp(atnr * atnr * 2);
  65
  66     for (int i = 0; i < atnr; ++i)
  67     {
  68         for (int j = 0; j < atnr; ++j)
  69         {
  70             const real c6i  = ffparams.iparams[i * (atnr + 1)].lj.c6;
  71             const real c12i = ffparams.iparams[i * (atnr + 1)].lj.c12;
  72             const real c6j  = ffparams.iparams[j * (atnr + 1)].lj.c6;
  73             const real c12j = ffparams.iparams[j * (atnr + 1)].lj.c12;
  74             real       c6   = std::sqrt(c6i * c6j);
  75             real       c12  = std::sqrt(c12i * c12j);
  76             if (comb_rule == CombinationRule::Arithmetic && !gmx_numzero(c6) && !gmx_numzero(c12))
  77             {
  78                 const real sigmai = gmx::sixthroot(c12i / c6i);
  79                 const real sigmaj = gmx::sixthroot(c12j / c6j);
  80                 const real epsi   = c6i * c6i / c12i;
  81                 const real epsj   = c6j * c6j / c12j;
  82                 const real sigma  = 0.5 * (sigmai + sigmaj);
  83                 const real eps    = std::sqrt(epsi * epsj);
  84                 c6                = eps * gmx::power6(sigma);
  85                 c12               = eps * gmx::power12(sigma);
  86             }
  87             C6(nbfp, atnr, i, j)  = c6 * 6.0;
  88             C12(nbfp, atnr, i, j) = c12 * 12.0;
  89         }
  90     }
  91
  92     return nbfp;
  93 }
  94
  95 /* Returns the A-topology atom type when aOrB=0, the B-topology atom type when aOrB=1 */
  96 static int atomtypeAOrB(const t_atom& atom, int aOrB)
  97 {
  98     if (aOrB == 0)
  99     {
 100         return atom.type;
 101     }
 102     else
 103     {
 104         return atom.typeB;
 105     }
 106 }
 107
 108 DispersionCorrection::TopologyParams::TopologyParams(const gmx_mtop_t&         mtop,
 109                                                      const t_inputrec&         inputrec,
 110                                                      bool                      useBuckingham,
 111                                                      int                       numAtomTypes,
 112                                                      gmx::ArrayRef<const real> nonbondedForceParameters)
 113 {
 114     const int      ntp   = numAtomTypes;
 115     const gmx_bool bBHAM = useBuckingham;
 116
 117     gmx::ArrayRef<const real> nbfp = nonbondedForceParameters;
 118     std::vector<real>         nbfp_comb;
 119     /* For LJ-PME, we want to correct for the difference between the
 120      * actual C6 values and the C6 values used by the LJ-PME based on
 121      * combination rules. */
 122     if (EVDW_PME(inputrec.vdwtype))
 123     {
 124         nbfp_comb = mk_nbfp_combination_rule(mtop.ffparams,
 125                                              (inputrec.ljpme_combination_rule == LongRangeVdW::LB)
 126                                                      ? CombinationRule::Arithmetic
 127                                                      : CombinationRule::Geometric);
 128         for (int tpi = 0; tpi < ntp; ++tpi)
 129         {
 130             for (int tpj = 0; tpj < ntp; ++tpj)
 131             {
 132                 C6(nbfp_comb, ntp, tpi, tpj) = C6(nbfp, ntp, tpi, tpj) - C6(nbfp_comb, ntp, tpi, tpj);
 133                 C12(nbfp_comb, ntp, tpi, tpj) = C12(nbfp, ntp, tpi, tpj);
 134             }
 135         }
 136         nbfp = nbfp_comb;
 137     }
 138
 139     for (int q = 0; q < (inputrec.efep == FreeEnergyPerturbationType::No ? 1 : 2); q++)
 140     {
 141         double  csix    = 0;
 142         double  ctwelve = 0;
 143         int64_t npair   = 0;
 144         int64_t nexcl   = 0;
 145         if (!EI_TPI(inputrec.eI))
 146         {
 147             numAtomsForDensity_ = mtop.natoms;
 148             numCorrections_     = 0.5 * mtop.natoms;
 149
 150             /* Count the types so we avoid natoms^2 operations */
 151             std::vector<int> typecount(ntp);
 152             gmx_mtop_count_atomtypes(mtop, q, typecount.data());
 153
 154             for (int tpi = 0; tpi < ntp; tpi++)
 155             {
 156                 for (int tpj = tpi; tpj < ntp; tpj++)
 157                 {
 158                     const int64_t iCount = typecount[tpi];
 159                     const int64_t jCount = typecount[tpj];
 160                     int64_t       npair_ij;
 161                     if (tpi != tpj)
 162                     {
 163                         npair_ij = iCount * jCount;
 164                     }
 165                     else
 166                     {
 167                         npair_ij = iCount * (iCount - 1) / 2;
 168                     }
 169                     if (bBHAM)
 170                     {
 171                         /* nbfp now includes the 6.0 derivative prefactor */
 172                         csix += npair_ij * BHAMC(nbfp, ntp, tpi, tpj) / 6.0;
 173                     }
 174                     else
 175                     {
 176                         /* nbfp now includes the 6.0/12.0 derivative prefactors */
 177                         csix += npair_ij * C6(nbfp, ntp, tpi, tpj) / 6.0;
 178                         ctwelve += npair_ij * C12(nbfp, ntp, tpi, tpj) / 12.0;
 179                     }
 180                     npair += npair_ij;
 181                 }
 182             }
 183             /* Subtract the excluded pairs.
 184              * The main reason for substracting exclusions is that in some cases
 185              * some combinations might never occur and the parameters could have
 186              * any value. These unused values should not influence the dispersion
 187              * correction.
 188              */
 189             for (const gmx_molblock_t& molb : mtop.molblock)
 190             {
 191                 const int      nmol  = molb.nmol;
 192                 const t_atoms* atoms = &mtop.moltype[molb.type].atoms;
 193                 const auto&    excl  = mtop.moltype[molb.type].excls;
 194                 for (int i = 0; (i < atoms->nr); i++)
 195                 {
 196                     const int tpi = atomtypeAOrB(atoms->atom[i], q);
 197                     for (const int k : excl[i])
 198                     {
 199                         if (k > i)
 200                         {
 201                             const int tpj = atomtypeAOrB(atoms->atom[k], q);
 202                             if (bBHAM)
 203                             {
 204                                 /* nbfp now includes the 6.0 derivative prefactor */
 205                                 csix -= nmol * BHAMC(nbfp, ntp, tpi, tpj) / 6.0;
 206                             }
 207                             else
 208                             {
 209                                 /* nbfp now includes the 6.0/12.0 derivative prefactors */
 210                                 csix -= nmol * C6(nbfp, ntp, tpi, tpj) / 6.0;
 211                                 ctwelve -= nmol * C12(nbfp, ntp, tpi, tpj) / 12.0;
 212                             }
 213                             nexcl += molb.nmol;
 214                         }
 215                     }
 216                 }
 217             }
 218         }
 219         else
 220         {
 221             const t_atoms& atoms_tpi = mtop.moltype[mtop.molblock.back().type].atoms;
 222
 223             /* Only correct for the interaction of the test particle
 224              * with the rest of the system.
 225              */
 226             numAtomsForDensity_ = mtop.natoms - atoms_tpi.nr;
 227             numCorrections_     = atoms_tpi.nr;
 228
 229             npair = 0;
 230             for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
 231             {
 232                 const gmx_molblock_t& molb  = mtop.molblock[mb];
 233                 const t_atoms&        atoms = mtop.moltype[molb.type].atoms;
 234                 for (int j = 0; j < atoms.nr; j++)
 235                 {
 236                     int nmolc = molb.nmol;
 237                     /* Remove the interaction of the test charge group
 238                      * with itself.
 239                      */
 240                     if (mb == mtop.molblock.size() - 1)
 241                     {
 242                         nmolc--;
 243
 244                         if (mb == 0 && molb.nmol == 1)
 245                         {
 246                             gmx_fatal(FARGS,
 247                                       "Old format tpr with TPI, please generate a new tpr file");
 248                         }
 249                     }
 250                     const int tpj = atomtypeAOrB(atoms.atom[j], q);
 251                     for (int i = 0; i < atoms_tpi.nr; i++)
 252                     {
 253                         const int tpi = atomtypeAOrB(atoms_tpi.atom[i], q);
 254                         if (bBHAM)
 255                         {
 256                             /* nbfp now includes the 6.0 derivative prefactor */
 257                             csix += nmolc * BHAMC(nbfp, ntp, tpi, tpj) / 6.0;
 258                         }
 259                         else
 260                         {
 261                             /* nbfp now includes the 6.0/12.0 derivative prefactors */
 262                             csix += nmolc * C6(nbfp, ntp, tpi, tpj) / 6.0;
 263                             ctwelve += nmolc * C12(nbfp, ntp, tpi, tpj) / 12.0;
 264                         }
 265                         npair += nmolc;
 266                     }
 267                 }
 268             }
 269         }
 270         if (npair - nexcl <= 0)
 271         {
 272             csix    = 0;
 273             ctwelve = 0;
 274         }
 275         else
 276         {
 277             csix /= npair - nexcl;
 278             ctwelve /= npair - nexcl;
 279         }
 280         if (debug)
 281         {
 282             fprintf(debug, "Counted %" PRId64 " exclusions\n", nexcl);
 283             fprintf(debug, "Average C6 parameter is: %10g\n", csix);
 284             fprintf(debug, "Average C12 parameter is: %10g\n", ctwelve);
 285         }
 286         avcsix_[q]    = csix;
 287         avctwelve_[q] = ctwelve;
 288     }
 289 }
 290
 291 static void integrate_table(const real vdwtab[],
 292                             const real scale,
 293                             const int  offstart,
 294                             const int  rstart,
 295                             const int  rend,
 296                             double*    enerout,
 297                             double*    virout)
 298 {
 299     const double invscale  = 1.0 / scale;
 300     const double invscale2 = invscale * invscale;
 301     const double invscale3 = invscale * invscale2;
 302
 303     /* Following summation derived from cubic spline definition,
 304      * Numerical Recipies in C, second edition, p. 113-116.  Exact for
 305      * the cubic spline.  We first calculate the negative of the
 306      * energy from rvdw to rvdw_switch, assuming that g(r)=1, and then
 307      * add the more standard, abrupt cutoff correction to that result,
 308      * yielding the long-range correction for a switched function.  We
 309      * perform both the pressure and energy loops at the same time for
 310      * simplicity, as the computational cost is low. */
 311
 312     /* Since the dispersion table has been scaled down a factor
 313      * 6.0 and the repulsion a factor 12.0 to compensate for the
 314      * c6/c12 parameters inside nbfp[] being scaled up (to save
 315      * flops in kernels), we need to correct for this.
 316      */
 317     const double tabfactor = (offstart == 0 ? 6.0 : 12.0);
 318
 319     double enersum = 0.0;
 320     double virsum  = 0.0;
 321     for (int ri = rstart; ri < rend; ++ri)
 322     {
 323         const double r  = ri * invscale;
 324         const double ea = invscale3;
 325         const double eb = 2.0 * invscale2 * r;
 326         const double ec = invscale * r * r;
 327
 328         const double pa = invscale3;
 329         const double pb = 3.0 * invscale2 * r;
 330         const double pc = 3.0 * invscale * r * r;
 331         const double pd = r * r * r;
 332
 333         /* this "8" is from the packing in the vdwtab array - perhaps
 334            should be defined? */
 335
 336         const int    offset = 8 * ri + offstart;
 337         const double y0     = vdwtab[offset];
 338         const double f      = vdwtab[offset + 1];
 339         const double g      = vdwtab[offset + 2];
 340         const double h      = vdwtab[offset + 3];
 341
 342         enersum += y0 * (ea / 3 + eb / 2 + ec) + f * (ea / 4 + eb / 3 + ec / 2)
 343                    + g * (ea / 5 + eb / 4 + ec / 3) + h * (ea / 6 + eb / 5 + ec / 4);
 344         virsum += f * (pa / 4 + pb / 3 + pc / 2 + pd) + 2 * g * (pa / 5 + pb / 4 + pc / 3 + pd / 2)
 345                   + 3 * h * (pa / 6 + pb / 5 + pc / 4 + pd / 3);
 346     }
 347     *enerout = 4.0 * M_PI * enersum * tabfactor;
 348     *virout  = 4.0 * M_PI * virsum * tabfactor;
 349 }
 350
 351 /* Struct for storing and passing energy or virial corrections */
 352 struct InteractionCorrection
 353 {
 354     real dispersion = 0;
 355     real repulsion  = 0;
 356 };
 357
 358 /* Adds the energy and virial corrections beyond the cut-off */
 359 static void addCorrectionBeyondCutoff(InteractionCorrection* energy,
 360                                       InteractionCorrection* virial,
 361                                       const double           cutoffDistance)
 362 {
 363     const double rc3 = cutoffDistance * cutoffDistance * cutoffDistance;
 364     const double rc9 = rc3 * rc3 * rc3;
 365
 366     energy->dispersion += -4.0 * M_PI / (3.0 * rc3);
 367     energy->repulsion += 4.0 * M_PI / (9.0 * rc9);
 368     virial->dispersion += 8.0 * M_PI / rc3;
 369     virial->repulsion += -16.0 * M_PI / (3.0 * rc9);
 370 }
 371
 372 void DispersionCorrection::setInteractionParameters(InteractionParams*         iParams,
 373                                                     const interaction_const_t& ic,
 374                                                     const char*                tableFileName)
 375 {
 376     /* We only need to set the tables at first call, i.e. tableFileName!=nullptr
 377      * or when we changed the cut-off with LJ-PME tuning.
 378      */
 379     if (tableFileName || EVDW_PME(ic.vdwtype))
 380     {
 381         iParams->dispersionCorrectionTable_ =
 382                 makeDispersionCorrectionTable(nullptr, &ic, ic.rvdw, tableFileName);
 383     }
 384
 385     InteractionCorrection energy;
 386     InteractionCorrection virial;
 387
 388     if ((ic.vdw_modifier == InteractionModifiers::PotShift)
 389         || (ic.vdw_modifier == InteractionModifiers::PotSwitch)
 390         || (ic.vdw_modifier == InteractionModifiers::ForceSwitch)
 391         || (ic.vdwtype == VanDerWaalsType::Shift) || (ic.vdwtype == VanDerWaalsType::Switch))
 392     {
 393         if (((ic.vdw_modifier == InteractionModifiers::PotSwitch)
 394              || (ic.vdw_modifier == InteractionModifiers::ForceSwitch)
 395              || (ic.vdwtype == VanDerWaalsType::Switch))
 396             && ic.rvdw_switch == 0)
 397         {
 398             gmx_fatal(FARGS,
 399                       "With dispersion correction rvdw-switch can not be zero "
 400                       "for vdw-type = %s",
 401                       enumValueToString(ic.vdwtype));
 402         }
 403
 404         GMX_ASSERT(iParams->dispersionCorrectionTable_, "We need an initialized table");
 405
 406         /* TODO This code depends on the logic in tables.c that
 407            constructs the table layout, which should be made
 408            explicit in future cleanup. */
 409         GMX_ASSERT(iParams->dispersionCorrectionTable_->interaction == TableInteraction::VdwRepulsionVdwDispersion,
 410                    "Dispersion-correction code needs a table with both repulsion and dispersion "
 411                    "terms");
 412         const real  scale  = iParams->dispersionCorrectionTable_->scale;
 413         const real* vdwtab = iParams->dispersionCorrectionTable_->data.data();
 414
 415         /* Round the cut-offs to exact table values for precision */
 416         int ri0 = static_cast<int>(std::floor(ic.rvdw_switch * scale));
 417         int ri1 = static_cast<int>(std::ceil(ic.rvdw * scale));
 418
 419         /* The code below has some support for handling force-switching, i.e.
 420          * when the force (instead of potential) is switched over a limited
 421          * region. This leads to a constant shift in the potential inside the
 422          * switching region, which we can handle by adding a constant energy
 423          * term in the force-switch case just like when we do potential-shift.
 424          *
 425          * For now this is not enabled, but to keep the functionality in the
 426          * code we check separately for switch and shift. When we do force-switch
 427          * the shifting point is rvdw_switch, while it is the cutoff when we
 428          * have a classical potential-shift.
 429          *
 430          * For a pure potential-shift the potential has a constant shift
 431          * all the way out to the cutoff, and that is it. For other forms
 432          * we need to calculate the constant shift up to the point where we
 433          * start modifying the potential.
 434          */
 435         ri0 = (ic.vdw_modifier == InteractionModifiers::PotShift) ? ri1 : ri0;
 436
 437         const double r0  = ri0 / scale;
 438         const double rc3 = r0 * r0 * r0;
 439         const double rc9 = rc3 * rc3 * rc3;
 440
 441         if ((ic.vdw_modifier == InteractionModifiers::ForceSwitch) || (ic.vdwtype == VanDerWaalsType::Shift))
 442         {
 443             /* Determine the constant energy shift below rvdw_switch.
 444              * Table has a scale factor since we have scaled it down to compensate
 445              * for scaling-up c6/c12 with the derivative factors to save flops in analytical kernels.
 446              */
 447             iParams->enershiftsix_ = static_cast<real>(-1.0 / (rc3 * rc3)) - 6.0 * vdwtab[8 * ri0];
 448             iParams->enershifttwelve_ = static_cast<real>(1.0 / (rc9 * rc3)) - 12.0 * vdwtab[8 * ri0 + 4];
 449         }
 450         else if (ic.vdw_modifier == InteractionModifiers::PotShift)
 451         {
 452             iParams->enershiftsix_    = static_cast<real>(-1.0 / (rc3 * rc3));
 453             iParams->enershifttwelve_ = static_cast<real>(1.0 / (rc9 * rc3));
 454         }
 455
 456         /* Add the constant part from 0 to rvdw_switch.
 457          * This integration from 0 to rvdw_switch overcounts the number
 458          * of interactions by 1, as it also counts the self interaction.
 459          * We will correct for this later.
 460          */
 461         energy.dispersion += 4.0 * M_PI * iParams->enershiftsix_ * rc3 / 3.0;
 462         energy.repulsion += 4.0 * M_PI * iParams->enershifttwelve_ * rc3 / 3.0;
 463
 464         /* Calculate the contribution in the range [r0,r1] where we
 465          * modify the potential. For a pure potential-shift modifier we will
 466          * have ri0==ri1, and there will not be any contribution here.
 467          */
 468         double enersum = 0;
 469         double virsum  = 0;
 470         integrate_table(vdwtab, scale, 0, ri0, ri1, &enersum, &virsum);
 471         energy.dispersion -= enersum;
 472         virial.dispersion -= virsum;
 473         integrate_table(vdwtab, scale, 4, ri0, ri1, &enersum, &virsum);
 474         energy.repulsion -= enersum;
 475         virial.repulsion -= virsum;
 476
 477
 478         /* Alright: Above we compensated by REMOVING the parts outside r0
 479          * corresponding to the ideal VdW 1/r6 and /r12 potentials.
 480          *
 481          * Regardless of whether r0 is the point where we start switching,
 482          * or the cutoff where we calculated the constant shift, we include
 483          * all the parts we are missing out to infinity from r0 by
 484          * calculating the analytical dispersion correction.
 485          */
 486         addCorrectionBeyondCutoff(&energy, &virial, r0);
 487     }
 488     else if (ic.vdwtype == VanDerWaalsType::Cut || EVDW_PME(ic.vdwtype)
 489              || ic.vdwtype == VanDerWaalsType::User)
 490     {
 491         /* Note that with LJ-PME, the dispersion correction is multiplied
 492          * by the difference between the actual C6 and the value of C6
 493          * that would produce the combination rule.
 494          * This means the normal energy and virial difference formulas
 495          * can be used here.
 496          */
 497
 498         const double rc3 = ic.rvdw * ic.rvdw * ic.rvdw;
 499         const double rc9 = rc3 * rc3 * rc3;
 500         if (ic.vdw_modifier == InteractionModifiers::PotShift)
 501         {
 502             /* Contribution within the cut-off */
 503             energy.dispersion += -4.0 * M_PI / (3.0 * rc3);
 504             energy.repulsion += 4.0 * M_PI / (3.0 * rc9);
 505         }
 506         /* Contribution beyond the cut-off */
 507         addCorrectionBeyondCutoff(&energy, &virial, ic.rvdw);
 508     }
 509     else
 510     {
 511         gmx_fatal(FARGS,
 512                   "Dispersion correction is not implemented for vdw-type = %s",
 513                   enumValueToString(ic.vdwtype));
 514     }
 515
 516     iParams->enerdiffsix_    = energy.dispersion;
 517     iParams->enerdifftwelve_ = energy.repulsion;
 518     /* The 0.5 is due to the Gromacs definition of the virial */
 519     iParams->virdiffsix_    = 0.5 * virial.dispersion;
 520     iParams->virdifftwelve_ = 0.5 * virial.repulsion;
 521 }
 522
 523 DispersionCorrection::DispersionCorrection(const gmx_mtop_t&          mtop,
 524                                            const t_inputrec&          inputrec,
 525                                            bool                       useBuckingham,
 526                                            int                        numAtomTypes,
 527                                            gmx::ArrayRef<const real>  nonbondedForceParameters,
 528                                            const interaction_const_t& ic,
 529                                            const char*                tableFileName) :
 530     eDispCorr_(inputrec.eDispCorr),
 531     vdwType_(inputrec.vdwtype),
 532     eFep_(inputrec.efep),
 533     topParams_(mtop, inputrec, useBuckingham, numAtomTypes, nonbondedForceParameters)
 534 {
 535     if (eDispCorr_ != DispersionCorrectionType::No)
 536     {
 537         GMX_RELEASE_ASSERT(tableFileName, "Need a table file name");
 538
 539         setInteractionParameters(&iParams_, ic, tableFileName);
 540     }
 541 }
 542
 543 bool DispersionCorrection::correctFullInteraction() const
 544 {
 545     return (eDispCorr_ == DispersionCorrectionType::AllEner
 546             || eDispCorr_ == DispersionCorrectionType::AllEnerPres);
 547 }
 548
 549 void DispersionCorrection::print(const gmx::MDLogger& mdlog) const
 550 {
 551     if (topParams_.avcsix_[0] == 0 && topParams_.avctwelve_[0] == 0)
 552     {
 553         GMX_LOG(mdlog.warning)
 554                 .asParagraph()
 555                 .appendText("WARNING: There are no atom pairs for dispersion correction");
 556     }
 557     else if (vdwType_ == VanDerWaalsType::User)
 558     {
 559         GMX_LOG(mdlog.warning)
 560                 .asParagraph()
 561                 .appendText("WARNING: using dispersion correction with user tables\n");
 562     }
 563
 564     std::string text = gmx::formatString("Long Range LJ corr.: <C6> %10.4e", topParams_.avcsix_[0]);
 565     if (correctFullInteraction())
 566     {
 567         text += gmx::formatString(" <C12> %10.4e", topParams_.avctwelve_[0]);
 568     }
 569     GMX_LOG(mdlog.info).appendText(text);
 570 }
 571
 572 void DispersionCorrection::setParameters(const interaction_const_t& ic)
 573 {
 574     if (eDispCorr_ != DispersionCorrectionType::No)
 575     {
 576         setInteractionParameters(&iParams_, ic, nullptr);
 577     }
 578 }
 579
 580 DispersionCorrection::Correction DispersionCorrection::calculate(const matrix box, const real lambda) const
 581
 582 {
 583     Correction corr;
 584
 585     if (eDispCorr_ == DispersionCorrectionType::No)
 586     {
 587         return corr;
 588     }
 589
 590     const bool bCorrAll  = correctFullInteraction();
 591     const bool bCorrPres = (eDispCorr_ == DispersionCorrectionType::EnerPres
 592                             || eDispCorr_ == DispersionCorrectionType::AllEnerPres);
 593
 594     const real invvol  = 1 / det(box);
 595     const real density = topParams_.numAtomsForDensity_ * invvol;
 596     const real numCorr = topParams_.numCorrections_;
 597
 598     real avcsix;
 599     real avctwelve;
 600     if (eFep_ == FreeEnergyPerturbationType::No)
 601     {
 602         avcsix    = topParams_.avcsix_[0];
 603         avctwelve = topParams_.avctwelve_[0];
 604     }
 605     else
 606     {
 607         avcsix    = (1 - lambda) * topParams_.avcsix_[0] + lambda * topParams_.avcsix_[1];
 608         avctwelve = (1 - lambda) * topParams_.avctwelve_[0] + lambda * topParams_.avctwelve_[1];
 609     }
 610
 611     const real enerdiff = numCorr * (density * iParams_.enerdiffsix_ - iParams_.enershiftsix_);
 612     corr.energy += avcsix * enerdiff;
 613     real dvdlambda = 0;
 614     if (eFep_ != FreeEnergyPerturbationType::No)
 615     {
 616         dvdlambda += (topParams_.avcsix_[1] - topParams_.avcsix_[0]) * enerdiff;
 617     }
 618     if (bCorrAll)
 619     {
 620         const real enerdiff = numCorr * (density * iParams_.enerdifftwelve_ - iParams_.enershifttwelve_);
 621         corr.energy += avctwelve * enerdiff;
 622         if (eFep_ != FreeEnergyPerturbationType::No)
 623         {
 624             dvdlambda += (topParams_.avctwelve_[1] - topParams_.avctwelve_[0]) * enerdiff;
 625         }
 626     }
 627
 628     if (bCorrPres)
 629     {
 630         corr.virial = numCorr * density * avcsix * iParams_.virdiffsix_ / 3.0;
 631         if (eDispCorr_ == DispersionCorrectionType::AllEnerPres)
 632         {
 633             corr.virial += numCorr * density * avctwelve * iParams_.virdifftwelve_ / 3.0;
 634         }
 635         /* The factor 2 is because of the Gromacs virial definition */
 636         corr.pressure = -2.0 * invvol * corr.virial * gmx::c_presfac;
 637     }
 638
 639     if (eFep_ != FreeEnergyPerturbationType::No)
 640     {
 641         corr.dvdl += dvdlambda;
 642     }
 643
 644     return corr;
 645 }