// These are the initialized values but we leave them here so that later
// these can become options.
- gmxForceCalculator_->interactionConst_->epsilon_r = 1.0;
- gmxForceCalculator_->interactionConst_->epsilon_rf = 1.0;
+ gmxForceCalculator_->interactionConst_->epsilon_r = 1.0;
+ gmxForceCalculator_->interactionConst_->reactionFieldPermitivity = 1.0;
/* Set the Coulomb energy conversion factor */
if (gmxForceCalculator_->interactionConst_->epsilon_r != 0)
calc_rffac(nullptr,
gmxForceCalculator_->interactionConst_->epsilon_r,
- gmxForceCalculator_->interactionConst_->epsilon_rf,
+ gmxForceCalculator_->interactionConst_->reactionFieldPermitivity,
gmxForceCalculator_->interactionConst_->rcoulomb,
- &gmxForceCalculator_->interactionConst_->k_rf,
- &gmxForceCalculator_->interactionConst_->c_rf);
+ &gmxForceCalculator_->interactionConst_->reactionFieldCoefficient,
+ &gmxForceCalculator_->interactionConst_->reactionFieldShift);
if (EEL_PME_EWALD(gmxForceCalculator_->interactionConst_->eeltype))
{
// Extract data from interaction_const_t
const real facel = ic->epsfac;
const real rCoulomb = ic->rcoulomb;
- const real krf = ic->k_rf;
- const real crf = ic->c_rf;
+ const real krf = ic->reactionFieldCoefficient;
+ const real crf = ic->reactionFieldShift;
const real shLjEwald = ic->sh_lj_ewald;
const real rVdw = ic->rvdw;
const real dispersionShift = ic->dispersion_shift.cpot;
}
else
{
- /* epsilon_rf = infinity */
+ /* reactionFieldPermitivity = infinity */
k_rf = 0.5 / gmx::power3(ir.rcoulomb);
}
}
if (EEL_RF(ic->eeltype))
{
GMX_RELEASE_ASSERT(ic->eeltype != eelGRF_NOTUSED, "GRF is no longer supported");
- ic->epsilon_rf = ir.epsilon_rf;
+ ic->reactionFieldPermitivity = ir.epsilon_rf;
- calc_rffac(fp, ic->epsilon_r, ic->epsilon_rf, ic->rcoulomb, &ic->k_rf, &ic->c_rf);
+ calc_rffac(fp,
+ ic->epsilon_r,
+ ic->reactionFieldPermitivity,
+ ic->rcoulomb,
+ &ic->reactionFieldCoefficient,
+ &ic->reactionFieldShift);
}
else
{
/* For plain cut-off we might use the reaction-field kernels */
- ic->epsilon_rf = ic->epsilon_r;
- ic->k_rf = 0;
+ ic->reactionFieldPermitivity = ic->epsilon_r;
+ ic->reactionFieldCoefficient = 0;
if (ir.coulomb_modifier == eintmodPOTSHIFT)
{
- ic->c_rf = 1 / ic->rcoulomb;
+ ic->reactionFieldShift = 1 / ic->rcoulomb;
}
else
{
- ic->c_rf = 0;
+ ic->reactionFieldShift = 0;
}
}
if (ic->eeltype == eelCUT)
{
- fprintf(fp, ", Coulomb %.e", -ic->c_rf);
+ fprintf(fp, ", Coulomb %.e", -ic->reactionFieldShift);
}
else if (EEL_PME(ic->eeltype))
{
for (int i = beginAtom; i < mdatoms.homenr; i++)
{
const real qi = mdatoms.chargeA[i];
- energy -= 0.5 * qi * qi * ic.c_rf;
+ energy -= 0.5 * qi * qi * ic.reactionFieldShift;
for (int j = i + 1; j < mdatoms.homenr; j++)
{
const real qj = mdatoms.chargeA[j];
const real rsq = distance2(x[i], x[j]);
- energy += qi * qj * (ic.k_rf * rsq - ic.c_rf);
+ energy += qi * qj * (ic.reactionFieldCoefficient * rsq - ic.reactionFieldShift);
}
}
* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 2012,2013,2014,2015,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, 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.
struct shift_consts_t dispersion_shift = { 0, 0, 0 };
struct shift_consts_t repulsion_shift = { 0, 0, 0 };
struct switch_consts_t vdw_switch = { 0, 0, 0 };
- gmx_bool useBuckingham = false;
+ bool useBuckingham = false;
real buckinghamBMax = 0;
/* type of electrostatics */
real epsfac = 1;
/* Constants for reaction-field or plain cut-off */
- real epsilon_rf = 1;
- real k_rf = 0;
- real c_rf = 0;
+ //! Dielectric constant for reaction field beyond the cutoff distance
+ real reactionFieldPermitivity = 1;
+ //! Coefficient for reaction field; scales relation between epsilon_r and reactionFieldPermitivity
+ real reactionFieldCoefficient = 0;
+ //! Constant shift to reaction field Coulomb interaction to make potential an integral of force
+ real reactionFieldShift = 0;
// Coulomb Ewald correction table
std::unique_ptr<EwaldCorrectionTables> coulombEwaldTables;
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020,2021, 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.
ic.coulomb_modifier = eintmodPOTSHIFT;
ic.rcoulomb = options.pairlistCutoff;
- // Reaction-field with epsilon_rf=inf
+ // Reaction-field with reactionFieldPermitivity=inf
// TODO: Replace by calc_rffac() after refactoring that
- ic.k_rf = 0.5 * std::pow(ic.rcoulomb, -3);
- ic.c_rf = 1 / ic.rcoulomb + ic.k_rf * ic.rcoulomb * ic.rcoulomb;
+ ic.reactionFieldCoefficient = 0.5 * std::pow(ic.rcoulomb, -3);
+ ic.reactionFieldShift = 1 / ic.rcoulomb + ic.reactionFieldCoefficient * ic.rcoulomb * ic.rcoulomb;
if (EEL_PME_EWALD(ic.eeltype))
{
* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
- * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2017,2018,2019,2020,2021, 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.
float beta = nbparam.ewald_beta;
float ewald_shift = nbparam.sh_ewald;
# else
- float c_rf = nbparam.c_rf;
+ float reactionFieldShift = nbparam.c_rf;
# endif /* EL_EWALD_ANY */
float* e_lj = atdat.e_lj;
float* e_el = atdat.e_el;
/* Correct for epsfac^2 due to adding qi^2 */
E_el /= nbparam.epsfac * c_clSize * NTHREAD_Z;
# if defined EL_RF || defined EL_CUTOFF
- E_el *= -0.5f * c_rf;
+ E_el *= -0.5f * reactionFieldShift;
# else
E_el *= -beta * M_FLOAT_1_SQRTPI; /* last factor 1/sqrt(pi) */
# endif
# ifdef CALC_ENERGIES
# ifdef EL_CUTOFF
- E_el += qi * qj_f * (int_bit * inv_r - c_rf);
+ E_el += qi * qj_f * (int_bit * inv_r - reactionFieldShift);
# endif
# ifdef EL_RF
- E_el += qi * qj_f * (int_bit * inv_r + 0.5f * two_k_rf * r2 - c_rf);
+ E_el += qi * qj_f
+ * (int_bit * inv_r + 0.5f * two_k_rf * r2 - reactionFieldShift);
# endif
# ifdef EL_EWALD_ANY
/* 1.0f - erff is faster than erfcf */
}
if (!bEwald)
{
- vctot *= -facel * 0.5 * iconst->c_rf;
+ vctot *= -facel * 0.5 * iconst->reactionFieldShift;
}
else
{
if (!bEwald)
{
/* Reaction-field */
- const real krsq = iconst->k_rf * rsq;
+ const real krsq = iconst->reactionFieldCoefficient * rsq;
fscal = qq * (int_bit * rinv - 2 * krsq) * rinvsq;
if (stepWork.computeEnergy)
{
- vcoul = qq * (int_bit * rinv + krsq - iconst->c_rf);
+ vcoul = qq * (int_bit * rinv + krsq - iconst->reactionFieldShift);
}
}
else
real fcoul = qq * (interact * rinv * rinvsq - k_rf2);
/* 4 flops for RF force */
# ifdef CALC_ENERGIES
- real vcoul = qq * (interact * rinv + k_rf * rsq - c_rf);
+ real vcoul = qq * (interact * rinv + reactionFieldCoefficient * rsq - reactionFieldShift);
/* 4 flops for RF energy */
# endif
# endif
#endif
#ifdef CALC_COUL_RF
- const real k_rf2 = 2 * ic->k_rf;
+ const real k_rf2 = 2 * ic->reactionFieldCoefficient;
# ifdef CALC_ENERGIES
- const real k_rf = ic->k_rf;
- const real c_rf = ic->c_rf;
+ const real reactionFieldCoefficient = ic->reactionFieldCoefficient;
+ const real reactionFieldShift = ic->reactionFieldShift;
# endif
#endif
#ifdef CALC_COUL_TAB
if (do_self)
{
# ifdef CALC_COUL_RF
- const real Vc_sub_self = 0.5 * c_rf;
+ const real Vc_sub_self = 0.5 * reactionFieldShift;
# endif
# ifdef CALC_COUL_TAB
# if GMX_DOUBLE
#ifdef CALC_COUL_RF
/* Reaction-field constants */
- mrc_3_S = SimdReal(-2 * ic->k_rf);
+ mrc_3_S = SimdReal(-2 * ic->reactionFieldCoefficient);
# ifdef CALC_ENERGIES
- hrc_3_S = SimdReal(ic->k_rf);
- moh_rc_S = SimdReal(-ic->c_rf);
+ hrc_3_S = SimdReal(ic->reactionFieldCoefficient);
+ moh_rc_S = SimdReal(-ic->reactionFieldShift);
# endif
#endif
if (do_coul)
{
# ifdef CALC_COUL_RF
- const real Vc_sub_self = 0.5 * ic->c_rf;
+ const real Vc_sub_self = 0.5 * ic->reactionFieldShift;
# endif
# ifdef CALC_COUL_TAB
# ifdef TAB_FDV0
#ifdef CALC_COUL_RF
/* Reaction-field constants */
- mrc_3_S = SimdReal(-2 * ic->k_rf);
+ mrc_3_S = SimdReal(-2 * ic->reactionFieldCoefficient);
# ifdef CALC_ENERGIES
- hrc_3_S = SimdReal(ic->k_rf);
- moh_rc_S = SimdReal(-ic->c_rf);
+ hrc_3_S = SimdReal(ic->reactionFieldCoefficient);
+ moh_rc_S = SimdReal(-ic->reactionFieldShift);
# endif
#endif
if (do_coul)
{
# ifdef CALC_COUL_RF
- const real Vc_sub_self = 0.5 * ic->c_rf;
+ const real Vc_sub_self = 0.5 * ic->reactionFieldShift;
# endif
# ifdef CALC_COUL_TAB
# ifdef TAB_FDV0
nbp->ewald_beta = ic->ewaldcoeff_q;
nbp->sh_ewald = ic->sh_ewald;
nbp->epsfac = ic->epsfac;
- nbp->two_k_rf = 2.0 * ic->k_rf;
- nbp->c_rf = ic->c_rf;
+ nbp->two_k_rf = 2.0 * ic->reactionFieldCoefficient;
+ nbp->c_rf = ic->reactionFieldShift;
nbp->rvdw_sq = ic->rvdw * ic->rvdw;
nbp->rcoulomb_sq = ic->rcoulomb * ic->rcoulomb;
nbp->rlistOuter_sq = listParams.rlistOuter * listParams.rlistOuter;
break;
case etabRF:
case etabRF_ZERO:
- Vtab = 1.0 / r + ic->k_rf * r2 - ic->c_rf;
- Ftab = 1.0 / r2 - 2 * ic->k_rf * r;
+ Vtab = 1.0 / r + ic->reactionFieldCoefficient * r2 - ic->reactionFieldShift;
+ Ftab = 1.0 / r2 - 2 * ic->reactionFieldCoefficient * r;
if (tp == etabRF_ZERO && r >= rc)
{
Vtab = 0;