[alexxy/gromacs.git] / src / gromacs / mdlib / dispersioncorrection.cpp

/*
 * This file is part of the GROMACS molecular simulation package.
 *
 * Copyright (c) 2019, by the GROMACS development team, led by
 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 * and including many others, as listed in the AUTHORS file in the
 * top-level source directory and at http://www.gromacs.org.
 *
 * GROMACS is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 *
 * GROMACS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with GROMACS; if not, see
 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 *
 * If you want to redistribute modifications to GROMACS, please
 * consider that scientific software is very special. Version
 * control is crucial - bugs must be traceable. We will be happy to
 * consider code for inclusion in the official distribution, but
 * derived work must not be called official GROMACS. Details are found
 * in the README & COPYING files - if they are missing, get the
 * official version at http://www.gromacs.org.
 *
 * To help us fund GROMACS development, we humbly ask that you cite
 * the research papers on the package. Check out http://www.gromacs.org.
 */
#include "gmxpre.h"

#include "dispersioncorrection.h"

#include <cstdio>

#include "gromacs/math/units.h"
#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdtypes/forcerec.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/nblist.h"
#include "gromacs/tables/forcetable.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/topology/topology.h"
#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/logger.h"

/* Implementation here to avoid other files needing to include the file that defines t_nblists */
DispersionCorrection::InteractionParams::~InteractionParams() = default;

/* Returns a matrix, as flat list, of combination rule combined LJ parameters */
static std::vector<real> mk_nbfp_combination_rule(const gmx_ffparams_t& ffparams, const int comb_rule)
{
    const int atnr = ffparams.atnr;

    std::vector<real> nbfp(atnr * atnr * 2);

    for (int i = 0; i < atnr; ++i)
    {
        for (int j = 0; j < atnr; ++j)
        {
            const real c6i  = ffparams.iparams[i * (atnr + 1)].lj.c6;
            const real c12i = ffparams.iparams[i * (atnr + 1)].lj.c12;
            const real c6j  = ffparams.iparams[j * (atnr + 1)].lj.c6;
            const real c12j = ffparams.iparams[j * (atnr + 1)].lj.c12;
            real       c6   = std::sqrt(c6i * c6j);
            real       c12  = std::sqrt(c12i * c12j);
            if (comb_rule == eCOMB_ARITHMETIC && !gmx_numzero(c6) && !gmx_numzero(c12))
            {
                const real sigmai = gmx::sixthroot(c12i / c6i);
                const real sigmaj = gmx::sixthroot(c12j / c6j);
                const real epsi   = c6i * c6i / c12i;
                const real epsj   = c6j * c6j / c12j;
                const real sigma  = 0.5 * (sigmai + sigmaj);
                const real eps    = std::sqrt(epsi * epsj);
                c6                = eps * gmx::power6(sigma);
                c12               = eps * gmx::power12(sigma);
            }
            C6(nbfp, atnr, i, j)  = c6 * 6.0;
            C12(nbfp, atnr, i, j) = c12 * 12.0;
        }
    }

    return nbfp;
}

/* Returns the A-topology atom type when aOrB=0, the B-topology atom type when aOrB=1 */
static int atomtypeAOrB(const t_atom& atom, int aOrB)
{
    if (aOrB == 0)
    {
        return atom.type;
    }
    else
    {
        return atom.typeB;
    }
}

