Move the PME Ewald coefficients to the PME structure

Ewald coefficients no longer get passed every step to gmx_pme_do().
Instead, they are passed to gmx_pme_(re)init,
as they only change together with the grid size, not every MD step.

Change-Id: Iea75063cf060317e3f0fcfd147431520b6c83a90

diff --git a/src/gromacs/ewald/pme-internal.h b/src/gromacs/ewald/pme-internal.h
index 8f9832c11387f433be0abf49791eb8ac2bc9f8c5..dfa7240e947e0563f6d9a88ed8d3803f1f2e45a3 100644 (file)
--- a/src/gromacs/ewald/pme-internal.h
+++ b/src/gromacs/ewald/pme-internal.h
@@ -94,12 +94,17 @@ static const real lb_scale_factor_symm[] = { 2.0/64, 12.0/64, 30.0/64, 20.0/64 }
  */
 #define PME_ORDER_MAX 12
 
-/*! \brief As gmx_pme_init, but takes most settings, except the grid, from pme_src */
+/*! \brief As gmx_pme_init, but takes most settings, except the grid/Ewald coefficients, from pme_src.
+ * This is only called when the PME cut-off/grid size changes.
+ */
 int gmx_pme_reinit(struct gmx_pme_t **pmedata,
                    t_commrec *        cr,
                    struct gmx_pme_t * pme_src,
                    const t_inputrec * ir,
-                   ivec               grid_size);
+                   ivec               grid_size,
+                   real               ewaldcoeff_q,
+                   real               ewaldcoeff_lj);
+
 
 /* The following three routines are for PME/PP node splitting in pme_pp.c */
 
@@ -252,6 +257,8 @@ typedef struct gmx_pme_t {
     int        nkx, nky, nkz; /* Grid dimensions */
     gmx_bool   bP3M;          /* Do P3M: optimize the influence function */
     int        pme_order;
+    real       ewaldcoeff_q;  /* Ewald splitting coefficient for Coulomb */
+    real       ewaldcoeff_lj; /* Ewald splitting coefficient for r^-6 */
     real       epsilon_r;
 
     int        ljpme_combination_rule;  /* Type of combination rule in LJ-PME */

diff --git a/src/gromacs/ewald/pme-load-balancing.cpp b/src/gromacs/ewald/pme-load-balancing.cpp
index c7f284e5ac05a2e1f02b404e2fe6e227712dce23..ed4e6e038dc6430ce10f6d3a0fece02f219a8adf 100644 (file)
--- a/src/gromacs/ewald/pme-load-balancing.cpp
+++ b/src/gromacs/ewald/pme-load-balancing.cpp
@@ -822,7 +822,7 @@ pme_load_balance(pme_load_balancing_t      *pme_lb,
              */
             gmx_pme_reinit(&set->pmedata,
                            cr, pme_lb->setup[0].pmedata, ir,
-                           set->grid);
+                           set->grid, set->ewaldcoeff_q, set->ewaldcoeff_lj);
         }
         *pmedata = set->pmedata;
     }

