Rename and add doxygen to reaction field params in interaction_const
authorJoe Jordan <ejjordan12@gmail.com>
Fri, 19 Feb 2021 07:01:02 +0000 (07:01 +0000)
committerPaul Bauer <paul.bauer.q@gmail.com>
Fri, 19 Feb 2021 07:01:02 +0000 (07:01 +0000)
This will make further changes to interaction_const RF params easier
to review.

15 files changed:
api/nblib/gmxsetup.cpp
src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp
src/gromacs/mdlib/calc_verletbuf.cpp
src/gromacs/mdlib/forcerec.cpp
src/gromacs/mdrun/tpi.cpp
src/gromacs/mdtypes/interaction_const.h
src/gromacs/nbnxm/benchmark/bench_setup.cpp
src/gromacs/nbnxm/cuda/nbnxm_cuda_kernel.cuh
src/gromacs/nbnxm/kernels_reference/kernel_gpu_ref.cpp
src/gromacs/nbnxm/kernels_reference/kernel_ref_inner.h
src/gromacs/nbnxm/kernels_reference/kernel_ref_outer.h
src/gromacs/nbnxm/kernels_simd_2xmm/kernel_outer.h
src/gromacs/nbnxm/kernels_simd_4xm/kernel_outer.h
src/gromacs/nbnxm/nbnxm_gpu_data_mgmt.cpp
src/gromacs/tables/forcetable.cpp

index 362105448a080b7a6a930756b61837ba54ff8bd7..b314ad7c0fad76e86f620a05b198ae1faeade9f0 100644 (file)
@@ -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))
     {
index c377221472d1afe5b4af9b89a2e137c73934ac22..fffb72a14cac2e310a23036623b7a48c2bf98c15 100644 (file)
@@ -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;
index cad2c62caac6c75d5801a1ef76e01be2c78fe032..3ff12f3342104ba246a099bce52b597f171c90c0 100644 (file)
@@ -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);
             }
         }
index bf8bb9be94dea60a6c793f2a1132d91685fe091a..cdc0b418bfb6269afc9cc399767b2396b4bc91ce 100644 (file)
@@ -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))
         {
index 7a41750ed420600fc377c25a3d023e9af4d710d8..2adb328745531b784b89204690e0408d3f9b4e64 100644 (file)
@@ -144,13 +144,13 @@ static real reactionFieldExclusionCorrection(gmx::ArrayRef<const gmx::RVec> 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);
         }
     }
 
index 9114297bab5e789f37299e823228083122b0e931..37812ee2b3a7e3064e42a83402eedda69fab142e 100644 (file)
@@ -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<EwaldCorrectionTables> coulombEwaldTables;
index 896ac888ab167be290e8ba1639e7ef52b69905cf..78c333b7b27ab5cf8a7159078d0c9e47420848bd 100644 (file)
@@ -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))
     {
index 94877a83c8533c00977ca734176e2aa2efec159a..0ff57b25b120cff746e1d06d00a2a1177eb729b4 100644 (file)
@@ -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 */
index cd0ee7f70ebc5ed8269f2c9ec2646de3d4a0a4e0..6bbef2c32f2ab688732b338f9bdf69a7bfad0fbf 100644 (file)
@@ -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
index 2e11df0728a765e99ea68adeaa3a65a290020d8d..1de0703a89d408f165cd58116d951b4713fd6311 100644 (file)
             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
index df34dbd49003579d80f4976d12a8be0e9f003a75..1203b8e4ba39ba0c8cc4bc8ab7d31e57f36b4d60 100644 (file)
@@ -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
index 4b918be9ca82d7f26f76f983ed859f976539ec55..f45fc6f84382931338b68a262d744c4ac0c36544 100644 (file)
 
 #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
index 187e01440029aef8af1033a506721ee2ee73d87d..ea0ec1fbace70c335ba504544b0be603e35b95d6 100644 (file)
 
 #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
index f811eee3e5ed0a23e952e8908fb1247882532085..344e69d54ed887a8813f44c00b0c555e7684147e 100644 (file)
@@ -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;
index ebd24fe6bb378179217615a50f44d9890842f752..3a8e1694833f7e391aa60c625ad4320910c3f0f4 100644 (file)
@@ -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;