DispersionCorrection::TopologyParams::TopologyParams(const gmx_mtop_t&         mtop,
                                                     const t_inputrec&         inputrec,
                                                     bool                      useBuckingham,
                                                     int                       numAtomTypes,
                                                     gmx::ArrayRef<const real> nonbondedForceParameters)
{
    const int      ntp   = numAtomTypes;
    const gmx_bool bBHAM = useBuckingham;

    gmx::ArrayRef<const real> nbfp = nonbondedForceParameters;
    std::vector<real>         nbfp_comb;
    /* For LJ-PME, we want to correct for the difference between the
     * actual C6 values and the C6 values used by the LJ-PME based on
     * combination rules. */
    if (EVDW_PME(inputrec.vdwtype))
    {
        nbfp_comb = mk_nbfp_combination_rule(
                mtop.ffparams,
                (inputrec.ljpme_combination_rule == eljpmeLB) ? eCOMB_ARITHMETIC : eCOMB_GEOMETRIC);
        for (int tpi = 0; tpi < ntp; ++tpi)
        {
            for (int tpj = 0; tpj < ntp; ++tpj)
            {
                C6(nbfp_comb, ntp, tpi, tpj) = C6(nbfp, ntp, tpi, tpj) - C6(nbfp_comb, ntp, tpi, tpj);
                C12(nbfp_comb, ntp, tpi, tpj) = C12(nbfp, ntp, tpi, tpj);
            }
        }
        nbfp = nbfp_comb;
    }

    for (int q = 0; q < (inputrec.efep == efepNO ? 1 : 2); q++)
    {
        double csix    = 0;
        double ctwelve = 0;
        int    npair   = 0;
        int    nexcl   = 0;
        if (!EI_TPI(inputrec.eI))
        {
            numAtomsForDensity_ = mtop.natoms;
            numCorrections_     = 0.5 * mtop.natoms;

            /* Count the types so we avoid natoms^2 operations */
            std::vector<int> typecount(ntp);
            gmx_mtop_count_atomtypes(&mtop, q, typecount.data());

            for (int tpi = 0; tpi < ntp; tpi++)
            {
                for (int tpj = tpi; tpj < ntp; tpj++)
                {
                    const int iCount = typecount[tpi];
                    const int jCount = typecount[tpj];
                    int       npair_ij;
                    if (tpi != tpj)
                    {
                        npair_ij = iCount * jCount;
                    }
                    else
                    {
                        npair_ij = iCount * (iCount - 1) / 2;
                    }
                    if (bBHAM)
                    {
                        /* nbfp now includes the 6.0 derivative prefactor */
                        csix += npair_ij * BHAMC(nbfp, ntp, tpi, tpj) / 6.0;
                    }
                    else
                    {
                        /* nbfp now includes the 6.0/12.0 derivative prefactors */
                        csix += npair_ij * C6(nbfp, ntp, tpi, tpj) / 6.0;
                        ctwelve += npair_ij * C12(nbfp, ntp, tpi, tpj) / 12.0;
                    }
                    npair += npair_ij;
                }
            }
            /* Subtract the excluded pairs.
             * The main reason for substracting exclusions is that in some cases
             * some combinations might never occur and the parameters could have
             * any value. These unused values should not influence the dispersion
             * correction.
             */
            for (const gmx_molblock_t& molb : mtop.molblock)
            {
                const int       nmol  = molb.nmol;
                const t_atoms*  atoms = &mtop.moltype[molb.type].atoms;
                const t_blocka* excl  = &mtop.moltype[molb.type].excls;
                for (int i = 0; (i < atoms->nr); i++)
                {
                    const int tpi = atomtypeAOrB(atoms->atom[i], q);
                    const int j1  = excl->index[i];
                    const int j2  = excl->index[i + 1];
                    for (int j = j1; j < j2; j++)
                    {
                        const int k = excl->a[j];
                        if (k > i)
                        {
                            const int tpj = atomtypeAOrB(atoms->atom[k], q);
                            if (bBHAM)
                            {
                                /* nbfp now includes the 6.0 derivative prefactor */
                                csix -= nmol * BHAMC(nbfp, ntp, tpi, tpj) / 6.0;
                            }
                            else
                            {
                                /* nbfp now includes the 6.0/12.0 derivative prefactors */
                                csix -= nmol * C6(nbfp, ntp, tpi, tpj) / 6.0;
                                ctwelve -= nmol * C12(nbfp, ntp, tpi, tpj) / 12.0;
                            }
                            nexcl += molb.nmol;
                        }
                    }
                }
            }
        }
        else
        {
            const t_atoms& atoms_tpi = mtop.moltype[mtop.molblock.back().type].atoms;

            /* Only correct for the interaction of the test particle
             * with the rest of the system.
             */
            numAtomsForDensity_ = mtop.natoms - atoms_tpi.nr;
            numCorrections_     = atoms_tpi.nr;

            npair = 0;
            for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
            {
                const gmx_molblock_t& molb  = mtop.molblock[mb];
                const t_atoms&        atoms = mtop.moltype[molb.type].atoms;
                for (int j = 0; j < atoms.nr; j++)
                {
                    int nmolc = molb.nmol;
                    /* Remove the interaction of the test charge group
                     * with itself.
                     */
                    if (mb == mtop.molblock.size() - 1)
                    {
                        nmolc--;

                        if (mb == 0 && molb.nmol == 1)
                        {
                            gmx_fatal(FARGS,
                                      "Old format tpr with TPI, please generate a new tpr file");
                        }
                    }
                    const int tpj = atomtypeAOrB(atoms.atom[j], q);
                    for (int i = 0; i < atoms_tpi.nr; i++)
                    {
                        const int tpi = atomtypeAOrB(atoms_tpi.atom[i], q);
                        if (bBHAM)
                        {
                            /* nbfp now includes the 6.0 derivative prefactor */
                            csix += nmolc * BHAMC(nbfp, ntp, tpi, tpj) / 6.0;
                        }
                        else
                        {
                            /* nbfp now includes the 6.0/12.0 derivative prefactors */
                            csix += nmolc * C6(nbfp, ntp, tpi, tpj) / 6.0;
                            ctwelve += nmolc * C12(nbfp, ntp, tpi, tpj) / 12.0;
                        }
                        npair += nmolc;
                    }
                }
            }
        }
        if (npair - nexcl <= 0)
        {
            csix    = 0;
            ctwelve = 0;
        }
        else
        {
            csix /= npair - nexcl;
            ctwelve /= npair - nexcl;
        }
        if (debug)
        {
            fprintf(debug, "Counted %d exclusions\n", nexcl);
            fprintf(debug, "Average C6 parameter is: %10g\n", csix);
            fprintf(debug, "Average C12 parameter is: %10g\n", ctwelve);
        }
        avcsix_[q]    = csix;
        avctwelve_[q] = ctwelve;
    }
}

