From: Berk Hess <hess@kth.se>
Date: Mon, 25 May 2015 13:54:19 +0000 (+0200)
Subject: Corrected nstcalclr/nstfep/nst_repl_ex checks
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=79703ee4eff93e10c5e3308ff82fc6eb659e7977;p=alexxy%2Fgromacs.git

Corrected nstcalclr/nstfep/nst_repl_ex checks

Several checks with twin-range electrostatics still checked for
multiples of nstlist instead of nstcalclr and nst_repl_ex was not
checked at all.
Free-energy with small nstdhdl triggered frequent, unnecessary
neighbor search. It seems this can't have caused isssues with
free-energy, because nstfep needs to be a multiple of nstcalcenergy
which again needs to be a multiple of nstcalclr (was nstlist).

Change-Id: I20b8e8ab8329a3aaf904ddfa21573aee8d58c0e4
---

diff --git a/src/gromacs/gmxpreprocess/readir.c b/src/gromacs/gmxpreprocess/readir.c
index cb01302bd2..52f408f21d 100644
--- a/src/gromacs/gmxpreprocess/readir.c
+++ b/src/gromacs/gmxpreprocess/readir.c
@@ -530,7 +530,7 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
         }
         if (IR_TWINRANGE(*ir))
         {
-            check_nst("nstlist", ir->nstlist,
+            check_nst("nstcalclr", ir->nstcalclr,
                       "nstcalcenergy", &ir->nstcalcenergy, wi);
             if (ir->epc != epcNO)
             {
@@ -757,6 +757,19 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
                 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
             }
         }
+
+        if (IR_TWINRANGE(*ir))
+        {
+            sprintf(err_buf, "nstdhdl must be divisible by nstcalclr");
+            CHECK(ir->fepvals->nstdhdl > 0 &&
+                  ir->fepvals->nstdhdl % ir->nstcalclr != 0);
+
+            if (ir->efep == efepEXPANDED)
+            {
+                sprintf(err_buf, "nstexpanded must be divisible by nstcalclr");
+                CHECK(ir->expandedvals->nstexpanded % ir->nstcalclr != 0);
+            }
+        }
     }
 
     if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
diff --git a/src/programs/mdrun/md.cpp b/src/programs/mdrun/md.cpp
index 060595cb18..505b9a7ca3 100644
--- a/src/programs/mdrun/md.cpp
+++ b/src/programs/mdrun/md.cpp
@@ -217,7 +217,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     matrix            lastbox;
     int               lamnew  = 0;
     /* for FEP */
-    int               nstfep;
+    int               nstfep = 0;
     double            cycles;
     real              saved_conserved_quantity = 0;
     real              last_ekin                = 0;
@@ -514,16 +514,37 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
     debug_gmx();
 
-    /* set free energy calculation frequency as the minimum
-       greatest common denominator of nstdhdl, nstexpanded, and repl_ex_nst*/
-    nstfep = ir->fepvals->nstdhdl;
-    if (ir->bExpanded)
+    if (IR_TWINRANGE(*ir) && repl_ex_nst % ir->nstcalclr != 0)
     {
-        nstfep = gmx_greatest_common_divisor(ir->fepvals->nstdhdl, nstfep);
+        /* We should exchange at nstcalclr steps to get correct integration */
+        gmx_fatal(FARGS, "The replica exchange period (%d) is not divisible by nstcalclr (%d)", repl_ex_nst, ir->nstcalclr);
     }
-    if (repl_ex_nst > 0)
+
+    if (ir->efep != efepNO)
     {
-        nstfep = gmx_greatest_common_divisor(repl_ex_nst, nstfep);
+        /* Set free energy calculation frequency as the greatest common
+         * denominator of nstdhdl and repl_ex_nst.
+         * Check for nstcalclr with twin-range, since we need the long-range
+         * contribution to the free-energy at the correct (nstcalclr) steps.
+         */
+        nstfep = ir->fepvals->nstdhdl;
+        if (ir->bExpanded)
+        {
+            if (IR_TWINRANGE(*ir) &&
+                ir->expandedvals->nstexpanded % ir->nstcalclr != 0)
+            {
+                gmx_fatal(FARGS, "nstexpanded should be divisible by nstcalclr");
+            }
+        }
+        if (repl_ex_nst > 0)
+        {
+            nstfep = gmx_greatest_common_divisor(repl_ex_nst, nstfep);
+        }
+        /* We checked divisibility of repl_ex_nst and nstcalclr above */
+        if (IR_TWINRANGE(*ir) && nstfep % ir->nstcalclr != 0)
+        {
+            gmx_incons("nstfep not divisible by nstcalclr");
+        }
     }
 
     /* Be REALLY careful about what flags you set here. You CANNOT assume
@@ -772,7 +793,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
             set_current_lambdas(step, ir->fepvals, bRerunMD, &rerun_fr, state_global, state, lam0);
             bDoDHDL      = do_per_step(step, ir->fepvals->nstdhdl);
-            bDoFEP       = (do_per_step(step, nstfep) && (ir->efep != efepNO));
+            bDoFEP       = ((ir->efep != efepNO) && do_per_step(step, nstfep));
             bDoExpanded  = (do_per_step(step, ir->expandedvals->nstexpanded)
                             && (ir->bExpanded) && (step > 0) && (!bStartingFromCpt));
         }
@@ -854,7 +875,8 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         }
         else
         {
-            bNS = (bFirstStep || bExchanged || bNeedRepartition || bNStList || bDoFEP);
+            /* Determine whether or not to do Neighbour Searching and LR */
+            bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
         }
 
         /* check whether we should stop because another simulation has