PME load balancing now works with LJ-PME

PME load balancing only works with LJ-PME when PME is also used
for electrostatics.

Change-Id: I7145b242fd5dd00690e6b789b9cc6f3e147d430b

diff --git a/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_data_mgmt.cu b/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_data_mgmt.cu
index dfdafbb2fd9335d21e1501cd095144bfc606c576..57d9def08ec0317bcbc6ba273a611652226dcf46 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_data_mgmt.cu
+++ b/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_data_mgmt.cu
@@ -259,18 +259,10 @@ static int pick_ewald_kernel_type(bool                   bTwinCut,
     return kernel_type;
 }
 
-
-/*! Initializes the nonbonded parameter data structure. */
-static void init_nbparam(cu_nbparam_t              *nbp,
-                         const interaction_const_t *ic,
-                         const nbnxn_atomdata_t    *nbat,
-                         const cuda_dev_info_t     *dev_info)
+/*! Copies all parameters related to the cut-off from ic to nbp */
+static void set_cutoff_parameters(cu_nbparam_t              *nbp,
+                                  const interaction_const_t *ic)
 {
-    cudaError_t stat;
-    int         ntypes, nnbfp, nnbfp_comb;
-
-    ntypes  = nbat->ntype;
-
     nbp->ewald_beta       = ic->ewaldcoeff_q;
     nbp->sh_ewald         = ic->sh_ewald;
     nbp->epsfac           = ic->epsfac;
@@ -287,6 +279,20 @@ static void init_nbparam(cu_nbparam_t              *nbp,
     nbp->dispersion_shift = ic->dispersion_shift;
     nbp->repulsion_shift  = ic->repulsion_shift;
     nbp->vdw_switch       = ic->vdw_switch;
+}
+
+/*! Initializes the nonbonded parameter data structure. */
+static void init_nbparam(cu_nbparam_t              *nbp,
+                         const interaction_const_t *ic,
+                         const nbnxn_atomdata_t    *nbat,
+                         const cuda_dev_info_t     *dev_info)
+{
+    cudaError_t stat;
+    int         ntypes, nnbfp, nnbfp_comb;
+
+    ntypes  = nbat->ntype;
+
+    set_cutoff_parameters(nbp, ic);
 
     if (ic->vdwtype == evdwCUT)
     {
@@ -424,9 +430,7 @@ void nbnxn_cuda_pme_loadbal_update_param(nbnxn_cuda_ptr_t           cu_nb,
 {
     cu_nbparam_t *nbp = cu_nb->nbparam;
 
-    nbp->rlist_sq       = ic->rlist * ic->rlist;
-    nbp->rcoulomb_sq    = ic->rcoulomb * ic->rcoulomb;
-    nbp->ewald_beta     = ic->ewaldcoeff_q;
+    set_cutoff_parameters(nbp, ic);
 
     nbp->eeltype        = pick_ewald_kernel_type(ic->rcoulomb != ic->rvdw,
                                                  cu_nb->dev_info);

diff --git a/src/programs/mdrun/md.c b/src/programs/mdrun/md.c
index f5c151202c50b0add3cc2c3990d6ae9144f00114..04c48de3fa58b55eef45b1c3c5750cf28b3ef101 100644 (file)
--- a/src/programs/mdrun/md.c
+++ b/src/programs/mdrun/md.c
@@ -458,11 +458,13 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                                         repl_ex_nst, repl_ex_nex, repl_ex_seed);
     }
 
-    /* PME tuning is only supported with GPUs or PME nodes and not with rerun or LJ-PME. */
+    /* PME tuning is only supported with GPUs or PME nodes and not with rerun.
+     * PME tuning is not supported with PME only for LJ and not for Coulomb.
+     */
     if ((Flags & MD_TUNEPME) &&
         EEL_PME(fr->eeltype) &&
         ( (fr->cutoff_scheme == ecutsVERLET && fr->nbv->bUseGPU) || !(cr->duty & DUTY_PME)) &&
-        !bRerunMD && !EVDW_PME(fr->vdwtype))
+        !bRerunMD)
     {
         pme_loadbal_init(&pme_loadbal, ir, state->box, fr->ic, fr->pmedata);
         cycles_pmes = 0;
@@ -1878,11 +1880,17 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                                          step);
 
                     /* Update constants in forcerec/inputrec to keep them in sync with fr->ic */
-                    fr->ewaldcoeff_q = fr->ic->ewaldcoeff_q;
-                    fr->rlist        = fr->ic->rlist;
-                    fr->rlistlong    = fr->ic->rlistlong;
-                    fr->rcoulomb     = fr->ic->rcoulomb;
-                    fr->rvdw         = fr->ic->rvdw;
+                    fr->ewaldcoeff_q  = fr->ic->ewaldcoeff_q;
+                    fr->ewaldcoeff_lj = fr->ic->ewaldcoeff_lj;
+                    fr->rlist         = fr->ic->rlist;
+                    fr->rlistlong     = fr->ic->rlistlong;
+                    fr->rcoulomb      = fr->ic->rcoulomb;
+                    fr->rvdw          = fr->ic->rvdw;
+
+                    if (ir->eDispCorr != edispcNO)
+                    {
+                        calc_enervirdiff(NULL, ir->eDispCorr, fr);
+                    }
                 }
                 cycles_pmes = 0;
             }

