Refactor scaled box handling for PME/Ewald

The box scaling code for Ewald and related logic is now encapsulated in
new box scaler class. This also allows box scaling moving box scaling
it further down the call chain so that logic is reduced in do_force.

Change-Id: Ic5071b825b9d36daca5f49d6f7c6c50261af1e1b

diff --git a/src/gromacs/ewald/ewald-utils.h b/src/gromacs/ewald/ewald-utils.h
index 078e1ef2d5fd0b18c54b4bdc5ad836e17e19d9c9..26b6c6b7a393a83da3cf7a06f787cf5e76ed28d9 100644 (file)
--- a/src/gromacs/ewald/ewald-utils.h
+++ b/src/gromacs/ewald/ewald-utils.h
@@ -36,18 +36,22 @@
  */
 /*! \libinternal \file
  *
- * \brief Declares functions for computing Ewald splitting coefficients
- *
- * These belong in the maths module because they do simple maths and
- * are used many parts of Gromacs.
+ * \brief Declares utility functions related to Ewald.
  *
  * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \author Szilárd Páll <pall.szilard@gmail.com>
+ *
  * \inlibraryapi
  * \ingroup module_ewald
  */
 #ifndef GMX_EWALD_UTILS_H
 #define GMX_EWALD_UTILS_H
 
+#include <assert.h>
+
+#include "gromacs/math/vec.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/pbcutil/pbc.h"
 #include "gromacs/utility/real.h"
 
 /*! \brief Computes the Ewald splitting coefficient for Coulomb
@@ -80,4 +84,59 @@ calc_ewaldcoeff_q(real rc, real rtol);
 real
 calc_ewaldcoeff_lj(real rc, real rtol);
 
+
+/*! \libinternal \brief Class to handle box scaling for Ewald and PME.
+ *
+ * At construction contents of inputrec determine whether scaling is necessary
+ * as well as the scaling factor used. Later, the scaleBox method can be used
+ * to apply the appropriate scaling (if needed) for Ewald-based methods.
+ *
+ */
+class EwaldBoxZScaler
+{
+
+    private:
+        bool scaleWithWalls_; /**< True if the simulation uses two walls and the box needs to be scaled in PME */
+        real scalingFactor_;  /**< Box The scaling factor PME uses with walls */
+
+    public:
+        EwaldBoxZScaler() = delete;
+
+        /*! \brief Constructor that takes the input record to initialize Ewald box scaling appropriately. */
+        EwaldBoxZScaler(const t_inputrec &ir)
+        {
+            if (inputrecPbcXY2Walls(&ir))
+            {
+                scaleWithWalls_ = true;
+                scalingFactor_  = ir.wall_ewald_zfac;
+            }
+            else
+            {
+                scaleWithWalls_ = false;
+                scalingFactor_  = 1;
+            }
+        }
+
+        /*! \brief Copy and scale the box for PME.
+         *
+         * When PME is used with 2D periodicity and two walls, the
+         * copy of the \p box passed is scaled with the Z scaling factor.
+         *
+         * \param[in] box        The current box matrix
+         * \param[out] scaledBox Scaled copy of the box matrix.
+         */
+        void scaleBox(const matrix box,
+                      matrix       scaledBox)
+        {
+            assert(box);
+            assert(scaledBox);
+
+            copy_mat(box, scaledBox);
+            if (scaleWithWalls_)
+            {
+                svmul(scalingFactor_, scaledBox[ZZ], scaledBox[ZZ]);
+            }
+        }
+};
+
 #endif

diff --git a/src/gromacs/ewald/ewald.cpp b/src/gromacs/ewald/ewald.cpp
index a628811aa906f0af385e465a4f08ca50b28d8b08..7471e5dd2524423c31d251abe8c24a9b9dca729d 100644 (file)
--- a/src/gromacs/ewald/ewald.cpp
+++ b/src/gromacs/ewald/ewald.cpp
@@ -56,6 +56,7 @@
 
 #include <algorithm>
 
+#include "gromacs/ewald/ewald-utils.h"
 #include "gromacs/math/functions.h"
 #include "gromacs/math/gmxcomplex.h"
 #include "gromacs/math/units.h"
