From: Berk Hess Date: Mon, 22 Jun 2020 10:56:00 +0000 (+0000) Subject: Move soft-core parameters to interaction_const_t X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=ea2b13ae71c57af94f15d3f30ba11f093928017c;p=alexxy%2Fgromacs.git Move soft-core parameters to interaction_const_t The soft-core free-energy parameters for use in the kernels have been put in a separate struct called SoftCoreParameters and have been moved from t_forcerec to interaction_const_t. --- diff --git a/docs/release-notes/2021/major/removed-functionality.rst b/docs/release-notes/2021/major/removed-functionality.rst index 4b1b189f6c..d631d1b6f9 100644 --- a/docs/release-notes/2021/major/removed-functionality.rst +++ b/docs/release-notes/2021/major/removed-functionality.rst @@ -7,3 +7,7 @@ Removed functionality Also, please use the syntax :issue:`number` to reference issues on GitLab, without the a space between the colon and number! +Removed GMX_SCSIGMA_MIN environment variable +"""""""""""""""""""""""""""""""""""""""""""" + +This was used to reproduce free-energy soft-core behavior of GROMACS versions before 4.5. diff --git a/docs/user-guide/environment-variables.rst b/docs/user-guide/environment-variables.rst index 4a19a2e208..8dbb431121 100644 --- a/docs/user-guide/environment-variables.rst +++ b/docs/user-guide/environment-variables.rst @@ -352,11 +352,6 @@ Performance and Run Control require the use of tabulated Coulombic and van der Waals interactions. -``GMX_SCSIGMA_MIN`` - the minimum value for soft-core sigma. **Note** that this value is set - using the :mdp:`sc-sigma` keyword in the :ref:`mdp` file, but this environment variable can be used - to reproduce pre-4.5 behavior with respect to this parameter. - ``GMX_TPIC_MASSES`` should contain multiple masses used for test particle insertion into a cavity. The center of mass of the last atoms is used for insertion into the cavity. diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp index e23eef5bac..53f3dbf9d7 100644 --- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp +++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp @@ -246,11 +246,12 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, const real lambda_coul = kernel_data->lambda[efptCOUL]; const real lambda_vdw = kernel_data->lambda[efptVDW]; real* dvdl = kernel_data->dvdl; - const real alpha_coul = fr->sc_alphacoul; - const real alpha_vdw = fr->sc_alphavdw; - const real lam_power = fr->sc_power; - const real sigma6_def = fr->sc_sigma6_def; - const real sigma6_min = fr->sc_sigma6_min; + const auto& scParams = *ic->softCoreParameters; + const real alpha_coul = scParams.alphaCoulomb; + const real alpha_vdw = scParams.alphaVdw; + const real lam_power = scParams.lambdaPower; + const real sigma6_def = scParams.sigma6WithInvalidSigma; + const real sigma6_min = scParams.sigma6Minimum; const bool doForces = ((kernel_data->flags & GMX_NONBONDED_DO_FORCE) != 0); const bool doShiftForces = ((kernel_data->flags & GMX_NONBONDED_DO_SHIFTFORCE) != 0); const bool doPotential = ((kernel_data->flags & GMX_NONBONDED_DO_POTENTIAL) != 0); @@ -938,14 +939,14 @@ static KernelFunction dispatchKernelOnScLambdasOrAlphasDifference(const bool scL } } -static KernelFunction dispatchKernel(const bool scLambdasOrAlphasDiffer, - const bool vdwInteractionTypeIsEwald, - const bool elecInteractionTypeIsEwald, - const bool vdwModifierIsPotSwitch, - const bool useSimd, - const t_forcerec* fr) +static KernelFunction dispatchKernel(const bool scLambdasOrAlphasDiffer, + const bool vdwInteractionTypeIsEwald, + const bool elecInteractionTypeIsEwald, + const bool vdwModifierIsPotSwitch, + const bool useSimd, + const interaction_const_t& ic) { - if (fr->sc_alphacoul == 0 && fr->sc_alphavdw == 0) + if (ic.softCoreParameters->alphaCoulomb == 0 && ic.softCoreParameters->alphaVdw == 0) { return (dispatchKernelOnScLambdasOrAlphasDifference( scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald, @@ -968,33 +969,33 @@ void gmx_nb_free_energy_kernel(const t_nblist* nlist, nb_kernel_data_t* kernel_data, t_nrnb* nrnb) { - GMX_ASSERT(EEL_PME_EWALD(fr->ic->eeltype) || fr->ic->eeltype == eelCUT || EEL_RF(fr->ic->eeltype), + const interaction_const_t& ic = *fr->ic; + GMX_ASSERT(EEL_PME_EWALD(ic.eeltype) || ic.eeltype == eelCUT || EEL_RF(ic.eeltype), "Unsupported eeltype with free energy"); + GMX_ASSERT(ic.softCoreParameters, "We need soft-core parameters"); - const bool vdwInteractionTypeIsEwald = (EVDW_PME(fr->ic->vdwtype)); - const bool elecInteractionTypeIsEwald = (EEL_PME_EWALD(fr->ic->eeltype)); - const bool vdwModifierIsPotSwitch = (fr->ic->vdw_modifier == eintmodPOTSWITCH); - bool scLambdasOrAlphasDiffer = true; - const bool useSimd = fr->use_simd_kernels; + const auto& scParams = *ic.softCoreParameters; + const bool vdwInteractionTypeIsEwald = (EVDW_PME(ic.vdwtype)); + const bool elecInteractionTypeIsEwald = (EEL_PME_EWALD(ic.eeltype)); + const bool vdwModifierIsPotSwitch = (ic.vdw_modifier == eintmodPOTSWITCH); + bool scLambdasOrAlphasDiffer = true; + const bool useSimd = fr->use_simd_kernels; - if (fr->sc_alphacoul == 0 && fr->sc_alphavdw == 0) + if (scParams.alphaCoulomb == 0 && scParams.alphaVdw == 0) { scLambdasOrAlphasDiffer = false; } - else if (fr->sc_r_power == 6.0_real) + else { - if (kernel_data->lambda[efptCOUL] == kernel_data->lambda[efptVDW] && fr->sc_alphacoul == fr->sc_alphavdw) + if (kernel_data->lambda[efptCOUL] == kernel_data->lambda[efptVDW] + && scParams.alphaCoulomb == scParams.alphaVdw) { scLambdasOrAlphasDiffer = false; } } - else - { - GMX_RELEASE_ASSERT(false, "Unsupported soft-core r-power"); - } KernelFunction kernelFunc; kernelFunc = dispatchKernel(scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, - elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, useSimd, fr); + elecInteractionTypeIsEwald, vdwModifierIsPotSwitch, useSimd, ic); kernelFunc(nlist, xx, ff, fr, mdatoms, kernel_data, nrnb); } diff --git a/src/gromacs/listed_forces/pairs.cpp b/src/gromacs/listed_forces/pairs.cpp index 7945d689f0..e41bac0f32 100644 --- a/src/gromacs/listed_forces/pairs.cpp +++ b/src/gromacs/listed_forces/pairs.cpp @@ -154,47 +154,44 @@ static real evaluate_single(real r2, return fscal; } +static inline real sixthRoot(const real r) +{ + return gmx::invsqrt(std::cbrt(r)); +} + /*! \brief Compute the energy and force for a single pair interaction under FEP */ -static real free_energy_evaluate_single(real r2, - real sc_r_power, - real alpha_coul, - real alpha_vdw, - real tabscale, - const real* vftab, - real tableStride, - real qqA, - real c6A, - real c12A, - real qqB, - real c6B, - real c12B, - const real LFC[2], - const real LFV[2], - const real DLF[2], - const real lfac_coul[2], - const real lfac_vdw[2], - const real dlfac_coul[2], - const real dlfac_vdw[2], - real sigma6_def, - real sigma6_min, - real sigma2_def, - real sigma2_min, - real* velectot, - real* vvdwtot, - real* dvdl) +static real free_energy_evaluate_single(real r2, + const interaction_const_t::SoftCoreParameters& scParams, + real tabscale, + const real* vftab, + real tableStride, + real qqA, + real c6A, + real c12A, + real qqB, + real c6B, + real c12B, + const real LFC[2], + const real LFV[2], + const real DLF[2], + const real lfac_coul[2], + const real lfac_vdw[2], + const real dlfac_coul[2], + const real dlfac_vdw[2], + real* velectot, + real* vvdwtot, + real* dvdl) { - real rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal; - real qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2]; + real rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal; + real qq[2], c6[2], c12[2], sigma6[2], sigma_pow[2]; real alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw; real rpinv, r_coul, r_vdw, velecsum, vvdwsum; real fscal_vdw[2], fscal_elec[2]; real velec[2], vvdw[2]; int i, ntab; - const real half = 0.5; - const real minusOne = -1.0; - const real one = 1.0; - const real two = 2.0; - const real six = 6.0; + const real half = 0.5_real; + const real one = 1.0_real; + const real two = 2.0_real; qq[0] = qqA; qq[1] = qqB; @@ -203,16 +200,8 @@ static real free_energy_evaluate_single(real r2, c12[0] = c12A; c12[1] = c12B; - if (sc_r_power == six) - { - rpm2 = r2 * r2; /* r4 */ - rp = rpm2 * r2; /* r6 */ - } - else - { - rp = std::pow(r2, half * sc_r_power); /* not currently supported as input, but can handle it */ - rpm2 = rp / r2; - } + const real rpm2 = r2 * r2; /* r4 */ + const real rp = rpm2 * r2; /* r6 */ /* Loop over state A(0) and B(1) */ for (i = 0; i < 2; i++) @@ -223,28 +212,16 @@ static real free_energy_evaluate_single(real r2, * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5. */ sigma6[i] = half * c12[i] / c6[i]; - sigma2[i] = std::cbrt(half * c12[i] / c6[i]); - /* should be able to get rid of this ^^^ internal pow call eventually. Will require agreement on - what data to store externally. Can't be fixed without larger scale changes, so not 5.0 */ - if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */ + if (sigma6[i] < scParams.sigma6Minimum) /* for disappearing coul and vdw with soft core at the same time */ { - sigma6[i] = sigma6_min; - sigma2[i] = sigma2_min; + sigma6[i] = scParams.sigma6Minimum; } } else { - sigma6[i] = sigma6_def; - sigma2[i] = sigma2_def; - } - if (sc_r_power == six) - { - sigma_pow[i] = sigma6[i]; - } - else - { /* not really supported as input, but in here for testing the general case*/ - sigma_pow[i] = std::pow(sigma2[i], sc_r_power / 2); + sigma6[i] = scParams.sigma6WithInvalidSigma; } + sigma_pow[i] = sigma6[i]; } /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/ @@ -255,8 +232,8 @@ static real free_energy_evaluate_single(real r2, } else { - alpha_vdw_eff = alpha_vdw; - alpha_coul_eff = alpha_coul; + alpha_vdw_eff = scParams.alphaVdw; + alpha_coul_eff = scParams.alphaCoulomb; } /* Loop over A and B states again */ @@ -272,7 +249,7 @@ static real free_energy_evaluate_single(real r2, { /* Coulomb */ rpinv = one / (alpha_coul_eff * lfac_coul[i] * sigma_pow[i] + rp); - r_coul = std::pow(rpinv, minusOne / sc_r_power); + r_coul = sixthRoot(rpinv); /* Electrostatics table lookup data */ rtab = r_coul * tabscale; @@ -293,7 +270,7 @@ static real free_energy_evaluate_single(real r2, /* Vdw */ rpinv = one / (alpha_vdw_eff * lfac_vdw[i] * sigma_pow[i] + rp); - r_vdw = std::pow(rpinv, minusOne / sc_r_power); + r_vdw = sixthRoot(rpinv); /* Vdw table lookup data */ rtab = r_vdw * tabscale; ntab = static_cast(rtab); @@ -380,9 +357,10 @@ static real do_pairs_general(int ftype, real* energygrp_vdw; static gmx_bool warned_rlimit = FALSE; /* Free energy stuff */ - gmx_bool bFreeEnergy; - real LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2]; - real qqB, c6B, c12B, sigma2_def, sigma2_min; + gmx_bool bFreeEnergy; + real LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2]; + real qqB, c6B, c12B; + const real oneSixth = 1.0_real / 6.0_real; switch (ftype) { @@ -413,24 +391,19 @@ static real do_pairs_general(int ftype, DLF[0] = -1; DLF[1] = 1; - /* precalculate */ - sigma2_def = std::cbrt(fr->sc_sigma6_def); - sigma2_min = std::cbrt(fr->sc_sigma6_min); + GMX_ASSERT(fr->ic->softCoreParameters, "We need soft-core parameters"); + const auto& scParams = *fr->ic->softCoreParameters; for (i = 0; i < 2; i++) { - lfac_coul[i] = (fr->sc_power == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i])); - dlfac_coul[i] = - DLF[i] * fr->sc_power / fr->sc_r_power * (fr->sc_power == 2 ? (1 - LFC[i]) : 1); - lfac_vdw[i] = (fr->sc_power == 2 ? (1 - LFV[i]) * (1 - LFV[i]) : (1 - LFV[i])); - dlfac_vdw[i] = - DLF[i] * fr->sc_power / fr->sc_r_power * (fr->sc_power == 2 ? (1 - LFV[i]) : 1); + lfac_coul[i] = (scParams.lambdaPower == 2 ? (1 - LFC[i]) * (1 - LFC[i]) : (1 - LFC[i])); + dlfac_coul[i] = DLF[i] * scParams.lambdaPower * oneSixth + * (scParams.lambdaPower == 2 ? (1 - LFC[i]) : 1); + lfac_vdw[i] = (scParams.lambdaPower == 2 ? (1 - LFV[i]) * (1 - LFV[i]) : (1 - LFV[i])); + dlfac_vdw[i] = DLF[i] * scParams.lambdaPower * oneSixth + * (scParams.lambdaPower == 2 ? (1 - LFV[i]) : 1); } } - else - { - sigma2_min = sigma2_def = 0; - } /* TODO This code depends on the logic in tables.c that constructs the table layout, which should be made explicit in future @@ -521,10 +494,9 @@ static real do_pairs_general(int ftype, c12B = iparams[itype].lj14.c12B * 12.0; fscal = free_energy_evaluate_single( - r2, fr->sc_r_power, fr->sc_alphacoul, fr->sc_alphavdw, fr->pairsTable->scale, - fr->pairsTable->data.data(), fr->pairsTable->stride, qq, c6, c12, qqB, c6B, c12B, - LFC, LFV, DLF, lfac_coul, lfac_vdw, dlfac_coul, dlfac_vdw, fr->sc_sigma6_def, - fr->sc_sigma6_min, sigma2_def, sigma2_min, &velec, &vvdw, dvdl); + r2, *fr->ic->softCoreParameters, fr->pairsTable->scale, fr->pairsTable->data.data(), + fr->pairsTable->stride, qq, c6, c12, qqB, c6B, c12B, LFC, LFV, DLF, lfac_coul, + lfac_vdw, dlfac_coul, dlfac_vdw, &velec, &vvdw, dvdl); } else { diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index e96755eeb6..dad37af31c 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -937,6 +937,12 @@ static void init_interaction_const(FILE* fp, fprintf(fp, "\n"); } + if (ir->efep != efepNO) + { + GMX_RELEASE_ASSERT(ir->fepvals, "ir->fepvals should be set wth free-energy"); + ic->softCoreParameters = std::make_unique(*ir->fepvals); + } + *interaction_const = ic; } @@ -1000,33 +1006,7 @@ void init_forcerec(FILE* fp, fr->fc_stepsize = ir->fc_stepsize; /* Free energy */ - fr->efep = ir->efep; - fr->sc_alphavdw = ir->fepvals->sc_alpha; - if (ir->fepvals->bScCoul) - { - fr->sc_alphacoul = ir->fepvals->sc_alpha; - fr->sc_sigma6_min = gmx::power6(ir->fepvals->sc_sigma_min); - } - else - { - fr->sc_alphacoul = 0; - fr->sc_sigma6_min = 0; /* only needed when bScCoul is on */ - } - fr->sc_power = ir->fepvals->sc_power; - fr->sc_r_power = ir->fepvals->sc_r_power; - fr->sc_sigma6_def = gmx::power6(ir->fepvals->sc_sigma); - - char* env = getenv("GMX_SCSIGMA_MIN"); - if (env != nullptr) - { - double dbl = 0; - sscanf(env, "%20lf", &dbl); - fr->sc_sigma6_min = gmx::power6(dbl); - if (fp) - { - fprintf(fp, "Setting the minimum soft core sigma to %g nm\n", dbl); - } - } + fr->efep = ir->efep; fr->bNonbonded = TRUE; if (getenv("GMX_NO_NONBONDED") != nullptr) diff --git a/src/gromacs/mdtypes/CMakeLists.txt b/src/gromacs/mdtypes/CMakeLists.txt index 9e54972b24..da3cecf437 100644 --- a/src/gromacs/mdtypes/CMakeLists.txt +++ b/src/gromacs/mdtypes/CMakeLists.txt @@ -37,6 +37,7 @@ file(GLOB MDTYPES_SOURCES group.cpp iforceprovider.cpp inputrec.cpp + interaction_const.cpp md_enums.cpp observableshistory.cpp state.cpp) diff --git a/src/gromacs/mdtypes/forcerec.h b/src/gromacs/mdtypes/forcerec.h index 81238d4b32..044f8d332c 100644 --- a/src/gromacs/mdtypes/forcerec.h +++ b/src/gromacs/mdtypes/forcerec.h @@ -229,13 +229,7 @@ struct t_forcerec t_forcetable* pairsTable = nullptr; /* for 1-4 interactions, [pairs] and [pairs_nb] */ /* Free energy */ - int efep = 0; - real sc_alphavdw = 0; - real sc_alphacoul = 0; - int sc_power = 0; - real sc_r_power = 0; - real sc_sigma6_def = 0; - real sc_sigma6_min = 0; + int efep = 0; /* Information about atom properties for the molecule blocks in the system */ std::vector cginfo_mb; diff --git a/src/gromacs/mdtypes/interaction_const.cpp b/src/gromacs/mdtypes/interaction_const.cpp new file mode 100644 index 0000000000..230dca2bff --- /dev/null +++ b/src/gromacs/mdtypes/interaction_const.cpp @@ -0,0 +1,53 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2020, 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 "interaction_const.h" + +#include + +#include "gromacs/math/functions.h" +#include "gromacs/mdtypes/inputrec.h" + +interaction_const_t::SoftCoreParameters::SoftCoreParameters(const t_lambda& fepvals) : + alphaVdw(fepvals.sc_alpha), + alphaCoulomb(fepvals.bScCoul ? fepvals.sc_alpha : 0), + lambdaPower(fepvals.sc_power), + sigma6WithInvalidSigma(gmx::power6(fepvals.sc_sigma)), + sigma6Minimum(fepvals.bScCoul ? gmx::power6(fepvals.sc_sigma_min) : 0) +{ + // This is checked during tpr reading, so we can assert here + GMX_RELEASE_ASSERT(fepvals.sc_r_power == 6.0, "We only support soft-core r-power 6"); +} diff --git a/src/gromacs/mdtypes/interaction_const.h b/src/gromacs/mdtypes/interaction_const.h index e9378fe833..8c4d3f373a 100644 --- a/src/gromacs/mdtypes/interaction_const.h +++ b/src/gromacs/mdtypes/interaction_const.h @@ -44,6 +44,8 @@ #include "gromacs/utility/basedefinitions.h" #include "gromacs/utility/real.h" +struct t_lambda; + /* Used with force switching or a constant potential shift: * rsw = max(r - r_switch, 0) * force/p = r^-(p+1) + c2*rsw^2 + c3*rsw^3 @@ -104,6 +106,28 @@ struct EwaldCorrectionTables */ struct interaction_const_t { + /* This struct contains the soft-core parameters from t_lambda, + * but processed for direct use in the kernels. + */ + struct SoftCoreParameters + { + // Constructor + SoftCoreParameters(const t_lambda& fepvals); + + // Alpha parameter for Van der Waals interactions + real alphaVdw; + // Alpha parameter for Coulomb interactions + real alphaCoulomb; + // Exponent for the dependence of the soft-core on lambda + int lambdaPower; + // Value for sigma^6 for LJ interaction with C6<=0 and/or C12<=0 + real sigma6WithInvalidSigma; + // Minimum value for sigma^6, used when soft-core is applied to Coulomb interactions + real sigma6Minimum; + }; + + // Cut-off scheme, only present for reading and (not) running old tpr files + // which still supported the group cutoff-scheme int cutoff_scheme = ecutsVERLET; /* VdW */ @@ -146,6 +170,9 @@ struct interaction_const_t std::unique_ptr coulombEwaldTables; // Van der Waals Ewald correction table std::unique_ptr vdwEwaldTables; + + // Free-energy parameters, only present when free-energy calculations are requested + std::unique_ptr softCoreParameters; }; #endif