Avoid PME tuning decreasing rcoulomb

[alexxy/gromacs.git] / src / programs / mdrun / pme_loadbal.c

diff --git a/src/programs/mdrun/pme_loadbal.c b/src/programs/mdrun/pme_loadbal.c
index 876c7881429e1983da00aabacba2b99c40f0d4c8..d6fd68f37269d10692d625172706519c053c2cf3 100644 (file)
--- a/src/programs/mdrun/pme_loadbal.c
+++ b/src/programs/mdrun/pme_loadbal.c
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2012,2013, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014, 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.
@@ -36,7 +36,8 @@
 #include <config.h>
 #endif
 
-#include "smalloc.h"
+#include "gromacs/utility/smalloc.h"
+#include "types/commrec.h"
 #include "network.h"
 #include "calcgrid.h"
 #include "pme.h"
@@ -260,6 +261,15 @@ static gmx_bool pme_loadbal_increase_cutoff(pme_load_balancing_t  pme_lb,
     while (sp <= 1.001*pme_lb->setup[pme_lb->cur].spacing || !grid_ok);
 
     set->rcut_coulomb = pme_lb->cut_spacing*sp;
+    if (set->rcut_coulomb < pme_lb->rcut_coulomb_start)
+    {
+        /* This is unlikely, but can happen when e.g. continuing from
+         * a checkpoint after equilibration where the box shrank a lot.
+         * We want to avoid rcoulomb getting smaller than rvdw
+         * and there might be more issues with decreasing rcoulomb.
+         */
+        set->rcut_coulomb = pme_lb->rcut_coulomb_start;
+    }
 
     if (pme_lb->cutoff_scheme == ecutsVERLET)
     {
@@ -317,6 +327,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;
@@ -373,7 +386,7 @@ static int pme_loadbal_end(pme_load_balancing_t pme_lb)
 }
 
 static void print_loadbal_limited(FILE *fp_err, FILE *fp_log,
-                                  gmx_large_int_t step,
+                                  gmx_int64_t step,
                                   pme_load_balancing_t pme_lb)
 {
     char buf[STRLEN], sbuf[22];
@@ -433,7 +446,7 @@ gmx_bool pme_load_balance(pme_load_balancing_t pme_lb,
                           interaction_const_t *ic,
                           nonbonded_verlet_t  *nbv,
                           gmx_pme_t           *pmedata,
-                          gmx_large_int_t      step)
+                          gmx_int64_t          step)
 {
     gmx_bool     OK;
     pme_setup_t *set;
@@ -662,6 +675,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 &&
@@ -688,17 +722,12 @@ gmx_bool pme_load_balance(pme_load_balancing_t pme_lb,
         }
 #endif  /* GMX_THREAD_MPI */
     }
-    else
-    {
-        init_interaction_const_tables(NULL, ic, bUsesSimpleTables,
-                                      rtab);
-    }
 
-    if (pme_lb->cutoff_scheme == ecutsVERLET && nbv->ngrp > 1)
-    {
-        init_interaction_const_tables(NULL, ic, bUsesSimpleTables,
-                                      rtab);
-    }
+    /* Usually we won't need the simple tables with GPUs.
+     * But we do with hybrid acceleration and with free energy.
+     * To avoid bugs, we always re-initialize the simple tables here.
+     */
+    init_interaction_const_tables(NULL, ic, bUsesSimpleTables, rtab);
 
     if (cr->duty & DUTY_PME)
     {
@@ -810,8 +839,8 @@ static void print_pme_loadbal_settings(pme_load_balancing_t pme_lb,
     {
         md_print_warn(cr, fplog,
                       "NOTE: PME load balancing increased the non-bonded workload by more than 50%%.\n"
-                      "      For better performance use (more) PME nodes (mdrun -npme),\n"
-                      "      or in case you are beyond the scaling limit, use less nodes in total.\n");
+                      "      For better performance, use (more) PME ranks (mdrun -npme),\n"
+                      "      or if you are beyond the scaling limit, use fewer total ranks (or nodes).\n");
     }
     else
     {