@@ -137,14 +138,13 @@ static void tabulateStructureFactors(int natom, rvec x[], int kmax, cvec **eir,
 real do_ewald(t_inputrec *ir,
               rvec x[],        rvec f[],
               real chargeA[],  real chargeB[],
-              rvec box,
+              matrix box,
               t_commrec *cr,   int natoms,
               matrix lrvir,    real ewaldcoeff,
               real lambda,     real *dvdlambda,
               struct gmx_ewald_tab_t *et)
 {
     real     factor     = -1.0/(4*ewaldcoeff*ewaldcoeff);
-    real     scaleRecip = 4.0*M_PI/(box[XX]*box[YY]*box[ZZ])*ONE_4PI_EPS0/ir->epsilon_r; /* 1/(Vol*e0) */
     real    *charge, energy_AB[2], energy;
     rvec     lll;
     int      lowiy, lowiz, ix, iy, iz, n, q;
@@ -159,6 +159,19 @@ real do_ewald(t_inputrec *ir,
         }
     }
 
+    /* Scale box with Ewald wall factor */
+    matrix          scaledBox;
+    EwaldBoxZScaler boxScaler(*ir);
+    boxScaler.scaleBox(box, scaledBox);
+
+    rvec boxDiag;
+    for (int i = 0; (i < DIM); i++)
+    {
+        boxDiag[i] = scaledBox[i][i];
+    }
+
+    /* 1/(Vol*e0) */
+    real scaleRecip = 4.0*M_PI/(boxDiag[XX]*boxDiag[YY]*boxDiag[ZZ])*ONE_4PI_EPS0/ir->epsilon_r;
 
     if (!et->eir) /* allocate if we need to */
     {
@@ -175,7 +188,7 @@ real do_ewald(t_inputrec *ir,
 
     clear_mat(lrvir);
 
-    calc_lll(box, lll);
+    calc_lll(boxDiag, lll);
     tabulateStructureFactors(natoms, x, et->kmax, et->eir, lll);
 
     for (q = 0; q < (bFreeEnergy ? 2 : 1); q++)

diff --git a/src/gromacs/ewald/ewald.h b/src/gromacs/ewald/ewald.h
index fac97d511fe767b525fb36f6c80385b73e8870ff..6d31ab8380cd8abfc27398dae617c7773c90db39 100644 (file)
--- a/src/gromacs/ewald/ewald.h
+++ b/src/gromacs/ewald/ewald.h
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2017, 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.
@@ -87,7 +87,7 @@ real
 do_ewald(t_inputrec *ir,
          rvec x[],        rvec f[],
          real chargeA[],  real chargeB[],
-         rvec box,
+         matrix box,
          t_commrec *cr,  int natoms,
          matrix lrvir,   real ewaldcoeff,
          real lambda,    real *dvdlambda,

diff --git a/src/gromacs/ewald/long-range-correction.cpp b/src/gromacs/ewald/long-range-correction.cpp
index 09438a88c061eb5cf49dd25ae45a0794791fb188..05afae2c056ae65e78ce9e948b771b5617249f49 100644 (file)
--- a/src/gromacs/ewald/long-range-correction.cpp
+++ b/src/gromacs/ewald/long-range-correction.cpp
@@ -40,16 +40,20 @@
 
 #include <cmath>
 
+#include "gromacs/ewald/ewald-utils.h"
 #include "gromacs/math/functions.h"
 #include "gromacs/math/units.h"
 #include "gromacs/math/utilities.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdtypes/commrec.h"
 #include "gromacs/mdtypes/forcerec.h"
+#include "gromacs/mdtypes/inputrec.h"
 #include "gromacs/mdtypes/md_enums.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/gmxassert.h"
 
+#include "pme-internal.h"
+
 /* There's nothing special to do here if just masses are perturbed,
  * but if either charge or type is perturbed then the implementation
  * requires that B states are defined for both charge and type, and
@@ -64,6 +68,7 @@ void ewald_LRcorrection(int numAtomsLocal,
                         t_commrec *cr,
                         int numThreads, int thread,
                         t_forcerec *fr,
+                        const t_inputrec *ir,
                         real *chargeA, real *chargeB,
                         real *C6A, real *C6B,
                         real *sigmaA, real *sigmaB,
@@ -105,13 +110,17 @@ void ewald_LRcorrection(int numAtomsLocal,
     real        c6Ai   = 0, c6Bi = 0, c6A = 0, c6B = 0, ewcdr2, ewcdr4, c6L = 0, rinv6;
     rvec        df, dx, mutot[2], dipcorrA, dipcorrB;
     tensor      dxdf_q = {{0}}, dxdf_lj = {{0}};
-    real        vol    = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
     real        L1_q, L1_lj, dipole_coeff, qqA, qqB, qqL, vr0_q, vr0_lj = 0;
     real        chargecorr[2] = { 0, 0 };
     gmx_bool    bMolPBC       = fr->bMolPBC;
     gmx_bool    bDoingLBRule  = (fr->ljpme_combination_rule == eljpmeLB);
     gmx_bool    bNeedLongRangeCorrection;
 
+    assert(ir);
+    matrix          scaledBox;
+    EwaldBoxZScaler boxScaler(*ir);
+    boxScaler.scaleBox(box, scaledBox);
+
     /* This routine can be made faster by using tables instead of analytical interactions
      * However, that requires a thorough verification that they are correct in all cases.
      */
@@ -145,13 +154,15 @@ void ewald_LRcorrection(int numAtomsLocal,
         dipcorrB[i] = 0;
     }
     dipole_coeff = 0;
+
+    real boxVolume = scaledBox[XX][XX]*scaledBox[YY][YY]*scaledBox[ZZ][ZZ];
     switch (ewald_geometry)
     {
         case eewg3D:
             if (epsilon_surface != 0)
             {
                 dipole_coeff =
-                    2*M_PI*ONE_4PI_EPS0/((2*epsilon_surface + fr->ic->epsilon_r)*vol);
+                    2*M_PI*ONE_4PI_EPS0/((2*epsilon_surface + fr->ic->epsilon_r)*boxVolume);
                 for (i = 0; (i < DIM); i++)
                 {
                     dipcorrA[i] = 2*dipole_coeff*mutot[0][i];
@@ -160,7 +171,7 @@ void ewald_LRcorrection(int numAtomsLocal,
             }
             break;
         case eewg3DC:
-            dipole_coeff  = 2*M_PI*one_4pi_eps/vol;
+            dipole_coeff  = 2*M_PI*one_4pi_eps/boxVolume;
             dipcorrA[ZZ]  = 2*dipole_coeff*mutot[0][ZZ];
             dipcorrB[ZZ]  = 2*dipole_coeff*mutot[1][ZZ];
             for (int q = 0; q < (bHaveChargeOrTypePerturbed ? 2 : 1); q++)

diff --git a/src/gromacs/ewald/long-range-correction.h b/src/gromacs/ewald/long-range-correction.h
index 37ca4fb4d0422eda75f69b14f09d6eb74a5fc4f9..026189262a4a71490ab6eedb87a4b6acfca2abf7 100644 (file)
--- a/src/gromacs/ewald/long-range-correction.h
+++ b/src/gromacs/ewald/long-range-correction.h
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017, 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.
@@ -57,6 +57,8 @@
 #include "gromacs/utility/basedefinitions.h"
 #include "gromacs/utility/real.h"
 
+struct t_inputrec;
+
 /*! \brief Calculate long-range Ewald correction terms.
  *
  * For the group cutoff scheme (only), calculates the correction to
@@ -70,6 +72,7 @@ ewald_LRcorrection(int numAtomsLocal,
                    t_commrec *cr,
                    int numThreads, int thread,
                    t_forcerec *fr,
+                   const t_inputrec *ir,
                    real *chargeA, real *chargeB,
                    real *C6A, real *C6B,
                    real *sigmaA, real *sigmaB,

diff --git a/src/gromacs/ewald/pme-internal.h b/src/gromacs/ewald/pme-internal.h
index 6aeb4b72de5ba5d998b584bca85b16d8909d794f..fbb0dfe37a828b175863244aded1553483a4b522 100644 (file)
--- a/src/gromacs/ewald/pme-internal.h
+++ b/src/gromacs/ewald/pme-internal.h
@@ -55,9 +55,11 @@
 
 #include "config.h"
 
+#include <assert.h>
 #include <stdio.h>
 
 #include "gromacs/math/gmxcomplex.h"
+#include "gromacs/math/vec.h"
 #include "gromacs/timing/wallcycle.h"
 #include "gromacs/timing/walltime_accounting.h"
 #include "gromacs/utility/gmxmpi.h"
@@ -262,6 +264,8 @@ typedef struct gmx_pme_t {
     real       ewaldcoeff_lj; /* Ewald splitting coefficient for r^-6 */
     real       epsilon_r;
 
+    class EwaldBoxZScaler *boxScaler;   /*! The scaling data Ewald uses with walls (set at pme_init constant for the entire run) */
+
     int        ljpme_combination_rule;  /* Type of combination rule in LJ-PME */
 
     int        ngrids;                  /* number of grids we maintain for pmegrid, (c)fftgrid and pfft_setups*/

diff --git a/src/gromacs/ewald/pme-load-balancing.cpp b/src/gromacs/ewald/pme-load-balancing.cpp
index 58caf85b458f45eca50a6d04476c42dab9fd644b..d42337c1533f38c234e7cabdccc307d09b7e53e6 100644 (file)
--- a/src/gromacs/ewald/pme-load-balancing.cpp
+++ b/src/gromacs/ewald/pme-load-balancing.cpp
@@ -56,6 +56,7 @@
 #include "gromacs/domdec/domdec.h"
 #include "gromacs/domdec/domdec_network.h"
 #include "gromacs/domdec/domdec_struct.h"
+#include "gromacs/ewald/ewald-utils.h"
 #include "gromacs/ewald/pme.h"
 #include "gromacs/fft/calcgrid.h"
 #include "gromacs/gmxlib/network.h"
@@ -200,11 +201,11 @@ void pme_loadbal_init(pme_load_balancing_t     **pme_lb_p,
     pme_lb->rbufInner_coulomb = listParams->rlistInner - ic->rcoulomb;
     pme_lb->rbufInner_vdw     = listParams->rlistInner - ic->rvdw;
 
-    copy_mat(box, pme_lb->box_start);
-    if (ir->ePBC == epbcXY && ir->nwall == 2)
-    {
-        svmul(ir->wall_ewald_zfac, pme_lb->box_start[ZZ], pme_lb->box_start[ZZ]);
-    }
+    /* 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);
+    boxScaler.scaleBox(box, pme_lb->box_start);
 
     pme_lb->n = 1;
     snew(pme_lb->setup, pme_lb->n);
@@ -222,7 +223,11 @@ void pme_loadbal_init(pme_load_balancing_t     **pme_lb_p,
     pme_lb->setup[0].ewaldcoeff_q    = ic->ewaldcoeff_q;
     pme_lb->setup[0].ewaldcoeff_lj   = ic->ewaldcoeff_lj;
 
-    pme_lb->setup[0].pmedata         = pmedata;
+    if (!pme_lb->bSepPMERanks)
+    {
+        GMX_RELEASE_ASSERT(pmedata, "On ranks doing both PP and PME we need a valid pmedata object");
+        pme_lb->setup[0].pmedata     = pmedata;
+    }
 
     spm = 0;
     for (d = 0; d < DIM; d++)

diff --git a/src/gromacs/ewald/pme.cpp b/src/gromacs/ewald/pme.cpp
index 0ea7ae8feccbfef9f006bf16f0aaedeacb0e43e0..fdd9b4a69b26e26bf90e153cd991103a7d109c14 100644 (file)
--- a/src/gromacs/ewald/pme.cpp
+++ b/src/gromacs/ewald/pme.cpp
@@ -81,6 +81,7 @@
 
 #include <algorithm>
 
+#include "gromacs/ewald/ewald-utils.h"
 #include "gromacs/fft/parallel_3dfft.h"
 #include "gromacs/fileio/pdbio.h"
 #include "gromacs/gmxlib/network.h"
@@ -664,6 +665,11 @@ int gmx_pme_init(struct gmx_pme_t **pmedata,
     /* Always constant LJ coefficients */
     pme->ljpme_combination_rule = ir->ljpme_combination_rule;
 
+    // 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);
+
     /* If we violate restrictions, generate a fatal error here */
     gmx_pme_check_restrictions(pme->pme_order,
                                pme->nkx, pme->nky, pme->nkz,
@@ -1032,7 +1038,10 @@ int gmx_pme_do(struct gmx_pme_t *pme,
         atc->f = f;
     }
 
-    gmx::invertBoxMatrix(box, pme->recipbox);
+    matrix scaledBox;
+    pme->boxScaler->scaleBox(box, scaledBox);
+
+    gmx::invertBoxMatrix(scaledBox, pme->recipbox);
     bFirst = TRUE;
 
     /* For simplicity, we construct the splines for all particles if
@@ -1191,7 +1200,7 @@ int gmx_pme_do(struct gmx_pme_t *pme,
                     {
                         loop_count =
                             solve_pme_yzx(pme, cfftgrid,
-                                          box[XX][XX]*box[YY][YY]*box[ZZ][ZZ],
+                                          scaledBox[XX][XX]*scaledBox[YY][YY]*scaledBox[ZZ][ZZ],
                                           bCalcEnerVir,
                                           pme->nthread, thread);
                     }
@@ -1199,7 +1208,7 @@ int gmx_pme_do(struct gmx_pme_t *pme,
                     {
                         loop_count =
                             solve_pme_lj_yzx(pme, &cfftgrid, FALSE,
-                                             box[XX][XX]*box[YY][YY]*box[ZZ][ZZ],
+                                             scaledBox[XX][XX]*scaledBox[YY][YY]*scaledBox[ZZ][ZZ],
                                              bCalcEnerVir,
                                              pme->nthread, thread);
                     }
@@ -1470,7 +1479,7 @@ int gmx_pme_do(struct gmx_pme_t *pme,
 
                         loop_count =
                             solve_pme_lj_yzx(pme, &pme->cfftgrid[2], TRUE,
-                                             box[XX][XX]*box[YY][YY]*box[ZZ][ZZ],
+                                             scaledBox[XX][XX]*scaledBox[YY][YY]*scaledBox[ZZ][ZZ],
                                              bCalcEnerVir,
                                              pme->nthread, thread);
                         if (thread == 0)
@@ -1681,6 +1690,8 @@ void gmx_pme_destroy(gmx_pme_t *pme)
         return;
     }
 
+    delete pme->boxScaler;
+
     sfree(pme->nnx);
     sfree(pme->nny);
     sfree(pme->nnz);

diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp
index bc47e738fda926156d870671b6cd7babd0f59464..2db01e8e5a355a5790cc9cf59cfdaddb2aaf5378 100644 (file)
--- a/src/gromacs/gmxpreprocess/grompp.cpp
+++ b/src/gromacs/gmxpreprocess/grompp.cpp
@@ -49,6 +49,7 @@
 #include <sys/types.h>
 
 #include "gromacs/commandline/pargs.h"
+#include "gromacs/ewald/ewald-utils.h"
 #include "gromacs/ewald/pme.h"
 #include "gromacs/fft/calcgrid.h"
 #include "gromacs/fileio/confio.h"
@@ -1773,7 +1774,6 @@ int gmx_grompp(int argc, char *argv[])
     gpp_atomtype_t     atype;
     int                nvsite, comb, mt;
     t_params          *plist;
-    matrix             box;
     real               fudgeQQ;
     double             reppow;
     const char        *mdparin;
@@ -2238,11 +2238,10 @@ int gmx_grompp(int argc, char *argv[])
     if (EEL_FULL(ir->coulombtype) || EVDW_PME(ir->vdwtype))
     {
         /* Calculate the optimal grid dimensions */
-        copy_mat(state.box, box);
-        if (ir->ePBC == epbcXY && ir->nwall == 2)
-        {
-            svmul(ir->wall_ewald_zfac, box[ZZ], box[ZZ]);
-        }
+        matrix          scaledBox;
+        EwaldBoxZScaler boxScaler(*ir);
+        boxScaler.scaleBox(state.box, scaledBox);
+
         if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0)
         {
             /* Mark fourier_spacing as not used */
@@ -2254,7 +2253,7 @@ int gmx_grompp(int argc, char *argv[])
             warning_error(wi, "Some of the Fourier grid sizes are set, but all of them need to be set.");
         }
         const int minGridSize = minimalPmeGridSize(ir->pme_order);
-        calcFftGrid(stdout, box, ir->fourier_spacing, minGridSize,
+        calcFftGrid(stdout, scaledBox, ir->fourier_spacing, minGridSize,
                     &(ir->nkx), &(ir->nky), &(ir->nkz));
         if (ir->nkx < minGridSize ||
             ir->nky < minGridSize ||

diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp
index 5e8c9257ea28bb0be47f410c2d177b1be7f48985..5f1224723a7522cc45a57b312045252e8c44e02f 100644 (file)
--- a/src/gromacs/mdlib/force.cpp
+++ b/src/gromacs/mdlib/force.cpp
@@ -149,10 +149,7 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
 {
     int         i, j;
     int         donb_flags;
-    gmx_bool    bSB;
     int         pme_flags;
-    matrix      boxs;
-    rvec        box_size;
     t_pbc       pbc;
     real        dvdl_dum[efptNR], dvdl_nb[efptNR];
 
@@ -169,12 +166,6 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
         dvdl_dum[i] = 0;
     }
 
-    /* Reset box */
-    for (i = 0; (i < DIM); i++)
-    {
-        box_size[i] = box[i][i];
-    }
-
     /* do QMMM first if requested */
     if (fr->bQMMM)
     {
@@ -380,14 +371,6 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
         real Vlr_q             = 0, Vlr_lj = 0, Vcorr_q = 0, Vcorr_lj = 0;
         real dvdl_long_range_q = 0, dvdl_long_range_lj = 0;
 
-        bSB = (ir->nwall == 2);
-        if (bSB)
-        {
-            copy_mat(box, boxs);
-            svmul(ir->wall_ewald_zfac, boxs[ZZ], boxs[ZZ]);
-            box_size[ZZ] *= ir->wall_ewald_zfac;
-        }
-
         if (EEL_PME_EWALD(fr->ic->eeltype) || EVDW_PME(fr->ic->vdwtype))
         {
             real dvdl_long_range_correction_q   = 0;
@@ -449,14 +432,14 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
                          * exclusion forces) are calculated, so we can store
                          * the forces in the normal, single fr->f_novirsum array.
                          */
-                        ewald_LRcorrection(md->homenr, cr, nthreads, t, fr,
+                        ewald_LRcorrection(md->homenr, cr, nthreads, t, fr, ir,
                                            md->chargeA, md->chargeB,
                                            md->sqrt_c6A, md->sqrt_c6B,
                                            md->sigmaA, md->sigmaB,
                                            md->sigma3A, md->sigma3B,
                                            md->nChargePerturbed || md->nTypePerturbed,
                                            ir->cutoff_scheme != ecutsVERLET,
-                                           excl, x, bSB ? boxs : box, mu_tot,
+                                           excl, x, box, mu_tot,
                                            ir->ewald_geometry,
                                            ir->epsilon_surface,
                                            as_rvec_array(fr->f_novirsum->data()),
@@ -529,7 +512,7 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
                                         md->chargeA, md->chargeB,
                                         md->sqrt_c6A, md->sqrt_c6B,
                                         md->sigmaA, md->sigmaB,
-                                        bSB ? boxs : box, cr,
+                                        box, cr,
                                         DOMAINDECOMP(cr) ? dd_pme_maxshift_x(cr->dd) : 0,
                                         DOMAINDECOMP(cr) ? dd_pme_maxshift_y(cr->dd) : 0,
                                         nrnb, wcycle,
@@ -572,7 +555,7 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
         {
             Vlr_q = do_ewald(ir, x, as_rvec_array(fr->f_novirsum->data()),
                              md->chargeA, md->chargeB,
-                             box_size, cr, md->homenr,
+                             box, cr, md->homenr,
                              fr->vir_el_recip, fr->ic->ewaldcoeff_q,
                              lambda[efptCOUL], &dvdl_long_range_q, fr->ewald_table);
         }

diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp
index 815b1d964d81dcfb0b8ccae37ee46ac068f5cb05..a29153b2792fa019c94fcc38642b2c92333db703 100644 (file)
--- a/src/gromacs/mdlib/sim_util.cpp
+++ b/src/gromacs/mdlib/sim_util.cpp
@@ -909,29 +909,16 @@ static void do_force_cutsVERLET(FILE *fplog, t_commrec *cr,
 #if GMX_MPI
     if (!(cr->duty & DUTY_PME))
     {
-        gmx_bool bBS;
-        matrix   boxs;
-
         /* Send particle coordinates to the pme nodes.
          * Since this is only implemented for domain decomposition
          * and domain decomposition does not use the graph,
          * we do not need to worry about shifting.
          */
-
         wallcycle_start(wcycle, ewcPP_PMESENDX);
-
-        bBS = (inputrec->nwall == 2);
-        if (bBS)
-        {
-            copy_mat(box, boxs);
-            svmul(inputrec->wall_ewald_zfac, boxs[ZZ], boxs[ZZ]);
-        }
-
-        gmx_pme_send_coordinates(cr, bBS ? boxs : box, x,
+        gmx_pme_send_coordinates(cr, box, x,
                                  lambda[efptCOUL], lambda[efptVDW],
                                  (flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY)),
                                  step);
-
         wallcycle_stop(wcycle, ewcPP_PMESENDX);
     }
 #endif /* GMX_MPI */
@@ -1591,29 +1578,16 @@ static void do_force_cutsGROUP(FILE *fplog, t_commrec *cr,
 #if GMX_MPI
     if (!(cr->duty & DUTY_PME))
     {
-        gmx_bool bBS;
-        matrix   boxs;
-
         /* Send particle coordinates to the pme nodes.
          * Since this is only implemented for domain decomposition
          * and domain decomposition does not use the graph,
          * we do not need to worry about shifting.
          */
-
         wallcycle_start(wcycle, ewcPP_PMESENDX);
-
-        bBS = (inputrec->nwall == 2);
-        if (bBS)
-        {
-            copy_mat(box, boxs);
-            svmul(inputrec->wall_ewald_zfac, boxs[ZZ], boxs[ZZ]);
-        }
-
-        gmx_pme_send_coordinates(cr, bBS ? boxs : box, x,
+        gmx_pme_send_coordinates(cr, box, x,
                                  lambda[efptCOUL], lambda[efptVDW],
                                  (flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY)),
                                  step);
-
         wallcycle_stop(wcycle, ewcPP_PMESENDX);
     }
 #endif /* GMX_MPI */

diff --git a/src/gromacs/mdtypes/inputrec.cpp b/src/gromacs/mdtypes/inputrec.cpp
index 0723baeb4eeedf17f0210dece3374a84fa6664b4..753165ab6d470da27c377f7cd2542cd9fa596606 100644 (file)
--- a/src/gromacs/mdtypes/inputrec.cpp
+++ b/src/gromacs/mdtypes/inputrec.cpp
@@ -1347,6 +1347,11 @@ gmx_bool inputrecNphTrotter(const t_inputrec *ir)
              (ir->epc == epcMTTK) && (ir->etc != etcNOSEHOOVER) );
 }
 
+bool inputrecPbcXY2Walls(const t_inputrec *ir)
+{
+    return (ir->ePBC == epbcXY && ir->nwall == 2);
+}
+
 bool integratorHasConservedEnergyQuantity(const t_inputrec *ir)
 {
     if (!EI_MD(ir->eI))

diff --git a/src/gromacs/mdtypes/inputrec.h b/src/gromacs/mdtypes/inputrec.h
index ec00b44d952e22d93abd0d00089fdb4381282589..f717846a610ac15ec815d8f392d6d90b3a00b389 100644 (file)
--- a/src/gromacs/mdtypes/inputrec.h
+++ b/src/gromacs/mdtypes/inputrec.h
@@ -459,6 +459,9 @@ gmx_bool inputrecNvtTrotter(const t_inputrec *ir);
 
 gmx_bool inputrecNphTrotter(const t_inputrec *ir);
 
+/*! \brief Return true if the simulation is 2D periodic with two walls. */
+bool     inputrecPbcXY2Walls(const t_inputrec *ir);
+
 /* Returns true for MD integator with T and/or P-coupling that supports
  * calculating the conserved energy quantity.
  */