Remove forcerec and inputrec from ewald_LRcorrection

Only pass needed parameters to ewald_LRcorrection, instead of
entire forcerec and inputrec. This required also refactoring
EwaldBoxZScaler to only take needed parameters instead of a
whole inputrec.

Also made some names more verbose.

diff --git a/src/gromacs/ewald/ewald.cpp b/src/gromacs/ewald/ewald.cpp
index b13e3e7c365e51a1673600236c774400464fb15f..4d04a4f512e059beaefc8dcd1fa39a5240b94c67 100644 (file)
--- a/src/gromacs/ewald/ewald.cpp
+++ b/src/gromacs/ewald/ewald.cpp
@@ -163,7 +163,7 @@ real do_ewald(const t_inputrec&              ir,
 
     /* Scale box with Ewald wall factor */
     matrix          scaledBox;
-    EwaldBoxZScaler boxScaler(ir);
+    EwaldBoxZScaler boxScaler(inputrecPbcXY2Walls(&ir), ir.wall_ewald_zfac);
     boxScaler.scaleBox(box, scaledBox);
 
     rvec boxDiag;

diff --git a/src/gromacs/ewald/ewald_utils.h b/src/gromacs/ewald/ewald_utils.h
index 4639d1571ba911c1d0bd36c009c0f6887edf3df6..20db56fcebda68f34a569b6527d74f3f646a6a1e 100644 (file)
--- a/src/gromacs/ewald/ewald_utils.h
+++ b/src/gromacs/ewald/ewald_utils.h
@@ -100,12 +100,12 @@ public:
     EwaldBoxZScaler() = delete;
 
     /*! \brief Constructor that takes the input record to initialize Ewald box scaling appropriately. */
-    EwaldBoxZScaler(const t_inputrec& ir)
+    EwaldBoxZScaler(bool havePbcXY2Walls, real wallEwaldZfac)
     {
-        if (inputrecPbcXY2Walls(&ir))
+        if (havePbcXY2Walls)
         {
             scaleWithWalls_ = true;
-            scalingFactor_  = ir.wall_ewald_zfac;
+            scalingFactor_  = wallEwaldZfac;
         }
         else
         {

diff --git a/src/gromacs/ewald/long_range_correction.cpp b/src/gromacs/ewald/long_range_correction.cpp
index b5791ead4746be2ef7935312c5cbc0e3a8f213ea..b3316285998c53aa9770827ab359afaed0ae033f 100644 (file)
--- a/src/gromacs/ewald/long_range_correction.cpp
+++ b/src/gromacs/ewald/long_range_correction.cpp
@@ -65,18 +65,22 @@
  * undefined when LJ-PME is not active. This works because
  * bHaveChargeOrTypePerturbed handles the control flow. */
 void ewald_LRcorrection(const int                      numAtomsLocal,
-                        const t_commrec*               cr,
+                        const t_commrec*               commrec,
                         int                            numThreads,
                         int                            thread,
-                        const t_forcerec&              fr,
-                        const t_inputrec&              ir,
+                        const real                     epsilonR,
+                        gmx::ArrayRef<const double>    qsum,
+                        EwaldGeometry                  ewaldGeometry,
+                        const real                     epsilonSurface,
+                        bool                           havePbcXY2Walls,
+                        real                           wallEwaldZfac,
                         gmx::ArrayRef<const real>      chargeA,
                         gmx::ArrayRef<const real>      chargeB,
                         bool                           bHaveChargePerturbed,
-                        gmx::ArrayRef<const gmx::RVec> x,
+                        gmx::ArrayRef<const gmx::RVec> coords,
                         const matrix                   box,
                         gmx::ArrayRef<const gmx::RVec> mu_tot,
-                        gmx::ArrayRef<gmx::RVec>       f,
+                        gmx::ArrayRef<gmx::RVec>       forces,
                         real*                          Vcorr_q,
                         real                           lambda_q,
                         real*                          dvdlambda_q)
@@ -95,10 +99,10 @@ void ewald_LRcorrection(const int                      numAtomsLocal,
 
     /* Scale the Ewald unit cell when dimension z is not periodic */
     matrix          scaledBox;
-    EwaldBoxZScaler boxScaler(ir);
+    EwaldBoxZScaler boxScaler(havePbcXY2Walls, wallEwaldZfac);
     boxScaler.scaleBox(box, scaledBox);
 
-    one_4pi_eps = gmx::c_one4PiEps0 / fr.ic->epsilon_r;
+    one_4pi_eps = gmx::c_one4PiEps0 / epsilonR;
     Vexcl_q     = 0;
     dvdl_excl_q = 0;
     Vdipole[0]  = 0;
@@ -117,13 +121,13 @@ void ewald_LRcorrection(const int                      numAtomsLocal,
     dipole_coeff = 0;
 
     real boxVolume = scaledBox[XX][XX] * scaledBox[YY][YY] * scaledBox[ZZ][ZZ];
-    switch (ir.ewald_geometry)
+    switch (ewaldGeometry)
     {
         case EwaldGeometry::ThreeD:
-            if (ir.epsilon_surface != 0)
+            if (epsilonSurface != 0)
             {
-                dipole_coeff = 2 * M_PI * gmx::c_one4PiEps0
-                               / ((2 * ir.epsilon_surface + fr.ic->epsilon_r) * boxVolume);
+                dipole_coeff =
+                        2 * M_PI * gmx::c_one4PiEps0 / ((2 * epsilonSurface + epsilonR) * boxVolume);
                 for (i = 0; (i < DIM); i++)
                 {
                     dipcorrA[i] = 2 * dipole_coeff * mutot[0][i];
@@ -138,19 +142,14 @@ void ewald_LRcorrection(const int                      numAtomsLocal,
             for (int q = 0; q < (bHaveChargePerturbed ? 2 : 1); q++)
             {
                 /* Avoid charge corrections with near-zero net charge */
-                if (fabs(fr.qsum[q]) > 1e-4)
+                if (fabs(qsum[q]) > 1e-4)
                 {
-                    chargecorr[q] = 2 * dipole_coeff * fr.qsum[q];
+                    chargecorr[q] = 2 * dipole_coeff * qsum[q];
                 }
             }
             break;
         default: gmx_incons("Unsupported Ewald geometry");
     }
-    if (debug)
-    {
-        fprintf(debug, "dipcorr = %8.3f  %8.3f  %8.3f\n", dipcorrA[XX], dipcorrA[YY], dipcorrA[ZZ]);
-        fprintf(debug, "mutot   = %8.3f  %8.3f  %8.3f\n", mutot[0][XX], mutot[0][YY], mutot[0][ZZ]);
-    }
     const bool bNeedLongRangeCorrection = (dipole_coeff != 0);
     if (bNeedLongRangeCorrection && !bHaveChargePerturbed)
     {
@@ -158,11 +157,11 @@ void ewald_LRcorrection(const int                      numAtomsLocal,
         {
             for (j = 0; (j < DIM); j++)
             {
-                f[i][j] -= dipcorrA[j] * chargeA[i];
+                forces[i][j] -= dipcorrA[j] * chargeA[i];
             }
             if (chargecorr[0] != 0)
             {
-                f[i][ZZ] += chargecorr[0] * chargeA[i] * x[i][ZZ];
+                forces[i][ZZ] += chargecorr[0] * chargeA[i] * coords[i][ZZ];
             }
         }
     }
@@ -172,11 +171,12 @@ void ewald_LRcorrection(const int                      numAtomsLocal,
         {
             for (j = 0; (j < DIM); j++)
             {
-                f[i][j] -= L1_q * dipcorrA[j] * chargeA[i] + lambda_q * dipcorrB[j] * chargeB[i];
+                forces[i][j] -= L1_q * dipcorrA[j] * chargeA[i] + lambda_q * dipcorrB[j] * chargeB[i];
             }
             if (chargecorr[0] != 0 || chargecorr[1] != 0)
             {
-                f[i][ZZ] += (L1_q * chargecorr[0] * chargeA[i] + lambda_q * chargecorr[1]) * x[i][ZZ];
+                forces[i][ZZ] +=
+                        (L1_q * chargecorr[0] * chargeA[i] + lambda_q * chargecorr[1]) * coords[i][ZZ];
             }
         }
     }
@@ -185,7 +185,7 @@ void ewald_LRcorrection(const int                      numAtomsLocal,
     Vself_q[1] = 0;
 
     /* Global corrections only on master process */
-    if (MASTER(cr) && thread == 0)
+    if (MASTER(commrec) && thread == 0)
     {
         for (q = 0; q < (bHaveChargePerturbed ? 2 : 1); q++)
         {
@@ -195,11 +195,11 @@ void ewald_LRcorrection(const int                      numAtomsLocal,
              */
             if (dipole_coeff != 0)
             {
-                if (ir.ewald_geometry == EwaldGeometry::ThreeD)
+                if (ewaldGeometry == EwaldGeometry::ThreeD)
                 {
                     Vdipole[q] = dipole_coeff * iprod(mutot[q], mutot[q]);
                 }
-                else if (ir.ewald_geometry == EwaldGeometry::ThreeDC)
+                else if (ewaldGeometry == EwaldGeometry::ThreeDC)
                 {
                     Vdipole[q] = dipole_coeff * mutot[q][ZZ] * mutot[q][ZZ];
 
@@ -215,10 +215,10 @@ void ewald_LRcorrection(const int                      numAtomsLocal,
                         real                      sumQZ2 = 0;
                         for (int i = 0; i < numAtomsLocal; i++)
                         {
-                            sumQZ2 += charge[i] * x[i][ZZ] * x[i][ZZ];
+                            sumQZ2 += charge[i] * coords[i][ZZ] * coords[i][ZZ];
                         }
-                        Vdipole[q] -= dipole_coeff * fr.qsum[q]
-                                      * (sumQZ2 + fr.qsum[q] * box[ZZ][ZZ] * box[ZZ][ZZ] / 12);
+                        Vdipole[q] -= dipole_coeff * qsum[q]
+                                      * (sumQZ2 + qsum[q] * box[ZZ][ZZ] * box[ZZ][ZZ] / 12);
                     }
                 }
             }
@@ -233,23 +233,4 @@ void ewald_LRcorrection(const int                      numAtomsLocal,
         *Vcorr_q = L1_q * (Vdipole[0] - Vself_q[0]) + lambda_q * (Vdipole[1] - Vself_q[1]) - Vexcl_q;
         *dvdlambda_q += Vdipole[1] - Vself_q[1] - (Vdipole[0] - Vself_q[0]) - dvdl_excl_q;
     }
-
-    if (debug)
-    {
-        fprintf(debug, "Long Range corrections for Ewald interactions:\n");
-        fprintf(debug,
-                "q2sum = %g, Vself_q=%g\n",
-                L1_q * fr.q2sum[0] + lambda_q * fr.q2sum[1],
-                L1_q * Vself_q[0] + lambda_q * Vself_q[1]);
-        fprintf(debug, "Electrostatic Long Range correction: Vexcl=%g\n", Vexcl_q);
-        if (MASTER(cr) && thread == 0)
-        {
-            if (ir.epsilon_surface > 0 || ir.ewald_geometry == EwaldGeometry::ThreeDC)
-            {
-                fprintf(debug,
-                        "Total dipole correction: Vdipole=%g\n",
-                        L1_q * Vdipole[0] + lambda_q * Vdipole[1]);
-            }
-        }
-    }
 }

diff --git a/src/gromacs/ewald/long_range_correction.h b/src/gromacs/ewald/long_range_correction.h
index f6fb0c5e8cc1d36c18e996a433cdfb573621d88d..b70e3b717d08d6f52aef9ed752bdd93c1f9f43fc 100644 (file)
--- a/src/gromacs/ewald/long_range_correction.h
+++ b/src/gromacs/ewald/long_range_correction.h
@@ -57,6 +57,7 @@
 struct t_commrec;
 struct t_forcerec;
 struct t_inputrec;
+enum class EwaldGeometry : int;
 
 namespace gmx
 {
@@ -69,18 +70,22 @@ class ArrayRef;
  * Calculate correction for electrostatic surface dipole terms.
  */
 void ewald_LRcorrection(int                            numAtomsLocal,
-                        const t_commrec*               cr,
+                        const t_commrec*               commrec,
                         int                            numThreads,
                         int                            thread,
-                        const t_forcerec&              fr,
-                        const t_inputrec&              ir,
+                        real                           epsilonR,
+                        gmx::ArrayRef<const double>    qsum,
+                        EwaldGeometry                  ewaldGeometry,
+                        real                           epsilonSurface,
+                        bool                           havePbcXY2Walls,
+                        real                           wallEwaldZfac,
                         gmx::ArrayRef<const real>      chargeA,
                         gmx::ArrayRef<const real>      chargeB,
                         bool                           bHaveChargePerturbed,
-                        gmx::ArrayRef<const gmx::RVec> x,
+                        gmx::ArrayRef<const gmx::RVec> coords,
                         const matrix                   box,
                         gmx::ArrayRef<const gmx::RVec> mu_tot,
-                        gmx::ArrayRef<gmx::RVec>       f,
+                        gmx::ArrayRef<gmx::RVec>       forces,
                         real*                          Vcorr_q,
                         real                           lambda_q,
                         real*                          dvdlambda_q);

diff --git a/src/gromacs/ewald/pme.cpp b/src/gromacs/ewald/pme.cpp
index 9044c496d21fc7a340ba8b75ca692d29a217b3c2..a3dd8357556c86b88c0f77ba32c2848d0be5b959 100644 (file)
--- a/src/gromacs/ewald/pme.cpp
+++ b/src/gromacs/ewald/pme.cpp
@@ -728,7 +728,7 @@ gmx_pme_t* gmx_pme_init(const t_commrec*     cr,
     // The box requires scaling with nwalls = 2, we store that condition as well
     // as the scaling factor
     delete pme->boxScaler;
-    pme->boxScaler = new EwaldBoxZScaler(*ir);
+    pme->boxScaler = new EwaldBoxZScaler(inputrecPbcXY2Walls(ir), ir->wall_ewald_zfac);
 
     /* If we violate restrictions, generate a fatal error here */
     gmx_pme_check_restrictions(

diff --git a/src/gromacs/ewald/pme_load_balancing.cpp b/src/gromacs/ewald/pme_load_balancing.cpp
index f554e43827c1e4e2eb6a2de614a3bfafdd76460d..52aaeb98621fec1c374449e8d60aa92cf1072448 100644 (file)
--- a/src/gromacs/ewald/pme_load_balancing.cpp
+++ b/src/gromacs/ewald/pme_load_balancing.cpp
@@ -251,7 +251,7 @@ void pme_loadbal_init(pme_load_balancing_t**     pme_lb_p,
     /* Scale box with Ewald wall factor; note that we pmedata->boxScaler
      * can't always usedd as it's not available with separate PME ranks.
      */
-    EwaldBoxZScaler boxScaler(ir);
+    EwaldBoxZScaler boxScaler(inputrecPbcXY2Walls(&ir), ir.wall_ewald_zfac);
     boxScaler.scaleBox(box, pme_lb->box_start);
 
     pme_lb->setup.resize(1);

diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp
index 63f93ae8db332692d91f8316fa314d7d56696a69..64d6de8c58bfc06a776db0ecabe2000de861df0f 100644 (file)
--- a/src/gromacs/gmxpreprocess/grompp.cpp
+++ b/src/gromacs/gmxpreprocess/grompp.cpp
@@ -2343,7 +2343,7 @@ int gmx_grompp(int argc, char* argv[])
     {
         /* Calculate the optimal grid dimensions */
         matrix          scaledBox;
-        EwaldBoxZScaler boxScaler(*ir);
+        EwaldBoxZScaler boxScaler(inputrecPbcXY2Walls(ir), ir->wall_ewald_zfac);
         boxScaler.scaleBox(state.box, scaledBox);
 
         if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0)

diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp
index 8f07ca623113671f55b66c0b7107b417e1bc00e2..11123448175767bc13e72cb3401d81aa8373e994 100644 (file)
--- a/src/gromacs/mdlib/force.cpp
+++ b/src/gromacs/mdlib/force.cpp
@@ -165,8 +165,12 @@ void calculateLongRangeNonbondeds(t_forcerec*                    fr,
                                 cr,
                                 nthreads,
                                 t,
-                                *fr,
-                                ir,
+                                fr->ic->epsilon_r,
+                                fr->qsum,
+                                ir.ewald_geometry,
+                                ir.epsilon_surface,
+                                inputrecPbcXY2Walls(&ir),
+                                ir.wall_ewald_zfac,
                                 md->chargeA ? gmx::constArrayRefFromArray(md->chargeA, md->nr)
                                             : gmx::ArrayRef<const real>{},
                                 md->chargeB ? gmx::constArrayRefFromArray(md->chargeB, md->nr)