Check for implicit solvent + Verlet scheme

[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / readir.c

diff --git a/src/gromacs/gmxpreprocess/readir.c b/src/gromacs/gmxpreprocess/readir.c
index 1cfdf9b690646d2da0a039fe3eaaf7f8815fb5a1..763b6d8e5da1001ebcde6dc9a9b4af91f2971788 100644 (file)
--- a/src/gromacs/gmxpreprocess/readir.c
+++ b/src/gromacs/gmxpreprocess/readir.c
@@ -370,6 +370,11 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
             warning_error(wi, warn_buf);
         }
 
+        if (ir->implicit_solvent != eisNO)
+        {
+            warning_error(wi, "Implicit solvent is not (yet) supported with the with Verlet lists.");
+        }
+
         if (ir->nstlist <= 0)
         {
             warning_error(wi, "With Verlet lists nstlist should be larger than 0");
@@ -575,6 +580,8 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
         CHECK(ir->nstlist <= 0);
         sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
         CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
+        sprintf(err_buf, "TPI does not work (yet) with the Verlet cut-off scheme");
+        CHECK(ir->cutoff_scheme == ecutsVERLET);
     }
 
     /* SHAKE / LINCS */
@@ -1003,7 +1010,7 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
         sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
         CHECK(ir->tau_p <= 0);
 
-        if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
+        if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
         {
             sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
                     EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
@@ -1329,11 +1336,9 @@ nd %s",
     /* IMPLICIT SOLVENT */
     if (ir->coulombtype == eelGB_NOTUSED)
     {
-        ir->coulombtype      = eelCUT;
-        ir->implicit_solvent = eisGBSA;
-        fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
-                "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
-                "setting implicit-solvent value to \"GBSA\" in input section.\n");
+        sprintf(warn_buf, "Invalid option %s for coulombtype",
+                eel_names[ir->coulombtype]);
+        warning_error(wi, warn_buf);
     }
 
     if (ir->sa_algorithm == esaSTILL)
@@ -2139,7 +2144,7 @@ void get_ir(const char *mdparin, const char *mdparout,
     STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], NULL);
     ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
     STYPE ("init-lambda-weights", is->lambda_weights, NULL);
-    EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
+    EETYPE("dhdl-print-energy", fep->edHdLPrintEnergy, edHdLPrintEnergy_names);
     RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
     ITYPE ("sc-power", fep->sc_power, 1);
     RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
@@ -2364,11 +2369,23 @@ void get_ir(const char *mdparin, const char *mdparout,
         }
     }
 
-    if (ir->bSimTemp)
+    if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
     {
-        fep->bPrintEnergy = TRUE;
-        /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
-           if the temperature is changing. */
+        fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
+        warning_note(wi, "Old option for dhdl-print-energy given: "
+                     "changing \"yes\" to \"total\"\n");
+    }
+
+    if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
+    {
+        /* always print out the energy to dhdl if we are doing
+           expanded ensemble, since we need the total energy for
+           analysis if the temperature is changing. In some
+           conditions one may only want the potential energy, so
+           we will allow that if the appropriate mdp setting has
+           been enabled. Otherwise, total it is:
+         */
+        fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
     }
 
     if ((ir->efep != efepNO) || ir->bSimTemp)
@@ -3266,7 +3283,7 @@ void do_index(const char* mdparin, const char *ndx,
         nstcmin = tcouple_min_integration_steps(ir->etc);
         if (nstcmin > 1)
         {
-            if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
+            if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
             {
                 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
                         ETCOUPLTYPE(ir->etc),
@@ -4116,11 +4133,12 @@ void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
         CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
     }
 
-    if (ir->eI == eiSD1 &&
-        (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
-         gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
+    if (ir->eI == eiSD2)
     {
-        sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
+        sprintf(warn_buf, "The stochastic dynamics integrator %s is deprecated, since\n"
+                "it is slower than integrator %s and is slightly less accurate\n"
+                "with constraints. Use the %s integrator.",
+                ei_names[ir->eI], ei_names[eiSD1], ei_names[eiSD1]);
         warning_note(wi, warn_buf);
     }