From aed54c2748a5d333f2eb32419a7fd24e7acce2ee Mon Sep 17 00:00:00 2001 From: Joe Jordan Date: Fri, 19 Feb 2021 07:01:02 +0000 Subject: [PATCH] Rename and add doxygen to reaction field params in interaction_const This will make further changes to interaction_const RF params easier to review. --- api/nblib/gmxsetup.cpp | 10 +++++----- .../gmxlib/nonbonded/nb_free_energy.cpp | 4 ++-- src/gromacs/mdlib/calc_verletbuf.cpp | 2 +- src/gromacs/mdlib/forcerec.cpp | 19 ++++++++++++------- src/gromacs/mdrun/tpi.cpp | 4 ++-- src/gromacs/mdtypes/interaction_const.h | 13 ++++++++----- src/gromacs/nbnxm/benchmark/bench_setup.cpp | 8 ++++---- src/gromacs/nbnxm/cuda/nbnxm_cuda_kernel.cuh | 11 ++++++----- .../kernels_reference/kernel_gpu_ref.cpp | 6 +++--- .../kernels_reference/kernel_ref_inner.h | 2 +- .../kernels_reference/kernel_ref_outer.h | 8 ++++---- .../nbnxm/kernels_simd_2xmm/kernel_outer.h | 8 ++++---- .../nbnxm/kernels_simd_4xm/kernel_outer.h | 8 ++++---- src/gromacs/nbnxm/nbnxm_gpu_data_mgmt.cpp | 4 ++-- src/gromacs/tables/forcetable.cpp | 4 ++-- 15 files changed, 60 insertions(+), 51 deletions(-) diff --git a/api/nblib/gmxsetup.cpp b/api/nblib/gmxsetup.cpp index 362105448a..b314ad7c0f 100644 --- a/api/nblib/gmxsetup.cpp +++ b/api/nblib/gmxsetup.cpp @@ -254,8 +254,8 @@ void NbvSetupUtil::setupInteractionConst(const NBKernelOptions& options) // 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) @@ -271,10 +271,10 @@ void NbvSetupUtil::setupInteractionConst(const NBKernelOptions& options) 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)) { diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp index c377221472..fffb72a14c 100644 --- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp +++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp @@ -261,8 +261,8 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist, // 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; diff --git a/src/gromacs/mdlib/calc_verletbuf.cpp b/src/gromacs/mdlib/calc_verletbuf.cpp index cad2c62caa..3ff12f3342 100644 --- a/src/gromacs/mdlib/calc_verletbuf.cpp +++ b/src/gromacs/mdlib/calc_verletbuf.cpp @@ -1004,7 +1004,7 @@ real calcVerletBufferSize(const gmx_mtop_t& mtop, } else { - /* epsilon_rf = infinity */ + /* reactionFieldPermitivity = infinity */ k_rf = 0.5 / gmx::power3(ir.rcoulomb); } } diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index bf8bb9be94..cdc0b418bf 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -903,22 +903,27 @@ static void init_interaction_const(FILE* fp, 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; } } @@ -937,7 +942,7 @@ static void init_interaction_const(FILE* fp, if (ic->eeltype == eelCUT) { - fprintf(fp, ", Coulomb %.e", -ic->c_rf); + fprintf(fp, ", Coulomb %.e", -ic->reactionFieldShift); } else if (EEL_PME(ic->eeltype)) { diff --git a/src/gromacs/mdrun/tpi.cpp b/src/gromacs/mdrun/tpi.cpp index 7a41750ed4..2adb328745 100644 --- a/src/gromacs/mdrun/tpi.cpp +++ b/src/gromacs/mdrun/tpi.cpp @@ -144,13 +144,13 @@ static real reactionFieldExclusionCorrection(gmx::ArrayRef x, 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); } } diff --git a/src/gromacs/mdtypes/interaction_const.h b/src/gromacs/mdtypes/interaction_const.h index 9114297bab..37812ee2b3 100644 --- a/src/gromacs/mdtypes/interaction_const.h +++ b/src/gromacs/mdtypes/interaction_const.h @@ -2,7 +2,7 @@ * 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. @@ -135,7 +135,7 @@ struct interaction_const_t 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 */ @@ -158,9 +158,12 @@ struct interaction_const_t 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 coulombEwaldTables; diff --git a/src/gromacs/nbnxm/benchmark/bench_setup.cpp b/src/gromacs/nbnxm/benchmark/bench_setup.cpp index 896ac888ab..78c333b7b2 100644 --- a/src/gromacs/nbnxm/benchmark/bench_setup.cpp +++ b/src/gromacs/nbnxm/benchmark/bench_setup.cpp @@ -1,7 +1,7 @@ /* * 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. @@ -154,10 +154,10 @@ static interaction_const_t setupInteractionConst(const KernelBenchOptions& optio 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)) { diff --git a/src/gromacs/nbnxm/cuda/nbnxm_cuda_kernel.cuh b/src/gromacs/nbnxm/cuda/nbnxm_cuda_kernel.cuh index 94877a83c8..0ff57b25b1 100644 --- a/src/gromacs/nbnxm/cuda/nbnxm_cuda_kernel.cuh +++ b/src/gromacs/nbnxm/cuda/nbnxm_cuda_kernel.cuh @@ -2,7 +2,7 @@ * 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. @@ -204,7 +204,7 @@ __launch_bounds__(THREADS_PER_BLOCK) 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; @@ -357,7 +357,7 @@ __launch_bounds__(THREADS_PER_BLOCK) /* 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 @@ -590,10 +590,11 @@ __launch_bounds__(THREADS_PER_BLOCK) # 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 */ diff --git a/src/gromacs/nbnxm/kernels_reference/kernel_gpu_ref.cpp b/src/gromacs/nbnxm/kernels_reference/kernel_gpu_ref.cpp index cd0ee7f70e..6bbef2c32f 100644 --- a/src/gromacs/nbnxm/kernels_reference/kernel_gpu_ref.cpp +++ b/src/gromacs/nbnxm/kernels_reference/kernel_gpu_ref.cpp @@ -135,7 +135,7 @@ void nbnxn_kernel_gpu_ref(const NbnxnPairlistGpu* nbl, } if (!bEwald) { - vctot *= -facel * 0.5 * iconst->c_rf; + vctot *= -facel * 0.5 * iconst->reactionFieldShift; } else { @@ -230,11 +230,11 @@ void nbnxn_kernel_gpu_ref(const NbnxnPairlistGpu* nbl, 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 diff --git a/src/gromacs/nbnxm/kernels_reference/kernel_ref_inner.h b/src/gromacs/nbnxm/kernels_reference/kernel_ref_inner.h index 2e11df0728..1de0703a89 100644 --- a/src/gromacs/nbnxm/kernels_reference/kernel_ref_inner.h +++ b/src/gromacs/nbnxm/kernels_reference/kernel_ref_inner.h @@ -251,7 +251,7 @@ 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 diff --git a/src/gromacs/nbnxm/kernels_reference/kernel_ref_outer.h b/src/gromacs/nbnxm/kernels_reference/kernel_ref_outer.h index df34dbd490..1203b8e4ba 100644 --- a/src/gromacs/nbnxm/kernels_reference/kernel_ref_outer.h +++ b/src/gromacs/nbnxm/kernels_reference/kernel_ref_outer.h @@ -140,10 +140,10 @@ void #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 @@ -238,7 +238,7 @@ void 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 diff --git a/src/gromacs/nbnxm/kernels_simd_2xmm/kernel_outer.h b/src/gromacs/nbnxm/kernels_simd_2xmm/kernel_outer.h index 4b918be9ca..f45fc6f843 100644 --- a/src/gromacs/nbnxm/kernels_simd_2xmm/kernel_outer.h +++ b/src/gromacs/nbnxm/kernels_simd_2xmm/kernel_outer.h @@ -197,10 +197,10 @@ #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 @@ -401,7 +401,7 @@ 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 diff --git a/src/gromacs/nbnxm/kernels_simd_4xm/kernel_outer.h b/src/gromacs/nbnxm/kernels_simd_4xm/kernel_outer.h index 187e014400..ea0ec1fbac 100644 --- a/src/gromacs/nbnxm/kernels_simd_4xm/kernel_outer.h +++ b/src/gromacs/nbnxm/kernels_simd_4xm/kernel_outer.h @@ -215,10 +215,10 @@ #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 @@ -425,7 +425,7 @@ 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 diff --git a/src/gromacs/nbnxm/nbnxm_gpu_data_mgmt.cpp b/src/gromacs/nbnxm/nbnxm_gpu_data_mgmt.cpp index f811eee3e5..344e69d54e 100644 --- a/src/gromacs/nbnxm/nbnxm_gpu_data_mgmt.cpp +++ b/src/gromacs/nbnxm/nbnxm_gpu_data_mgmt.cpp @@ -152,8 +152,8 @@ void set_cutoff_parameters(NBParamGpu* nbp, const interaction_const_t* ic, const 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; diff --git a/src/gromacs/tables/forcetable.cpp b/src/gromacs/tables/forcetable.cpp index ebd24fe6bb..3a8e169483 100644 --- a/src/gromacs/tables/forcetable.cpp +++ b/src/gromacs/tables/forcetable.cpp @@ -1042,8 +1042,8 @@ static void fill_table(t_tabledata* td, int tp, const interaction_const_t* ic, g 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; -- 2.22.0