diff --git a/src/programs/mdrun/pme_loadbal.c b/src/programs/mdrun/pme_loadbal.c
index 5ea760ec3ae69fd144b3c795f59a9209e60885db..63efb07eea468fddd884f05295c33e71daf4fc58 100644 (file)
--- a/src/programs/mdrun/pme_loadbal.c
+++ b/src/programs/mdrun/pme_loadbal.c
@@ -317,6 +317,9 @@ static gmx_bool pme_loadbal_increase_cutoff(pme_load_balancing_t  pme_lb,
     /* The Ewald coefficient is inversly proportional to the cut-off */
     set->ewaldcoeff_q =
         pme_lb->setup[0].ewaldcoeff_q*pme_lb->setup[0].rcut_coulomb/set->rcut_coulomb;
+    /* We set ewaldcoeff_lj in set, even when LJ-PME is not used */
+    set->ewaldcoeff_lj =
+        pme_lb->setup[0].ewaldcoeff_lj*pme_lb->setup[0].rcut_coulomb/set->rcut_coulomb;
 
     set->count   = 0;
     set->cycles  = 0;
@@ -662,10 +665,27 @@ gmx_bool pme_load_balance(pme_load_balancing_t pme_lb,
     ic->rlistlong    = set->rlistlong;
     ir->nstcalclr    = set->nstcalclr;
     ic->ewaldcoeff_q = set->ewaldcoeff_q;
+    /* TODO: centralize the code that sets the potentials shifts */
     if (ic->coulomb_modifier == eintmodPOTSHIFT)
     {
         ic->sh_ewald = gmx_erfc(ic->ewaldcoeff_q*ic->rcoulomb);
     }
+    if (EVDW_PME(ic->vdwtype))
+    {
+        /* We have PME for both Coulomb and VdW, set rvdw equal to rcoulomb */
+        ic->rvdw            = set->rcut_coulomb;
+        ic->ewaldcoeff_lj   = set->ewaldcoeff_lj;
+        if (ic->vdw_modifier == eintmodPOTSHIFT)
+        {
+            real crc2;
+
+            ic->dispersion_shift.cpot = -pow(ic->rvdw, -6.0);
+            ic->repulsion_shift.cpot  = -pow(ic->rvdw, -12.0);
+            ic->sh_invrc6             = -ic->dispersion_shift.cpot;
+            crc2                      = sqr(ic->ewaldcoeff_lj*ic->rvdw);
+            ic->sh_lj_ewald           = (exp(-crc2)*(1 + crc2 + 0.5*crc2*crc2) - 1)*pow(ic->rvdw, -6.0);
+        }
+    }
 
     bUsesSimpleTables = uses_simple_tables(ir->cutoff_scheme, nbv, 0);
     if (pme_lb->cutoff_scheme == ecutsVERLET &&