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.
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.
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);
}
}
-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<false>(
scLambdasOrAlphasDiffer, vdwInteractionTypeIsEwald, elecInteractionTypeIsEwald,
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);
}
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;
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++)
* 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!*/
}
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 */
{
/* 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;
/* 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<int>(rtab);
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)
{
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
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
{
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<interaction_const_t::SoftCoreParameters>(*ir->fepvals);
+ }
+
*interaction_const = ic;
}
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)
group.cpp
iforceprovider.cpp
inputrec.cpp
+ interaction_const.cpp
md_enums.cpp
observableshistory.cpp
state.cpp)
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_t> cginfo_mb;
--- /dev/null
+/*
+ * 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 <cstdio>
+
+#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");
+}
#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
*/
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 */
std::unique_ptr<EwaldCorrectionTables> coulombEwaldTables;
// Van der Waals Ewald correction table
std::unique_ptr<EwaldCorrectionTables> vdwEwaldTables;
+
+ // Free-energy parameters, only present when free-energy calculations are requested
+ std::unique_ptr<SoftCoreParameters> softCoreParameters;
};
#endif
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff