Remove inputrec from do_ewald
authorejjordan <ejjordan@kth.se>
Wed, 19 May 2021 16:28:07 +0000 (18:28 +0200)
committerPaul Bauer <paul.bauer.q@gmail.com>
Thu, 20 May 2021 09:15:25 +0000 (09:15 +0000)
Also rename some parameters more verbosely and remove a dead code
block.

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

index 4d04a4f512e059beaefc8dcd1fa39a5240b94c67..65e31644c3c71587179efd8b7c30bf1a55e2b5e0 100644 (file)
@@ -131,13 +131,16 @@ static void tabulateStructureFactors(int natom, gmx::ArrayRef<const gmx::RVec> x
     }
 }
 
-real do_ewald(const t_inputrec&              ir,
-              gmx::ArrayRef<const gmx::RVec> x,
-              gmx::ArrayRef<gmx::RVec>       f,
+real do_ewald(bool                           havePbcXY2Walls,
+              real                           wallEwaldZfac,
+              real                           epsilonR,
+              FreeEnergyPerturbationType     freeEnergyPerturbationType,
+              gmx::ArrayRef<const gmx::RVec> coords,
+              gmx::ArrayRef<gmx::RVec>       forces,
               gmx::ArrayRef<const real>      chargeA,
               gmx::ArrayRef<const real>      chargeB,
               const matrix                   box,
-              const t_commrec*               cr,
+              const t_commrec*               commrec,
               int                            natoms,
               matrix                         lrvir,
               real                           ewaldcoeff,
@@ -153,9 +156,9 @@ real do_ewald(const t_inputrec&              ir,
     cvec** eir;
     bool   bFreeEnergy;
 
-    if (cr != nullptr)
+    if (commrec != nullptr)
     {
-        if (PAR(cr))
+        if (PAR(commrec))
         {
             gmx_fatal(FARGS, "No parallel Ewald. Use PME instead.\n");
         }
@@ -163,7 +166,7 @@ real do_ewald(const t_inputrec&              ir,
 
     /* Scale box with Ewald wall factor */
     matrix          scaledBox;
-    EwaldBoxZScaler boxScaler(inputrecPbcXY2Walls(&ir), ir.wall_ewald_zfac);
+    EwaldBoxZScaler boxScaler(havePbcXY2Walls, wallEwaldZfac);
     boxScaler.scaleBox(box, scaledBox);
 
     rvec boxDiag;
@@ -173,8 +176,7 @@ real do_ewald(const t_inputrec&              ir,
     }
 
     /* 1/(Vol*e0) */
-    real scaleRecip =
-            4.0 * M_PI / (boxDiag[XX] * boxDiag[YY] * boxDiag[ZZ]) * gmx::c_one4PiEps0 / ir.epsilon_r;
+    real scaleRecip = 4.0 * M_PI / (boxDiag[XX] * boxDiag[YY] * boxDiag[ZZ]) * gmx::c_one4PiEps0 / epsilonR;
 
     snew(eir, et->kmax);
     for (n = 0; n < et->kmax; n++)
@@ -184,12 +186,12 @@ real do_ewald(const t_inputrec&              ir,
     et->tab_xy.resize(natoms);
     et->tab_qxyz.resize(natoms);
 
-    bFreeEnergy = (ir.efep != FreeEnergyPerturbationType::No);
+    bFreeEnergy = (freeEnergyPerturbationType != FreeEnergyPerturbationType::No);
 
     clear_mat(lrvir);
 
     calc_lll(boxDiag, lll);
-    tabulateStructureFactors(natoms, x, et->kmax, eir, lll);
+    tabulateStructureFactors(natoms, coords, et->kmax, eir, lll);
 
     gmx::ArrayRef<const real> charge;
     for (q = 0; q < (bFreeEnergy ? 2 : 1); q++)
@@ -272,14 +274,9 @@ real do_ewald(const t_inputrec&              ir,
                     {
                         /*tmp=scale*ak*(cs*tab_qxyz[n].im-ss*tab_qxyz[n].re);*/
                         tmp = scale * ak * (cs * et->tab_qxyz[n].im - ss * et->tab_qxyz[n].re);
-                        f[n][XX] += tmp * mx * 2 * scaleRecip;
-                        f[n][YY] += tmp * my * 2 * scaleRecip;
-                        f[n][ZZ] += tmp * mz * 2 * scaleRecip;
-#if 0
-                        f[n][XX] += tmp*mx;
-                        f[n][YY] += tmp*my;
-                        f[n][ZZ] += tmp*mz;
-#endif
+                        forces[n][XX] += tmp * mx * 2 * scaleRecip;
+                        forces[n][YY] += tmp * my * 2 * scaleRecip;
+                        forces[n][ZZ] += tmp * mz * 2 * scaleRecip;
                     }
                     lowiz = 1 - et->nz;
                 }
index cfb6cd2a866e1f2466faabc143ab9f7625f17db5..809bc8f4f78f981b079098d21cb6b55e0994ac9e 100644 (file)
@@ -77,6 +77,7 @@ struct t_commrec;
 struct t_forcerec;
 struct t_inputrec;
 struct t_complex;
+enum class FreeEnergyPerturbationType : int;
 
 namespace gmx
 {
@@ -100,13 +101,16 @@ struct gmx_ewald_tab_t
 };
 
 /*! \brief Do the long-ranged part of an Ewald calculation */
-real do_ewald(const t_inputrec&              ir,
-              gmx::ArrayRef<const gmx::RVec> x,
-              gmx::ArrayRef<gmx::RVec>       f,
+real do_ewald(bool                           havePbcXY2Walls,
+              real                           wallEwaldZfac,
+              real                           epsilonR,
+              FreeEnergyPerturbationType     freeEnergyPerturbationType,
+              gmx::ArrayRef<const gmx::RVec> coords,
+              gmx::ArrayRef<gmx::RVec>       forces,
               gmx::ArrayRef<const real>      chargeA,
               gmx::ArrayRef<const real>      chargeB,
               const matrix                   box,
-              const t_commrec*               cr,
+              const t_commrec*               commrec,
               int                            natoms,
               matrix                         lrvir,
               real                           ewaldcoeff,
index 11123448175767bc13e72cb3401d81aa8373e994..3ac2484a2b37e21ca42acbecbf1e78bfc179aae3 100644 (file)
@@ -281,7 +281,10 @@ void calculateLongRangeNonbondeds(t_forcerec*                    fr,
 
         if (fr->ic->eeltype == CoulombInteractionType::Ewald)
         {
-            Vlr_q = do_ewald(ir,
+            Vlr_q = do_ewald(inputrecPbcXY2Walls(&ir),
+                             ir.wall_ewald_zfac,
+                             ir.epsilon_r,
+                             ir.efep,
                              coordinates,
                              forceWithVirial->force_,
                              gmx::arrayRefFromArray(md->chargeA, md->nr),