static void integrate_table(const real vdwtab[],
                            const real scale,
                            const int  offstart,
                            const int  rstart,
                            const int  rend,
                            double*    enerout,
                            double*    virout)
{
    const double invscale  = 1.0 / scale;
    const double invscale2 = invscale * invscale;
    const double invscale3 = invscale * invscale2;

    /* Following summation derived from cubic spline definition,
     * Numerical Recipies in C, second edition, p. 113-116.  Exact for
     * the cubic spline.  We first calculate the negative of the
     * energy from rvdw to rvdw_switch, assuming that g(r)=1, and then
     * add the more standard, abrupt cutoff correction to that result,
     * yielding the long-range correction for a switched function.  We
     * perform both the pressure and energy loops at the same time for
     * simplicity, as the computational cost is low. */

    /* Since the dispersion table has been scaled down a factor
     * 6.0 and the repulsion a factor 12.0 to compensate for the
     * c6/c12 parameters inside nbfp[] being scaled up (to save
     * flops in kernels), we need to correct for this.
     */
    const double tabfactor = (offstart == 0 ? 6.0 : 12.0);

    double enersum = 0.0;
    double virsum  = 0.0;
    for (int ri = rstart; ri < rend; ++ri)
    {
        const double r  = ri * invscale;
        const double ea = invscale3;
        const double eb = 2.0 * invscale2 * r;
        const double ec = invscale * r * r;

        const double pa = invscale3;
        const double pb = 3.0 * invscale2 * r;
        const double pc = 3.0 * invscale * r * r;
        const double pd = r * r * r;

        /* this "8" is from the packing in the vdwtab array - perhaps
           should be defined? */

        const int    offset = 8 * ri + offstart;
        const double y0     = vdwtab[offset];
        const double f      = vdwtab[offset + 1];
        const double g      = vdwtab[offset + 2];
        const double h      = vdwtab[offset + 3];

        enersum += y0 * (ea / 3 + eb / 2 + ec) + f * (ea / 4 + eb / 3 + ec / 2)
                   + g * (ea / 5 + eb / 4 + ec / 3) + h * (ea / 6 + eb / 5 + ec / 4);
        virsum += f * (pa / 4 + pb / 3 + pc / 2 + pd) + 2 * g * (pa / 5 + pb / 4 + pc / 3 + pd / 2)
                  + 3 * h * (pa / 6 + pb / 5 + pc / 4 + pd / 3);
    }
    *enerout = 4.0 * M_PI * enersum * tabfactor;
    *virout  = 4.0 * M_PI * virsum * tabfactor;
}

