Remove forcerec from ewald_charge_correction
authorejjordan <ejjordan@kth.se>
Tue, 18 May 2021 18:35:53 +0000 (20:35 +0200)
committerPaul Bauer <paul.bauer.q@gmail.com>
Wed, 19 May 2021 08:42:46 +0000 (08:42 +0000)
Also renamed commrec to be more verbose.

src/gromacs/ewald/ewald.cpp
src/gromacs/ewald/ewald.h
src/gromacs/mdlib/force.cpp

index 3a62f306d1f906baf7ea398209fc52d02a0357e4..b13e3e7c365e51a1673600236c774400464fb15f 100644 (file)
@@ -314,47 +314,38 @@ real do_ewald(const t_inputrec&              ir,
     return energy;
 }
 
-real ewald_charge_correction(const t_commrec*  cr,
-                             const t_forcerec* fr,
-                             const real        lambda,
-                             const matrix      box,
-                             real*             dvdlambda,
-                             tensor            vir)
+real ewald_charge_correction(const t_commrec*            commrec,
+                             const real                  epsilonR,
+                             const real                  ewaldcoeffQ,
+                             gmx::ArrayRef<const double> qsum,
+                             const real                  lambda,
+                             const matrix                box,
+                             real*                       dvdlambda,
+                             tensor                      vir)
 
 {
-    real vol, fac, qs2A, qs2B, vc, enercorr;
-    int  d;
+    real enercorr = 0;
 
-    if (MASTER(cr))
+    if (MASTER(commrec))
     {
         /* Apply charge correction */
-        vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
+        real vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
 
-        fac = M_PI * gmx::c_one4PiEps0
-              / (fr->ic->epsilon_r * 2.0 * vol * vol * gmx::square(fr->ic->ewaldcoeff_q));
+        real fac = M_PI * gmx::c_one4PiEps0 / (epsilonR * 2.0 * vol * vol * gmx::square(ewaldcoeffQ));
 
-        qs2A = fr->qsum[0] * fr->qsum[0];
-        qs2B = fr->qsum[1] * fr->qsum[1];
+        real qs2A = qsum[0] * qsum[0];
+        real qs2B = qsum[1] * qsum[1];
 
-        vc = (qs2A * (1 - lambda) + qs2B * lambda) * fac;
+        real vc = (qs2A * (1 - lambda) + qs2B * lambda) * fac;
 
         enercorr = -vol * vc;
 
         *dvdlambda += -vol * (qs2B - qs2A) * fac;
 
-        for (d = 0; d < DIM; d++)
+        for (int d = 0; d < DIM; d++)
         {
             vir[d][d] += vc;
         }
-
-        if (debug)
-        {
-            fprintf(debug, "Total charge correction: Vcharge=%g\n", enercorr);
-        }
-    }
-    else
-    {
-        enercorr = 0;
     }
 
     return enercorr;
index 2d4d20bc98c8ded83c0253c6693eb2963d27c4b0..cfb6cd2a866e1f2466faabc143ab9f7625f17db5 100644 (file)
@@ -118,11 +118,13 @@ real do_ewald(const t_inputrec&              ir,
  * charge.
  *
  * Should only be called on one thread. */
-real ewald_charge_correction(const t_commrec*  cr,
-                             const t_forcerec* fr,
-                             real              lambda,
-                             const matrix      box,
-                             real*             dvdlambda,
-                             tensor            vir);
+real ewald_charge_correction(const t_commrec*            commrec,
+                             real                        epsilonR,
+                             real                        ewaldcoeffQ,
+                             gmx::ArrayRef<const double> qsum,
+                             real                        lambda,
+                             const matrix                box,
+                             real*                       dvdlambda,
+                             tensor                      vir);
 
 #endif
index a939841283a504e4231d7bf914b51128adb81464..8f07ca623113671f55b66c0b7107b417e1bc00e2 100644 (file)
@@ -195,7 +195,9 @@ void calculateLongRangeNonbondeds(t_forcerec*                    fr,
                    negligible and constant-sized amount of time */
                 ewaldOutput.Vcorr_q += ewald_charge_correction(
                         cr,
-                        fr,
+                        fr->ic->epsilon_r,
+                        fr->ic->ewaldcoeff_q,
+                        fr->qsum,
                         lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)],
                         box,
                         &ewaldOutput.dvdl[FreeEnergyPerturbationCouplingType::Coul],