diff --git a/src/gromacs/ewald/pme-only.cpp b/src/gromacs/ewald/pme-only.cpp
index acb78c6812d827e1cbf04330f478080b7867b53d..6d8c0984357626a346079437b562ad4a3adcfcf2 100644 (file)
--- a/src/gromacs/ewald/pme-only.cpp
+++ b/src/gromacs/ewald/pme-only.cpp
@@ -107,6 +107,7 @@ static void reset_pmeonly_counters(gmx_wallcycle_t wcycle,
 
 static void gmx_pmeonly_switch(int *npmedata, struct gmx_pme_t ***pmedata,
                                ivec grid_size,
+                               real ewaldcoeff_q, real ewaldcoeff_lj,
                                t_commrec *cr, t_inputrec *ir,
                                struct gmx_pme_t **pme_ret)
 {
@@ -133,7 +134,7 @@ static void gmx_pmeonly_switch(int *npmedata, struct gmx_pme_t ***pmedata,
     srenew(*pmedata, *npmedata);
 
     /* Generate a new PME data structure, copying part of the old pointers */
-    gmx_pme_reinit(&((*pmedata)[ind]), cr, pme, ir, grid_size);
+    gmx_pme_reinit(&((*pmedata)[ind]), cr, pme, ir, grid_size, ewaldcoeff_q, ewaldcoeff_lj);
 
     *pme_ret = (*pmedata)[ind];
 }
@@ -197,7 +198,7 @@ int gmx_pmeonly(struct gmx_pme_t *pme,
             if (ret == pmerecvqxSWITCHGRID)
             {
                 /* Switch the PME grid to grid_switch */
-                gmx_pmeonly_switch(&npmedata, &pmedata, grid_switch, cr, ir, &pme);
+                gmx_pmeonly_switch(&npmedata, &pmedata, grid_switch, ewaldcoeff_q, ewaldcoeff_lj, cr, ir, &pme);
             }
 
             if (ret == pmerecvqxRESETCOUNTERS)
@@ -230,7 +231,7 @@ int gmx_pmeonly(struct gmx_pme_t *pme,
         gmx_pme_do(pme, 0, natoms, x_pp, f_pp,
                    chargeA, chargeB, c6A, c6B, sigmaA, sigmaB, box,
                    cr, maxshift_x, maxshift_y, mynrnb, wcycle,
-                   vir_q, ewaldcoeff_q, vir_lj, ewaldcoeff_lj,
+                   vir_q, vir_lj,
                    &energy_q, &energy_lj, lambda_q, lambda_lj, &dvdlambda_q, &dvdlambda_lj,
                    GMX_PME_DO_ALL_F | (bEnerVir ? GMX_PME_CALC_ENER_VIR : 0));
 

diff --git a/src/gromacs/ewald/pme-solve.cpp b/src/gromacs/ewald/pme-solve.cpp
index b9bef2570f6003c4b92de590e3f70e048e97cd0a..71553ad7a2d51ee41a2b502a3d898649df4c0afb 100644 (file)
--- a/src/gromacs/ewald/pme-solve.cpp
+++ b/src/gromacs/ewald/pme-solve.cpp
@@ -287,8 +287,7 @@ gmx_inline static void calc_exponentials_lj(int start, int end, real *r, real *t
 }
 #endif
 
-int solve_pme_yzx(struct gmx_pme_t *pme, t_complex *grid,
-                  real ewaldcoeff, real vol,
+int solve_pme_yzx(struct gmx_pme_t *pme, t_complex *grid, real vol,
                   gmx_bool bEnerVir,
                   int nthread, int thread)
 {
@@ -298,7 +297,8 @@ int solve_pme_yzx(struct gmx_pme_t *pme, t_complex *grid,
     int                      kx, ky, kz, maxkx, maxky;
     int                      nx, ny, nz, iyz0, iyz1, iyz, iy, iz, kxstart, kxend;
     real                     mx, my, mz;
-    real                     factor = M_PI*M_PI/(ewaldcoeff*ewaldcoeff);
+    real                     ewaldcoeff = pme->ewaldcoeff_q;
+    real                     factor     = M_PI*M_PI/(ewaldcoeff*ewaldcoeff);
     real                     ets2, struct2, vfactor, ets2vf;
     real                     d1, d2, energy = 0;
     real                     by, bz;
@@ -538,8 +538,7 @@ int solve_pme_yzx(struct gmx_pme_t *pme, t_complex *grid,
     return local_ndata[YY]*local_ndata[XX];
 }
 
-int solve_pme_lj_yzx(struct gmx_pme_t *pme, t_complex **grid, gmx_bool bLB,
-                     real ewaldcoeff, real vol,
+int solve_pme_lj_yzx(struct gmx_pme_t *pme, t_complex **grid, gmx_bool bLB, real vol,
                      gmx_bool bEnerVir, int nthread, int thread)
 {
     /* do recip sum over local cells in grid */
@@ -548,7 +547,8 @@ int solve_pme_lj_yzx(struct gmx_pme_t *pme, t_complex **grid, gmx_bool bLB,
     int                      kx, ky, kz, maxkx, maxky;
     int                      nx, ny, nz, iy, iyz0, iyz1, iyz, iz, kxstart, kxend;
     real                     mx, my, mz;
-    real                     factor = M_PI*M_PI/(ewaldcoeff*ewaldcoeff);
+    real                     ewaldcoeff = pme->ewaldcoeff_lj;
+    real                     factor     = M_PI*M_PI/(ewaldcoeff*ewaldcoeff);
     real                     ets2, ets2vf;
     real                     eterm, vterm, d1, d2, energy = 0;
     real                     by, bz;

diff --git a/src/gromacs/ewald/pme-solve.h b/src/gromacs/ewald/pme-solve.h
index 6b26a58f1bc37df9b0ac46c5475d0b663b22ae5a..33c828a39d58db951039b6dca8a3945739d699ad 100644 (file)
--- a/src/gromacs/ewald/pme-solve.h
+++ b/src/gromacs/ewald/pme-solve.h
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016, 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.
@@ -70,12 +70,12 @@ void get_pme_ener_vir_lj(struct pme_solve_work_t *work, int nthread,
                          real *mesh_energy, matrix vir);
 
 int solve_pme_yzx(struct gmx_pme_t *pme, t_complex *grid,
-                  real ewaldcoeff, real vol,
+                  real vol,
                   gmx_bool bEnerVir,
                   int nthread, int thread);
 
 int solve_pme_lj_yzx(struct gmx_pme_t *pme, t_complex **grid, gmx_bool bLB,
-                     real ewaldcoeff, real vol,
+                     real vol,
                      gmx_bool bEnerVir, int nthread, int thread);
 
 #endif

diff --git a/src/gromacs/ewald/pme.cpp b/src/gromacs/ewald/pme.cpp
index 859d5256ab4bc03b4a0663b7cdbb33d270ccf712..fa65128965d588c9370a18250a75f934c27c5a41 100644 (file)
--- a/src/gromacs/ewald/pme.cpp
+++ b/src/gromacs/ewald/pme.cpp
@@ -479,6 +479,8 @@ int gmx_pme_init(struct gmx_pme_t **pmedata,
                  gmx_bool           bFreeEnergy_q,
                  gmx_bool           bFreeEnergy_lj,
                  gmx_bool           bReproducible,
+                 real               ewaldcoeff_q,
+                 real               ewaldcoeff_lj,
                  int                nthread)
 {
     struct gmx_pme_t *pme = NULL;
@@ -609,19 +611,21 @@ int gmx_pme_init(struct gmx_pme_t **pmedata,
      * not calculating free-energy for Coulomb and/or LJ while gmx_pme_init()
      * configures with free-energy, but that has never been tested.
      */
-    pme->doCoulomb   = EEL_PME(ir->coulombtype);
-    pme->doLJ        = EVDW_PME(ir->vdwtype);
-    pme->bFEP_q      = ((ir->efep != efepNO) && bFreeEnergy_q);
-    pme->bFEP_lj     = ((ir->efep != efepNO) && bFreeEnergy_lj);
-    pme->bFEP        = (pme->bFEP_q || pme->bFEP_lj);
-    pme->nkx         = ir->nkx;
-    pme->nky         = ir->nky;
-    pme->nkz         = ir->nkz;
-    pme->bP3M        = (ir->coulombtype == eelP3M_AD || getenv("GMX_PME_P3M") != NULL);
-    pme->pme_order   = ir->pme_order;
+    pme->doCoulomb     = EEL_PME(ir->coulombtype);
+    pme->doLJ          = EVDW_PME(ir->vdwtype);
+    pme->bFEP_q        = ((ir->efep != efepNO) && bFreeEnergy_q);
+    pme->bFEP_lj       = ((ir->efep != efepNO) && bFreeEnergy_lj);
+    pme->bFEP          = (pme->bFEP_q || pme->bFEP_lj);
+    pme->nkx           = ir->nkx;
+    pme->nky           = ir->nky;
+    pme->nkz           = ir->nkz;
+    pme->bP3M          = (ir->coulombtype == eelP3M_AD || getenv("GMX_PME_P3M") != NULL);
+    pme->pme_order     = ir->pme_order;
+    pme->ewaldcoeff_q  = ewaldcoeff_q;
+    pme->ewaldcoeff_lj = ewaldcoeff_lj;
 
     /* Always constant electrostatics coefficients */
-    pme->epsilon_r   = ir->epsilon_r;
+    pme->epsilon_r     = ir->epsilon_r;
 
     /* Always constant LJ coefficients */
     pme->ljpme_combination_rule = ir->ljpme_combination_rule;
@@ -815,7 +819,9 @@ int gmx_pme_reinit(struct gmx_pme_t **pmedata,
                    t_commrec *        cr,
                    struct gmx_pme_t * pme_src,
                    const t_inputrec * ir,
-                   ivec               grid_size)
+                   ivec               grid_size,
+                   real               ewaldcoeff_q,
+                   real               ewaldcoeff_lj)
 {
     t_inputrec irc;
     int        homenr;
@@ -836,7 +842,7 @@ int gmx_pme_reinit(struct gmx_pme_t **pmedata,
     }
 
     ret = gmx_pme_init(pmedata, cr, pme_src->nnodes_major, pme_src->nnodes_minor,
-                       &irc, homenr, pme_src->bFEP_q, pme_src->bFEP_lj, FALSE, pme_src->nthread);
+                       &irc, homenr, pme_src->bFEP_q, pme_src->bFEP_lj, FALSE, ewaldcoeff_q, ewaldcoeff_lj, pme_src->nthread);
 
     if (ret == 0)
     {
@@ -921,8 +927,7 @@ int gmx_pme_do(struct gmx_pme_t *pme,
                matrix box,      t_commrec *cr,
                int  maxshift_x, int maxshift_y,
                t_nrnb *nrnb,    gmx_wallcycle_t wcycle,
-               matrix vir_q,    real ewaldcoeff_q,
-               matrix vir_lj,   real ewaldcoeff_lj,
+               matrix vir_q,    matrix vir_lj,
                real *energy_q,  real *energy_lj,
                real lambda_q,   real lambda_lj,
                real *dvdlambda_q, real *dvdlambda_lj,
@@ -1139,7 +1144,7 @@ int gmx_pme_do(struct gmx_pme_t *pme,
                     if (grid_index < DO_Q)
                     {
                         loop_count =
-                            solve_pme_yzx(pme, cfftgrid, ewaldcoeff_q,
+                            solve_pme_yzx(pme, cfftgrid,
                                           box[XX][XX]*box[YY][YY]*box[ZZ][ZZ],
                                           bCalcEnerVir,
                                           pme->nthread, thread);
@@ -1147,7 +1152,7 @@ int gmx_pme_do(struct gmx_pme_t *pme,
                     else
                     {
                         loop_count =
-                            solve_pme_lj_yzx(pme, &cfftgrid, FALSE, ewaldcoeff_lj,
+                            solve_pme_lj_yzx(pme, &cfftgrid, FALSE,
                                              box[XX][XX]*box[YY][YY]*box[ZZ][ZZ],
                                              bCalcEnerVir,
                                              pme->nthread, thread);
@@ -1418,7 +1423,7 @@ int gmx_pme_do(struct gmx_pme_t *pme,
                         }
 
                         loop_count =
-                            solve_pme_lj_yzx(pme, &pme->cfftgrid[2], TRUE, ewaldcoeff_lj,
+                            solve_pme_lj_yzx(pme, &pme->cfftgrid[2], TRUE,
                                              box[XX][XX]*box[YY][YY]*box[ZZ][ZZ],
                                              bCalcEnerVir,
                                              pme->nthread, thread);

diff --git a/src/gromacs/ewald/pme.h b/src/gromacs/ewald/pme.h
index 11d843070e7979b50fda2836da3ef5a95ee4fc83..2712cc7e56bdbdc6fc8f381897c14642e6b55220 100644 (file)
--- a/src/gromacs/ewald/pme.h
+++ b/src/gromacs/ewald/pme.h
@@ -74,7 +74,9 @@ int gmx_pme_init(struct gmx_pme_t **pmedata, struct t_commrec *cr,
                  int nnodes_major, int nnodes_minor,
                  t_inputrec *ir, int homenr,
                  gmx_bool bFreeEnergy_q, gmx_bool bFreeEnergy_lj,
-                 gmx_bool bReproducible, int nthread);
+                 gmx_bool bReproducible,
+                 real ewaldcoeff_q, real ewaldcoeff_lj,
+                 int nthread);
 
 /*! \brief Destroy the PME data structures respectively.
  *
@@ -113,8 +115,7 @@ int gmx_pme_do(struct gmx_pme_t *pme,
                matrix box,      t_commrec *cr,
                int  maxshift_x, int maxshift_y,
                t_nrnb *nrnb,    gmx_wallcycle_t wcycle,
-               matrix vir_q,    real ewaldcoeff_q,
-               matrix vir_lj,   real ewaldcoeff_lj,
+               matrix vir_q,    matrix vir_lj,
                real *energy_q,  real *energy_lj,
                real lambda_q,   real lambda_lj,
                real *dvdlambda_q, real *dvdlambda_lj,

diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp
index c55a0d967538fde60fc80554768178320497f126..cf842a9a0c1c99e3f7990ea27a1abd549d9cfd69 100644 (file)
--- a/src/gromacs/mdlib/force.cpp
+++ b/src/gromacs/mdlib/force.cpp
@@ -521,8 +521,7 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
                                         DOMAINDECOMP(cr) ? dd_pme_maxshift_x(cr->dd) : 0,
                                         DOMAINDECOMP(cr) ? dd_pme_maxshift_y(cr->dd) : 0,
                                         nrnb, wcycle,
-                                        fr->vir_el_recip, fr->ewaldcoeff_q,
-                                        fr->vir_lj_recip, fr->ewaldcoeff_lj,
+                                        fr->vir_el_recip, fr->vir_lj_recip,
                                         &Vlr_q, &Vlr_lj,
                                         lambda[efptCOUL], lambda[efptVDW],
                                         &dvdl_long_range_q, &dvdl_long_range_lj, pme_flags);

diff --git a/src/programs/mdrun/runner.cpp b/src/programs/mdrun/runner.cpp
index 638561f1ba30a7abb619bf93503bb2b6cf933785..8fa892bfffd3764dbc82bab1f4797711a0f4832d 100644 (file)
--- a/src/programs/mdrun/runner.cpp
+++ b/src/programs/mdrun/runner.cpp
@@ -1311,7 +1311,9 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
         {
             status = gmx_pme_init(pmedata, cr, npme_major, npme_minor, inputrec,
                                   mtop ? mtop->natoms : 0, nChargePerturbed, nTypePerturbed,
-                                  (Flags & MD_REPRODUCIBLE), nthreads_pme);
+                                  (Flags & MD_REPRODUCIBLE),
+                                  ewaldcoeff_q, ewaldcoeff_lj,
+                                  nthreads_pme);
             if (status != 0)
             {
                 gmx_fatal(FARGS, "Error %d initializing PME", status);