/* Struct for storing and passing energy or virial corrections */
struct InteractionCorrection
{
    real dispersion = 0;
    real repulsion  = 0;
};

/* Adds the energy and virial corrections beyond the cut-off */
static void addCorrectionBeyondCutoff(InteractionCorrection* energy,
                                      InteractionCorrection* virial,
                                      const double           cutoffDistance)
{
    const double rc3 = cutoffDistance * cutoffDistance * cutoffDistance;
    const double rc9 = rc3 * rc3 * rc3;

    energy->dispersion += -4.0 * M_PI / (3.0 * rc3);
    energy->repulsion += 4.0 * M_PI / (9.0 * rc9);
    virial->dispersion += 8.0 * M_PI / rc3;
    virial->repulsion += -16.0 * M_PI / (3.0 * rc9);
}

void DispersionCorrection::setInteractionParameters(InteractionParams*         iParams,
                                                    const interaction_const_t& ic,
                                                    const char*                tableFileName)
{
    /* We only need to set the tables at first call, i.e. tableFileName!=nullptr
     * or when we changed the cut-off with LJ-PME tuning.
     */
    if (tableFileName || EVDW_PME(ic.vdwtype))
    {
        iParams->dispersionCorrectionTable_ =
                makeDispersionCorrectionTable(nullptr, &ic, ic.rvdw, tableFileName);
    }

    InteractionCorrection energy;
    InteractionCorrection virial;

    if ((ic.vdw_modifier == eintmodPOTSHIFT) || (ic.vdw_modifier == eintmodPOTSWITCH)
        || (ic.vdw_modifier == eintmodFORCESWITCH) || (ic.vdwtype == evdwSHIFT)
        || (ic.vdwtype == evdwSWITCH))
    {
        if (((ic.vdw_modifier == eintmodPOTSWITCH) || (ic.vdw_modifier == eintmodFORCESWITCH)
             || (ic.vdwtype == evdwSWITCH))
            && ic.rvdw_switch == 0)
        {
            gmx_fatal(FARGS,
                      "With dispersion correction rvdw-switch can not be zero "
                      "for vdw-type = %s",
                      evdw_names[ic.vdwtype]);
        }

        GMX_ASSERT(iParams->dispersionCorrectionTable_, "We need an initialized table");

        /* TODO This code depends on the logic in tables.c that
           constructs the table layout, which should be made
           explicit in future cleanup. */
        GMX_ASSERT(iParams->dispersionCorrectionTable_->interaction == GMX_TABLE_INTERACTION_VDWREP_VDWDISP,
                   "Dispersion-correction code needs a table with both repulsion and dispersion "
                   "terms");
        const real  scale  = iParams->dispersionCorrectionTable_->scale;
        const real* vdwtab = iParams->dispersionCorrectionTable_->data;

        /* Round the cut-offs to exact table values for precision */
        int ri0 = static_cast<int>(std::floor(ic.rvdw_switch * scale));
        int ri1 = static_cast<int>(std::ceil(ic.rvdw * scale));

        /* The code below has some support for handling force-switching, i.e.
         * when the force (instead of potential) is switched over a limited
         * region. This leads to a constant shift in the potential inside the
         * switching region, which we can handle by adding a constant energy
         * term in the force-switch case just like when we do potential-shift.
         *
         * For now this is not enabled, but to keep the functionality in the
         * code we check separately for switch and shift. When we do force-switch
         * the shifting point is rvdw_switch, while it is the cutoff when we
         * have a classical potential-shift.
         *
         * For a pure potential-shift the potential has a constant shift
         * all the way out to the cutoff, and that is it. For other forms
         * we need to calculate the constant shift up to the point where we
         * start modifying the potential.
         */
        ri0 = (ic.vdw_modifier == eintmodPOTSHIFT) ? ri1 : ri0;

        const double r0  = ri0 / scale;
        const double rc3 = r0 * r0 * r0;
        const double rc9 = rc3 * rc3 * rc3;

        if ((ic.vdw_modifier == eintmodFORCESWITCH) || (ic.vdwtype == evdwSHIFT))
        {
            /* Determine the constant energy shift below rvdw_switch.
             * Table has a scale factor since we have scaled it down to compensate
             * for scaling-up c6/c12 with the derivative factors to save flops in analytical kernels.
             */
            iParams->enershiftsix_ = static_cast<real>(-1.0 / (rc3 * rc3)) - 6.0 * vdwtab[8 * ri0];
            iParams->enershifttwelve_ = static_cast<real>(1.0 / (rc9 * rc3)) - 12.0 * vdwtab[8 * ri0 + 4];
        }
        else if (ic.vdw_modifier == eintmodPOTSHIFT)
        {
            iParams->enershiftsix_    = static_cast<real>(-1.0 / (rc3 * rc3));
            iParams->enershifttwelve_ = static_cast<real>(1.0 / (rc9 * rc3));
        }

        /* Add the constant part from 0 to rvdw_switch.
         * This integration from 0 to rvdw_switch overcounts the number
         * of interactions by 1, as it also counts the self interaction.
         * We will correct for this later.
         */
        energy.dispersion += 4.0 * M_PI * iParams->enershiftsix_ * rc3 / 3.0;
        energy.repulsion += 4.0 * M_PI * iParams->enershifttwelve_ * rc3 / 3.0;

        /* Calculate the contribution in the range [r0,r1] where we
         * modify the potential. For a pure potential-shift modifier we will
         * have ri0==ri1, and there will not be any contribution here.
         */
        double enersum = 0;
        double virsum  = 0;
        integrate_table(vdwtab, scale, 0, ri0, ri1, &enersum, &virsum);
        energy.dispersion -= enersum;
        virial.dispersion -= virsum;
        integrate_table(vdwtab, scale, 4, ri0, ri1, &enersum, &virsum);
        energy.repulsion -= enersum;
        virial.repulsion -= virsum;


        /* Alright: Above we compensated by REMOVING the parts outside r0
         * corresponding to the ideal VdW 1/r6 and /r12 potentials.
         *
         * Regardless of whether r0 is the point where we start switching,
         * or the cutoff where we calculated the constant shift, we include
         * all the parts we are missing out to infinity from r0 by
         * calculating the analytical dispersion correction.
         */
        addCorrectionBeyondCutoff(&energy, &virial, r0);
    }
    else if (ic.vdwtype == evdwCUT || EVDW_PME(ic.vdwtype) || ic.vdwtype == evdwUSER)
    {
        /* Note that with LJ-PME, the dispersion correction is multiplied
         * by the difference between the actual C6 and the value of C6
         * that would produce the combination rule.
         * This means the normal energy and virial difference formulas
         * can be used here.
         */

        const double rc3 = ic.rvdw * ic.rvdw * ic.rvdw;
        const double rc9 = rc3 * rc3 * rc3;
        if (ic.vdw_modifier == eintmodPOTSHIFT)
        {
            /* Contribution within the cut-off */
            energy.dispersion += -4.0 * M_PI / (3.0 * rc3);
            energy.repulsion += 4.0 * M_PI / (3.0 * rc9);
        }
        /* Contribution beyond the cut-off */
        addCorrectionBeyondCutoff(&energy, &virial, ic.rvdw);
    }
    else
    {
        gmx_fatal(FARGS, "Dispersion correction is not implemented for vdw-type = %s",
                  evdw_names[ic.vdwtype]);
    }

    iParams->enerdiffsix_    = energy.dispersion;
    iParams->enerdifftwelve_ = energy.repulsion;
    /* The 0.5 is due to the Gromacs definition of the virial */
    iParams->virdiffsix_    = 0.5 * virial.dispersion;
    iParams->virdifftwelve_ = 0.5 * virial.repulsion;
}

