2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 #include "dispersioncorrection.h"
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/mdtypes/nblist.h"
49 #include "gromacs/tables/forcetable.h"
50 #include "gromacs/topology/mtop_util.h"
51 #include "gromacs/topology/topology.h"
52 #include "gromacs/utility/arrayref.h"
53 #include "gromacs/utility/fatalerror.h"
54 #include "gromacs/utility/logger.h"
56 /* Implementation here to avoid other files needing to include the file that defines t_nblists */
57 DispersionCorrection::InteractionParams::~InteractionParams() = default;
59 /* Returns a matrix, as flat list, of combination rule combined LJ parameters */
60 static std::vector<real> mk_nbfp_combination_rule(const gmx_ffparams_t& ffparams,
61 const CombinationRule comb_rule)
63 const int atnr = ffparams.atnr;
65 std::vector<real> nbfp(atnr * atnr * 2);
67 for (int i = 0; i < atnr; ++i)
69 for (int j = 0; j < atnr; ++j)
71 const real c6i = ffparams.iparams[i * (atnr + 1)].lj.c6;
72 const real c12i = ffparams.iparams[i * (atnr + 1)].lj.c12;
73 const real c6j = ffparams.iparams[j * (atnr + 1)].lj.c6;
74 const real c12j = ffparams.iparams[j * (atnr + 1)].lj.c12;
75 real c6 = std::sqrt(c6i * c6j);
76 real c12 = std::sqrt(c12i * c12j);
77 if (comb_rule == CombinationRule::Arithmetic && !gmx_numzero(c6) && !gmx_numzero(c12))
79 const real sigmai = gmx::sixthroot(c12i / c6i);
80 const real sigmaj = gmx::sixthroot(c12j / c6j);
81 const real epsi = c6i * c6i / c12i;
82 const real epsj = c6j * c6j / c12j;
83 const real sigma = 0.5 * (sigmai + sigmaj);
84 const real eps = std::sqrt(epsi * epsj);
85 c6 = eps * gmx::power6(sigma);
86 c12 = eps * gmx::power12(sigma);
88 C6(nbfp, atnr, i, j) = c6 * 6.0;
89 C12(nbfp, atnr, i, j) = c12 * 12.0;
96 /* Returns the A-topology atom type when aOrB=0, the B-topology atom type when aOrB=1 */
97 static int atomtypeAOrB(const t_atom& atom, int aOrB)
109 DispersionCorrection::TopologyParams::TopologyParams(const gmx_mtop_t& mtop,
110 const t_inputrec& inputrec,
113 gmx::ArrayRef<const real> nonbondedForceParameters)
115 const int ntp = numAtomTypes;
116 const gmx_bool bBHAM = useBuckingham;
118 gmx::ArrayRef<const real> nbfp = nonbondedForceParameters;
119 std::vector<real> nbfp_comb;
120 /* For LJ-PME, we want to correct for the difference between the
121 * actual C6 values and the C6 values used by the LJ-PME based on
122 * combination rules. */
123 if (EVDW_PME(inputrec.vdwtype))
125 nbfp_comb = mk_nbfp_combination_rule(mtop.ffparams,
126 (inputrec.ljpme_combination_rule == LongRangeVdW::LB)
127 ? CombinationRule::Arithmetic
128 : CombinationRule::Geometric);
129 for (int tpi = 0; tpi < ntp; ++tpi)
131 for (int tpj = 0; tpj < ntp; ++tpj)
133 C6(nbfp_comb, ntp, tpi, tpj) = C6(nbfp, ntp, tpi, tpj) - C6(nbfp_comb, ntp, tpi, tpj);
134 C12(nbfp_comb, ntp, tpi, tpj) = C12(nbfp, ntp, tpi, tpj);
140 for (int q = 0; q < (inputrec.efep == FreeEnergyPerturbationType::No ? 1 : 2); q++)
146 if (!EI_TPI(inputrec.eI))
148 numAtomsForDensity_ = mtop.natoms;
149 numCorrections_ = 0.5 * mtop.natoms;
151 /* Count the types so we avoid natoms^2 operations */
152 std::vector<int> typecount(ntp);
153 gmx_mtop_count_atomtypes(mtop, q, typecount.data());
155 for (int tpi = 0; tpi < ntp; tpi++)
157 for (int tpj = tpi; tpj < ntp; tpj++)
159 const int64_t iCount = typecount[tpi];
160 const int64_t jCount = typecount[tpj];
164 npair_ij = iCount * jCount;
168 npair_ij = iCount * (iCount - 1) / 2;
172 /* nbfp now includes the 6.0 derivative prefactor */
173 csix += npair_ij * BHAMC(nbfp, ntp, tpi, tpj) / 6.0;
177 /* nbfp now includes the 6.0/12.0 derivative prefactors */
178 csix += npair_ij * C6(nbfp, ntp, tpi, tpj) / 6.0;
179 ctwelve += npair_ij * C12(nbfp, ntp, tpi, tpj) / 12.0;
184 /* Subtract the excluded pairs.
185 * The main reason for substracting exclusions is that in some cases
186 * some combinations might never occur and the parameters could have
187 * any value. These unused values should not influence the dispersion
190 for (const gmx_molblock_t& molb : mtop.molblock)
192 const int nmol = molb.nmol;
193 const t_atoms* atoms = &mtop.moltype[molb.type].atoms;
194 const auto& excl = mtop.moltype[molb.type].excls;
195 for (int i = 0; (i < atoms->nr); i++)
197 const int tpi = atomtypeAOrB(atoms->atom[i], q);
198 for (const int k : excl[i])
202 const int tpj = atomtypeAOrB(atoms->atom[k], q);
205 /* nbfp now includes the 6.0 derivative prefactor */
206 csix -= nmol * BHAMC(nbfp, ntp, tpi, tpj) / 6.0;
210 /* nbfp now includes the 6.0/12.0 derivative prefactors */
211 csix -= nmol * C6(nbfp, ntp, tpi, tpj) / 6.0;
212 ctwelve -= nmol * C12(nbfp, ntp, tpi, tpj) / 12.0;
222 const t_atoms& atoms_tpi = mtop.moltype[mtop.molblock.back().type].atoms;
224 /* Only correct for the interaction of the test particle
225 * with the rest of the system.
227 numAtomsForDensity_ = mtop.natoms - atoms_tpi.nr;
228 numCorrections_ = atoms_tpi.nr;
231 for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
233 const gmx_molblock_t& molb = mtop.molblock[mb];
234 const t_atoms& atoms = mtop.moltype[molb.type].atoms;
235 for (int j = 0; j < atoms.nr; j++)
237 int nmolc = molb.nmol;
238 /* Remove the interaction of the test charge group
241 if (mb == mtop.molblock.size() - 1)
245 if (mb == 0 && molb.nmol == 1)
248 "Old format tpr with TPI, please generate a new tpr file");
251 const int tpj = atomtypeAOrB(atoms.atom[j], q);
252 for (int i = 0; i < atoms_tpi.nr; i++)
254 const int tpi = atomtypeAOrB(atoms_tpi.atom[i], q);
257 /* nbfp now includes the 6.0 derivative prefactor */
258 csix += nmolc * BHAMC(nbfp, ntp, tpi, tpj) / 6.0;
262 /* nbfp now includes the 6.0/12.0 derivative prefactors */
263 csix += nmolc * C6(nbfp, ntp, tpi, tpj) / 6.0;
264 ctwelve += nmolc * C12(nbfp, ntp, tpi, tpj) / 12.0;
271 if (npair - nexcl <= 0)
278 csix /= npair - nexcl;
279 ctwelve /= npair - nexcl;
283 fprintf(debug, "Counted %" PRId64 " exclusions\n", nexcl);
284 fprintf(debug, "Average C6 parameter is: %10g\n", csix);
285 fprintf(debug, "Average C12 parameter is: %10g\n", ctwelve);
288 avctwelve_[q] = ctwelve;
292 static void integrate_table(const real vdwtab[],
300 const double invscale = 1.0 / scale;
301 const double invscale2 = invscale * invscale;
302 const double invscale3 = invscale * invscale2;
304 /* Following summation derived from cubic spline definition,
305 * Numerical Recipies in C, second edition, p. 113-116. Exact for
306 * the cubic spline. We first calculate the negative of the
307 * energy from rvdw to rvdw_switch, assuming that g(r)=1, and then
308 * add the more standard, abrupt cutoff correction to that result,
309 * yielding the long-range correction for a switched function. We
310 * perform both the pressure and energy loops at the same time for
311 * simplicity, as the computational cost is low. */
313 /* Since the dispersion table has been scaled down a factor
314 * 6.0 and the repulsion a factor 12.0 to compensate for the
315 * c6/c12 parameters inside nbfp[] being scaled up (to save
316 * flops in kernels), we need to correct for this.
318 const double tabfactor = (offstart == 0 ? 6.0 : 12.0);
320 double enersum = 0.0;
322 for (int ri = rstart; ri < rend; ++ri)
324 const double r = ri * invscale;
325 const double ea = invscale3;
326 const double eb = 2.0 * invscale2 * r;
327 const double ec = invscale * r * r;
329 const double pa = invscale3;
330 const double pb = 3.0 * invscale2 * r;
331 const double pc = 3.0 * invscale * r * r;
332 const double pd = r * r * r;
334 /* this "8" is from the packing in the vdwtab array - perhaps
335 should be defined? */
337 const int offset = 8 * ri + offstart;
338 const double y0 = vdwtab[offset];
339 const double f = vdwtab[offset + 1];
340 const double g = vdwtab[offset + 2];
341 const double h = vdwtab[offset + 3];
343 enersum += y0 * (ea / 3 + eb / 2 + ec) + f * (ea / 4 + eb / 3 + ec / 2)
344 + g * (ea / 5 + eb / 4 + ec / 3) + h * (ea / 6 + eb / 5 + ec / 4);
345 virsum += f * (pa / 4 + pb / 3 + pc / 2 + pd) + 2 * g * (pa / 5 + pb / 4 + pc / 3 + pd / 2)
346 + 3 * h * (pa / 6 + pb / 5 + pc / 4 + pd / 3);
348 *enerout = 4.0 * M_PI * enersum * tabfactor;
349 *virout = 4.0 * M_PI * virsum * tabfactor;
352 /* Struct for storing and passing energy or virial corrections */
353 struct InteractionCorrection
359 /* Adds the energy and virial corrections beyond the cut-off */
360 static void addCorrectionBeyondCutoff(InteractionCorrection* energy,
361 InteractionCorrection* virial,
362 const double cutoffDistance)
364 const double rc3 = cutoffDistance * cutoffDistance * cutoffDistance;
365 const double rc9 = rc3 * rc3 * rc3;
367 energy->dispersion += -4.0 * M_PI / (3.0 * rc3);
368 energy->repulsion += 4.0 * M_PI / (9.0 * rc9);
369 virial->dispersion += 8.0 * M_PI / rc3;
370 virial->repulsion += -16.0 * M_PI / (3.0 * rc9);
373 void DispersionCorrection::setInteractionParameters(InteractionParams* iParams,
374 const interaction_const_t& ic,
375 const char* tableFileName)
377 /* We only need to set the tables at first call, i.e. tableFileName!=nullptr
378 * or when we changed the cut-off with LJ-PME tuning.
380 if (tableFileName || EVDW_PME(ic.vdwtype))
382 iParams->dispersionCorrectionTable_ =
383 makeDispersionCorrectionTable(nullptr, &ic, ic.rvdw, tableFileName);
386 InteractionCorrection energy;
387 InteractionCorrection virial;
389 if ((ic.vdw_modifier == InteractionModifiers::PotShift)
390 || (ic.vdw_modifier == InteractionModifiers::PotSwitch)
391 || (ic.vdw_modifier == InteractionModifiers::ForceSwitch)
392 || (ic.vdwtype == VanDerWaalsType::Shift) || (ic.vdwtype == VanDerWaalsType::Switch))
394 if (((ic.vdw_modifier == InteractionModifiers::PotSwitch)
395 || (ic.vdw_modifier == InteractionModifiers::ForceSwitch)
396 || (ic.vdwtype == VanDerWaalsType::Switch))
397 && ic.rvdw_switch == 0)
400 "With dispersion correction rvdw-switch can not be zero "
402 enumValueToString(ic.vdwtype));
405 GMX_ASSERT(iParams->dispersionCorrectionTable_, "We need an initialized table");
407 /* TODO This code depends on the logic in tables.c that
408 constructs the table layout, which should be made
409 explicit in future cleanup. */
410 GMX_ASSERT(iParams->dispersionCorrectionTable_->interaction == GMX_TABLE_INTERACTION_VDWREP_VDWDISP,
411 "Dispersion-correction code needs a table with both repulsion and dispersion "
413 const real scale = iParams->dispersionCorrectionTable_->scale;
414 const real* vdwtab = iParams->dispersionCorrectionTable_->data.data();
416 /* Round the cut-offs to exact table values for precision */
417 int ri0 = static_cast<int>(std::floor(ic.rvdw_switch * scale));
418 int ri1 = static_cast<int>(std::ceil(ic.rvdw * scale));
420 /* The code below has some support for handling force-switching, i.e.
421 * when the force (instead of potential) is switched over a limited
422 * region. This leads to a constant shift in the potential inside the
423 * switching region, which we can handle by adding a constant energy
424 * term in the force-switch case just like when we do potential-shift.
426 * For now this is not enabled, but to keep the functionality in the
427 * code we check separately for switch and shift. When we do force-switch
428 * the shifting point is rvdw_switch, while it is the cutoff when we
429 * have a classical potential-shift.
431 * For a pure potential-shift the potential has a constant shift
432 * all the way out to the cutoff, and that is it. For other forms
433 * we need to calculate the constant shift up to the point where we
434 * start modifying the potential.
436 ri0 = (ic.vdw_modifier == InteractionModifiers::PotShift) ? ri1 : ri0;
438 const double r0 = ri0 / scale;
439 const double rc3 = r0 * r0 * r0;
440 const double rc9 = rc3 * rc3 * rc3;
442 if ((ic.vdw_modifier == InteractionModifiers::ForceSwitch) || (ic.vdwtype == VanDerWaalsType::Shift))
444 /* Determine the constant energy shift below rvdw_switch.
445 * Table has a scale factor since we have scaled it down to compensate
446 * for scaling-up c6/c12 with the derivative factors to save flops in analytical kernels.
448 iParams->enershiftsix_ = static_cast<real>(-1.0 / (rc3 * rc3)) - 6.0 * vdwtab[8 * ri0];
449 iParams->enershifttwelve_ = static_cast<real>(1.0 / (rc9 * rc3)) - 12.0 * vdwtab[8 * ri0 + 4];
451 else if (ic.vdw_modifier == InteractionModifiers::PotShift)
453 iParams->enershiftsix_ = static_cast<real>(-1.0 / (rc3 * rc3));
454 iParams->enershifttwelve_ = static_cast<real>(1.0 / (rc9 * rc3));
457 /* Add the constant part from 0 to rvdw_switch.
458 * This integration from 0 to rvdw_switch overcounts the number
459 * of interactions by 1, as it also counts the self interaction.
460 * We will correct for this later.
462 energy.dispersion += 4.0 * M_PI * iParams->enershiftsix_ * rc3 / 3.0;
463 energy.repulsion += 4.0 * M_PI * iParams->enershifttwelve_ * rc3 / 3.0;
465 /* Calculate the contribution in the range [r0,r1] where we
466 * modify the potential. For a pure potential-shift modifier we will
467 * have ri0==ri1, and there will not be any contribution here.
471 integrate_table(vdwtab, scale, 0, ri0, ri1, &enersum, &virsum);
472 energy.dispersion -= enersum;
473 virial.dispersion -= virsum;
474 integrate_table(vdwtab, scale, 4, ri0, ri1, &enersum, &virsum);
475 energy.repulsion -= enersum;
476 virial.repulsion -= virsum;
479 /* Alright: Above we compensated by REMOVING the parts outside r0
480 * corresponding to the ideal VdW 1/r6 and /r12 potentials.
482 * Regardless of whether r0 is the point where we start switching,
483 * or the cutoff where we calculated the constant shift, we include
484 * all the parts we are missing out to infinity from r0 by
485 * calculating the analytical dispersion correction.
487 addCorrectionBeyondCutoff(&energy, &virial, r0);
489 else if (ic.vdwtype == VanDerWaalsType::Cut || EVDW_PME(ic.vdwtype)
490 || ic.vdwtype == VanDerWaalsType::User)
492 /* Note that with LJ-PME, the dispersion correction is multiplied
493 * by the difference between the actual C6 and the value of C6
494 * that would produce the combination rule.
495 * This means the normal energy and virial difference formulas
499 const double rc3 = ic.rvdw * ic.rvdw * ic.rvdw;
500 const double rc9 = rc3 * rc3 * rc3;
501 if (ic.vdw_modifier == InteractionModifiers::PotShift)
503 /* Contribution within the cut-off */
504 energy.dispersion += -4.0 * M_PI / (3.0 * rc3);
505 energy.repulsion += 4.0 * M_PI / (3.0 * rc9);
507 /* Contribution beyond the cut-off */
508 addCorrectionBeyondCutoff(&energy, &virial, ic.rvdw);
513 "Dispersion correction is not implemented for vdw-type = %s",
514 enumValueToString(ic.vdwtype));
517 iParams->enerdiffsix_ = energy.dispersion;
518 iParams->enerdifftwelve_ = energy.repulsion;
519 /* The 0.5 is due to the Gromacs definition of the virial */
520 iParams->virdiffsix_ = 0.5 * virial.dispersion;
521 iParams->virdifftwelve_ = 0.5 * virial.repulsion;
524 DispersionCorrection::DispersionCorrection(const gmx_mtop_t& mtop,
525 const t_inputrec& inputrec,
528 gmx::ArrayRef<const real> nonbondedForceParameters,
529 const interaction_const_t& ic,
530 const char* tableFileName) :
531 eDispCorr_(inputrec.eDispCorr),
532 vdwType_(inputrec.vdwtype),
533 eFep_(inputrec.efep),
534 topParams_(mtop, inputrec, useBuckingham, numAtomTypes, nonbondedForceParameters)
536 if (eDispCorr_ != DispersionCorrectionType::No)
538 GMX_RELEASE_ASSERT(tableFileName, "Need a table file name");
540 setInteractionParameters(&iParams_, ic, tableFileName);
544 bool DispersionCorrection::correctFullInteraction() const
546 return (eDispCorr_ == DispersionCorrectionType::AllEner
547 || eDispCorr_ == DispersionCorrectionType::AllEnerPres);
550 void DispersionCorrection::print(const gmx::MDLogger& mdlog) const
552 if (topParams_.avcsix_[0] == 0 && topParams_.avctwelve_[0] == 0)
554 GMX_LOG(mdlog.warning)
556 .appendText("WARNING: There are no atom pairs for dispersion correction");
558 else if (vdwType_ == VanDerWaalsType::User)
560 GMX_LOG(mdlog.warning)
562 .appendText("WARNING: using dispersion correction with user tables\n");
565 std::string text = gmx::formatString("Long Range LJ corr.: <C6> %10.4e", topParams_.avcsix_[0]);
566 if (correctFullInteraction())
568 text += gmx::formatString(" <C12> %10.4e", topParams_.avctwelve_[0]);
570 GMX_LOG(mdlog.info).appendText(text);
573 void DispersionCorrection::setParameters(const interaction_const_t& ic)
575 if (eDispCorr_ != DispersionCorrectionType::No)
577 setInteractionParameters(&iParams_, ic, nullptr);
581 DispersionCorrection::Correction DispersionCorrection::calculate(const matrix box, const real lambda) const
586 if (eDispCorr_ == DispersionCorrectionType::No)
591 const bool bCorrAll = correctFullInteraction();
592 const bool bCorrPres = (eDispCorr_ == DispersionCorrectionType::EnerPres
593 || eDispCorr_ == DispersionCorrectionType::AllEnerPres);
595 const real invvol = 1 / det(box);
596 const real density = topParams_.numAtomsForDensity_ * invvol;
597 const real numCorr = topParams_.numCorrections_;
601 if (eFep_ == FreeEnergyPerturbationType::No)
603 avcsix = topParams_.avcsix_[0];
604 avctwelve = topParams_.avctwelve_[0];
608 avcsix = (1 - lambda) * topParams_.avcsix_[0] + lambda * topParams_.avcsix_[1];
609 avctwelve = (1 - lambda) * topParams_.avctwelve_[0] + lambda * topParams_.avctwelve_[1];
612 const real enerdiff = numCorr * (density * iParams_.enerdiffsix_ - iParams_.enershiftsix_);
613 corr.energy += avcsix * enerdiff;
615 if (eFep_ != FreeEnergyPerturbationType::No)
617 dvdlambda += (topParams_.avcsix_[1] - topParams_.avcsix_[0]) * enerdiff;
621 const real enerdiff = numCorr * (density * iParams_.enerdifftwelve_ - iParams_.enershifttwelve_);
622 corr.energy += avctwelve * enerdiff;
623 if (eFep_ != FreeEnergyPerturbationType::No)
625 dvdlambda += (topParams_.avctwelve_[1] - topParams_.avctwelve_[0]) * enerdiff;
631 corr.virial = numCorr * density * avcsix * iParams_.virdiffsix_ / 3.0;
632 if (eDispCorr_ == DispersionCorrectionType::AllEnerPres)
634 corr.virial += numCorr * density * avctwelve * iParams_.virdifftwelve_ / 3.0;
636 /* The factor 2 is because of the Gromacs virial definition */
637 corr.pressure = -2.0 * invvol * corr.virial * gmx::c_presfac;
640 if (eFep_ != FreeEnergyPerturbationType::No)
642 corr.dvdl += dvdlambda;