Avoid FP exception with empty .mdp file

The default value for nsteps is zero, which used to lead to division
by zero in the above case. Refactored to avoid that.

Added a basic test case for grompp with an empty .mdp file (which
causes all default settings to be used). It is not optimal that the
mdrun integration test machinery is recycled for this, but
re-implementing or refactoring that to avoid the ugliness is not great
for a release branch.

Noted some TODOs. Renamed some variables for clarity.

Change-Id: I5c98c4f460f658277490b11f7fb5629d6fa84a08

diff --git a/src/gromacs/gmxpreprocess/grompp.c b/src/gromacs/gmxpreprocess/grompp.c
index 2f976b8d2eef1c8667d524a3ac287ca999f9aef2..bc98819dabf49d01d39e0648970616e02ad5723a 100644 (file)
--- a/src/gromacs/gmxpreprocess/grompp.c
+++ b/src/gromacs/gmxpreprocess/grompp.c
@@ -1913,25 +1913,29 @@ int gmx_grompp(int argc, char *argv[])
             }
             else
             {
-                /* We warn for NVE simulations with >1(.1)% drift tolerance */
-                const real drift_tol = 0.01;
-                real       ener_runtime;
-
-                /* We use 2 DOF per atom = 2kT pot+kin energy, to be on
-                 * the safe side with constraints (without constraints: 3 DOF).
+                /* We warn for NVE simulations with a drift tolerance that
+                 * might result in a 1(.1)% drift over the total run-time.
+                 * Note that we can't warn when nsteps=0, since we don't
+                 * know how many steps the user intends to run.
                  */
-                ener_runtime = 2*BOLTZ*buffer_temp/(ir->nsteps*ir->delta_t);
-
                 if (EI_MD(ir->eI) && ir->etc == etcNO && ir->nstlist > 1 &&
-                    ir->nsteps > 0 &&
-                    ir->verletbuf_tol > 1.1*drift_tol*ener_runtime)
+                    ir->nsteps > 0)
                 {
-                    sprintf(warn_buf, "You are using a Verlet buffer tolerance of %g kJ/mol/ps for an NVE simulation of length %g ps, which can give a final drift of %d%%. For conserving energy to %d%%, you might need to set verlet-buffer-tolerance to %.1e.",
-                            ir->verletbuf_tol, ir->nsteps*ir->delta_t,
-                            (int)(ir->verletbuf_tol/ener_runtime*100 + 0.5),
-                            (int)(100*drift_tol + 0.5),
-                            drift_tol*ener_runtime);
-                    warning_note(wi, warn_buf);
+                    const real driftTolerance = 0.01;
+                    /* We use 2 DOF per atom = 2kT pot+kin energy,
+                     * to be on the safe side with constraints.
+                     */
+                    const real totalEnergyDriftPerAtomPerPicosecond = 2*BOLTZ*buffer_temp/(ir->nsteps*ir->delta_t);
+
+                    if (ir->verletbuf_tol > 1.1*driftTolerance*totalEnergyDriftPerAtomPerPicosecond)
+                    {
+                        sprintf(warn_buf, "You are using a Verlet buffer tolerance of %g kJ/mol/ps for an NVE simulation of length %g ps, which can give a final drift of %d%%. For conserving energy to %d%% when using constraints, you might need to set verlet-buffer-tolerance to %.1e.",
+                                ir->verletbuf_tol, ir->nsteps*ir->delta_t,
+                                (int)(ir->verletbuf_tol/totalEnergyDriftPerAtomPerPicosecond*100 + 0.5),
+                                (int)(100*driftTolerance + 0.5),
+                                driftTolerance*totalEnergyDriftPerAtomPerPicosecond);
+                        warning_note(wi, warn_buf);
+                    }
                 }
 
                 set_verlet_buffer(sys, ir, buffer_temp, state->box, wi);

diff --git a/src/programs/mdrun/tests/CMakeLists.txt b/src/programs/mdrun/tests/CMakeLists.txt
index c5eb3730d479c25160aa7f8837100ff43dc61ddb..5e47c91afe7343291cb8f99e181af95445ba4c06 100644 (file)
--- a/src/programs/mdrun/tests/CMakeLists.txt
+++ b/src/programs/mdrun/tests/CMakeLists.txt
@@ -39,6 +39,7 @@ gmx_build_unit_test(
     ${testname}
     ${exename}
     # files with code for tests
+    grompp.cpp
     rerun.cpp
     trajectory_writing.cpp
     compressed_x_output.cpp

diff --git a/src/programs/mdrun/tests/grompp.cpp b/src/programs/mdrun/tests/grompp.cpp
new file mode 100644 (file)
index 0000000..afd9317
--- /dev/null
+++ b/src/programs/mdrun/tests/grompp.cpp
@@ -0,0 +1,78 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2015, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ * \brief
+ * Tests for basic grompp functionality
+ *
+ * \todo Refactor SimulationRunner to split off SimulationPreparer, so
+ * that integration tests of grompp can stand apart from tests of
+ * mdrun.
+ *
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \ingroup module_mdrun_integration_tests
+ */
+#include "gmxpre.h"
+
+#include <gtest/gtest.h>
+
+#include "moduletest.h"
+
+namespace
+{
+
+//! Test fixture for grompp
+class GromppTest :
+    public gmx::test::MdrunTestFixture
+{
+    public:
+        //! Execute the trajectory writing test
+        void runTest()
+        {
+            runner_.useTopGroAndNdxFromDatabase("spc-and-methanol");
+            EXPECT_EQ(0, runner_.callGrompp());
+        }
+};
+
+/* This test ensures that an empty .mdp file (ie. all default values) works. */
+TEST_F(GromppTest, EmptyMdpFileWorks)
+{
+    /* TODO Now that Verlet is the default, change the implementation
+       of useEmptyMdpFile() to do that. */
+    runner_.useStringAsMdpFile("");
+    runTest();
+}
+
+} // namespace