DispersionCorrection::DispersionCorrection(const gmx_mtop_t&          mtop,
                                           const t_inputrec&          inputrec,
                                           bool                       useBuckingham,
                                           int                        numAtomTypes,
                                           gmx::ArrayRef<const real>  nonbondedForceParameters,
                                           const interaction_const_t& ic,
                                           const char*                tableFileName) :
    eDispCorr_(inputrec.eDispCorr),
    vdwType_(inputrec.vdwtype),
    eFep_(inputrec.efep),
    topParams_(mtop, inputrec, useBuckingham, numAtomTypes, nonbondedForceParameters)
{
    if (eDispCorr_ != edispcNO)
    {
        GMX_RELEASE_ASSERT(tableFileName, "Need a table file name");

        setInteractionParameters(&iParams_, ic, tableFileName);
    }
}

bool DispersionCorrection::correctFullInteraction() const
{
    return (eDispCorr_ == edispcAllEner || eDispCorr_ == edispcAllEnerPres);
}

void DispersionCorrection::print(const gmx::MDLogger& mdlog) const
{
    if (topParams_.avcsix_[0] == 0 && topParams_.avctwelve_[0] == 0)
    {
        GMX_LOG(mdlog.warning)
                .asParagraph()
                .appendText("WARNING: There are no atom pairs for dispersion correction");
    }
    else if (vdwType_ == evdwUSER)
    {
        GMX_LOG(mdlog.warning)
                .asParagraph()
                .appendText("WARNING: using dispersion correction with user tables\n");
    }

    std::string text = gmx::formatString("Long Range LJ corr.: <C6> %10.4e", topParams_.avcsix_[0]);
    if (correctFullInteraction())
    {
        text += gmx::formatString(" <C12> %10.4e", topParams_.avctwelve_[0]);
    }
    GMX_LOG(mdlog.info).appendText(text);
}

void DispersionCorrection::setParameters(const interaction_const_t& ic)
{
    if (eDispCorr_ != edispcNO)
    {
        setInteractionParameters(&iParams_, ic, nullptr);
    }
}

DispersionCorrection::Correction DispersionCorrection::calculate(const matrix box, const real lambda) const

{
    Correction corr;

    if (eDispCorr_ == edispcNO)
    {
        return corr;
    }

    const bool bCorrAll  = correctFullInteraction();
    const bool bCorrPres = (eDispCorr_ == edispcEnerPres || eDispCorr_ == edispcAllEnerPres);

    const real invvol  = 1 / det(box);
    const real density = topParams_.numAtomsForDensity_ * invvol;
    const real numCorr = topParams_.numCorrections_;

    real avcsix;
    real avctwelve;
    if (eFep_ == efepNO)
    {
        avcsix    = topParams_.avcsix_[0];
        avctwelve = topParams_.avctwelve_[0];
    }
    else
    {
        avcsix    = (1 - lambda) * topParams_.avcsix_[0] + lambda * topParams_.avcsix_[1];
        avctwelve = (1 - lambda) * topParams_.avctwelve_[0] + lambda * topParams_.avctwelve_[1];
    }

    const real enerdiff = numCorr * (density * iParams_.enerdiffsix_ - iParams_.enershiftsix_);
    corr.energy += avcsix * enerdiff;
    real dvdlambda = 0;
    if (eFep_ != efepNO)
    {
        dvdlambda += (topParams_.avcsix_[1] - topParams_.avcsix_[0]) * enerdiff;
    }
    if (bCorrAll)
    {
        const real enerdiff = numCorr * (density * iParams_.enerdifftwelve_ - iParams_.enershifttwelve_);
        corr.energy += avctwelve * enerdiff;
        if (eFep_ != efepNO)
        {
            dvdlambda += (topParams_.avctwelve_[1] - topParams_.avctwelve_[0]) * enerdiff;
        }
    }

    if (bCorrPres)
    {
        corr.virial = numCorr * density * avcsix * iParams_.virdiffsix_ / 3.0;
        if (eDispCorr_ == edispcAllEnerPres)
        {
            corr.virial += numCorr * density * avctwelve * iParams_.virdifftwelve_ / 3.0;
        }
        /* The factor 2 is because of the Gromacs virial definition */
        corr.pressure = -2.0 * invvol * corr.virial * PRESFAC;
    }

    if (eFep_ != efepNO)
    {
        corr.dvdl += dvdlambda;
    }

    return corr;
}