Reworked convert-tpr runtime extension
authorEliane Briand <eliane.briand@mpibpc.mpg.de>
Fri, 28 May 2021 07:11:10 +0000 (07:11 +0000)
committerPaul Bauer <paul.bauer.q@gmail.com>
Fri, 28 May 2021 07:11:10 +0000 (07:11 +0000)
src/gromacs/tools/convert_tpr.cpp
src/gromacs/tools/tests/CMakeLists.txt
src/gromacs/tools/tests/convert-tpr.cpp [new file with mode: 0644]
src/gromacs/tools/tests/refdata/HelpwritingTest_ConvertTprWritesHelp.xml

index a1dd5ef858d9421a19f75a3e626e59274b1bf49b..9d93a0d082ac7cc229304ac22ff7c173e87d3bc8 100644 (file)
@@ -283,6 +283,19 @@ static void zeroq(const int index[], gmx_mtop_t* mtop)
     }
 }
 
+static void print_runtime_info(t_inputrec* ir)
+{
+    char buf[STEPSTRSIZE];
+
+    printf("  Run start step                %22s     \n", gmx_step_str(ir->init_step, buf));
+    printf("  Run start time                %22g ps  \n", ir->init_step * ir->delta_t + ir->init_t);
+    printf("  Step to be made during run    %22s     \n", gmx_step_str(ir->nsteps, buf));
+    printf("  Runtime for the run           %22g ps  \n", ir->nsteps * ir->delta_t);
+    printf("  Run end step                  %22s     \n", gmx_step_str(ir->init_step + ir->nsteps, buf));
+    printf("  Run end time                  %22g ps  \n\n",
+           (ir->init_step + ir->nsteps) * ir->delta_t + ir->init_t);
+}
+
 namespace gmx
 {
 
@@ -376,8 +389,10 @@ void ConvertTpr::initOptions(IOptionsContainer* options, ICommandLineOptionsModu
                                .storeIsSet(&runToMaxTimeIsSet_)
                                .timeValue()
                                .description("Extend runtime until this ending time (ps)"));
-    options->addOption(
-            Int64Option("nsteps").store(&maxSteps_).storeIsSet(&maxStepsIsSet_).description("Change the number of steps"));
+    options->addOption(Int64Option("nsteps")
+                               .store(&maxSteps_)
+                               .storeIsSet(&maxStepsIsSet_)
+                               .description("Change the number of steps remaining to be made"));
     options->addOption(
             BooleanOption("zeroq").store(&zeroQIsSet_).description("Set the charges of a group (from the index) to zero"));
 }
@@ -389,119 +404,124 @@ int ConvertTpr::run()
     gmx_mtop_t mtop;
     t_atoms    atoms;
     t_state    state;
-    char       buf[200], buf2[200];
+    char       buf[STEPSTRSIZE];
+
+    if (static_cast<int>(extendTimeIsSet_) + static_cast<int>(maxStepsIsSet_) + static_cast<int>(runToMaxTimeIsSet_)
+        > 1)
+    {
+        printf("Multiple runtime modification operations cannot be done in a single call.\n");
+        return 1;
+    }
+
+    if ((extendTimeIsSet_ || maxStepsIsSet_ || runToMaxTimeIsSet_) && (zeroQIsSet_ || haveReadIndexFile_))
+    {
+        printf("Cannot do both runtime modification and charge-zeroing/index group extraction in a "
+               "single call.\n");
+        return 1;
+    }
+
+    if (zeroQIsSet_ && !haveReadIndexFile_)
+    {
+        printf("Charge zeroing need an index file.\n");
+        return 1;
+    }
 
-    fprintf(stderr, "Reading toplogy and stuff from %s\n", inputTprFileName_.c_str());
 
     t_inputrec  irInstance;
     t_inputrec* ir = &irInstance;
     read_tpx_state(inputTprFileName_.c_str(), ir, &state, &mtop);
-    int64_t currentMaxStep    = ir->init_step;
-    double  currentRunTime    = ir->init_step * ir->delta_t + ir->init_t;
-    real    currentMaxRunTime = 0.0;
 
-    if (maxStepsIsSet_)
-    {
-        fprintf(stderr, "Setting nsteps to %s\n", gmx_step_str(maxSteps_, buf));
-        ir->nsteps = maxSteps_;
-    }
-    else
+
+    if (extendTimeIsSet_ || maxStepsIsSet_ || runToMaxTimeIsSet_)
     {
-        /* Determine total number of steps remaining */
-        if (extendTimeIsSet_)
-        {
-            ir->nsteps = ir->nsteps - (currentMaxStep - ir->init_step)
-                         + gmx::roundToInt64(extendTime_ / ir->delta_t);
-            printf("Extending remaining runtime of by %g ps (now %s steps)\n",
-                   extendTime_,
-                   gmx_step_str(ir->nsteps, buf));
-        }
-        else if (runToMaxTimeIsSet_)
+        // Doing runtime modification
+
+        const double  inputTimeAtStartOfRun = ir->init_step * ir->delta_t + ir->init_t;
+        const int64_t inputStepAtEndOfRun   = ir->init_step + ir->nsteps;
+        const double  inputTimeAtEndOfRun   = inputStepAtEndOfRun * ir->delta_t + ir->init_t;
+
+        printf("Input file:\n");
+        print_runtime_info(ir);
+
+        if (maxStepsIsSet_)
         {
-            printf("nsteps = %s, run_step = %s, current_t = %g, until = %g\n",
-                   gmx_step_str(ir->nsteps, buf),
-                   gmx_step_str(currentMaxStep, buf2),
-                   currentRunTime,
-                   runToMaxTime_);
-            ir->nsteps = gmx::roundToInt64((currentMaxRunTime - currentRunTime) / ir->delta_t);
-            printf("Extending remaining runtime until %g ps (now %s steps)\n",
-                   currentMaxRunTime,
-                   gmx_step_str(ir->nsteps, buf));
+            fprintf(stderr, "Setting nsteps to %s\n", gmx_step_str(maxSteps_, buf));
+            ir->nsteps = maxSteps_;
         }
-        else
+        else if (extendTimeIsSet_)
         {
-            ir->nsteps -= currentMaxStep - ir->init_step;
-            /* Print message */
-            printf("%s steps (%g ps) remaining from first run.\n",
-                   gmx_step_str(ir->nsteps, buf),
-                   ir->nsteps * ir->delta_t);
+            printf("Extending remaining runtime by %g ps\n", extendTime_);
+            ir->nsteps += gmx::roundToInt64(extendTime_ / ir->delta_t);
         }
-    }
-
-    if (maxStepsIsSet_ || zeroQIsSet_ || (ir->nsteps > 0))
-    {
-        ir->init_step = currentMaxStep;
-
-        if (haveReadIndexFile_ || !(maxStepsIsSet_ || extendTimeIsSet_ || runToMaxTimeIsSet_))
+        else if (runToMaxTimeIsSet_)
         {
-            atoms         = gmx_mtop_global_atoms(mtop);
-            int   gnx     = 0;
-            int*  index   = nullptr;
-            char* grpname = nullptr;
-            get_index(&atoms, inputIndexFileName_.c_str(), 1, &gnx, &index, &grpname);
-            bool bSel = false;
-            if (!zeroQIsSet_)
+            if (runToMaxTime_ <= inputTimeAtStartOfRun)
             {
-                bSel = (gnx != state.natoms);
-                for (int i = 0; ((i < gnx) && (!bSel)); i++)
-                {
-                    bSel = (i != index[i]);
-                }
+                printf("The requested run end time is at/before the run start time.\n");
+                return 1;
             }
-            else
+            if (runToMaxTime_ < inputTimeAtEndOfRun)
             {
-                bSel = false;
-            }
-            if (bSel)
-            {
-                fprintf(stderr,
-                        "Will write subset %s of original tpx containing %d "
-                        "atoms\n",
-                        grpname,
-                        gnx);
-                reduce_topology_x(gnx, index, &mtop, state.x.rvec_array(), state.v.rvec_array());
-                state.natoms = gnx;
-            }
-            else if (zeroQIsSet_)
-            {
-                zeroq(index, &mtop);
-                fprintf(stderr, "Zero-ing charges for group %s\n", grpname);
+                printf("The requested run end time is before the original run end time.\n");
+                printf("Reducing remaining runtime to %g ps\n", runToMaxTime_);
             }
             else
             {
-                fprintf(stderr, "Will write full tpx file (no selection)\n");
+                printf("Extending remaining runtime to %g ps\n", runToMaxTime_);
             }
+            ir->nsteps = gmx::roundToInt64((runToMaxTime_ - inputTimeAtStartOfRun) / ir->delta_t);
         }
 
-        double stateTime = ir->init_t + ir->init_step * ir->delta_t;
-        sprintf(buf,
-                "Writing statusfile with starting step %s%s and length %s%s steps...\n",
-                "%10",
-                PRId64,
-                "%10",
-                PRId64);
-        fprintf(stderr, buf, ir->init_step, ir->nsteps);
-        fprintf(stderr,
-                "                                 time %10.3f and length %10.3f ps\n",
-                stateTime,
-                ir->nsteps * ir->delta_t);
-        write_tpx_state(outputTprFileName_.c_str(), ir, &state, mtop);
+        printf("\nOutput file:\n");
+        print_runtime_info(ir);
     }
     else
     {
-        printf("You've simulated long enough. Not writing tpr file\n");
+        // If zeroQIsSet_, then we are doing charge zero-ing; otherwise index group extraction
+        // In both cases an index filename has been provided
+
+        atoms         = gmx_mtop_global_atoms(mtop);
+        int   gnx     = 0;
+        int*  index   = nullptr;
+        char* grpname = nullptr;
+        get_index(&atoms, inputIndexFileName_.c_str(), 1, &gnx, &index, &grpname);
+        bool bSel = false;
+        if (!zeroQIsSet_)
+        {
+            bSel = (gnx != state.natoms);
+            for (int i = 0; ((i < gnx) && (!bSel)); i++)
+            {
+                bSel = (i != index[i]);
+            }
+        }
+        else
+        {
+            bSel = false;
+        }
+        if (bSel)
+        {
+            fprintf(stderr,
+                    "Will write subset %s of original tpx containing %d "
+                    "atoms\n",
+                    grpname,
+                    gnx);
+            reduce_topology_x(gnx, index, &mtop, state.x.rvec_array(), state.v.rvec_array());
+            state.natoms = gnx;
+        }
+        else if (zeroQIsSet_)
+        {
+            zeroq(index, &mtop);
+            fprintf(stderr, "Zero-ing charges for group %s\n", grpname);
+        }
+        else
+        {
+            fprintf(stderr, "Will write full tpx file (no selection)\n");
+        }
     }
 
+    // Writing output tpx regardless of the operation performed
+    write_tpx_state(outputTprFileName_.c_str(), ir, &state, mtop);
+
     return 0;
 }
 
index 22f15ad5853e93c25d5d67da986651619f5f4ab2..4066c68aa8cde703edc888bf1f0b5382bc6d54c8 100644 (file)
@@ -1,7 +1,7 @@
 #
 # This file is part of the GROMACS molecular simulation package.
 #
-# Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+# Copyright (c) 2018,2019,2020,2021, 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.
@@ -38,5 +38,6 @@ gmx_add_gtest_executable(tool-test
         helpwriting.cpp
         report_methods.cpp
         trjconv.cpp
+        convert-tpr.cpp
         )
 gmx_register_gtest_test(ToolUnitTests tool-test SLOW_TEST)
diff --git a/src/gromacs/tools/tests/convert-tpr.cpp b/src/gromacs/tools/tests/convert-tpr.cpp
new file mode 100644 (file)
index 0000000..ef15481
--- /dev/null
@@ -0,0 +1,202 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, 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 gmx convert-tpr.
+ *
+ * \author Eliane Briand <eliane@br.iand.fr>
+ */
+#include "gmxpre.h"
+
+#include "gromacs/math/functions.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/state.h"
+#include "gromacs/tools/convert_tpr.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/fileio/tpxio.h"
+
+#include "testutils/cmdlinetest.h"
+#include "testutils/simulationdatabase.h"
+#include "testutils/tprfilegenerator.h"
+
+namespace gmx
+{
+namespace test
+{
+namespace
+{
+
+class ConvertTprTest : public ::testing::Test
+{
+protected:
+    ConvertTprTest() : tprFileHandle("lysozyme") {}
+
+    //! Storage for opened file handles.
+    TprAndFileManager tprFileHandle;
+};
+
+
+void readTprInput(const char* filename, gmx_mtop_t* mtop, t_inputrec* ir)
+{
+    // read tpr into variables needed for output
+    {
+        t_state state;
+        read_tpx_state(filename, ir, &state, mtop);
+    }
+}
+
+TEST_F(ConvertTprTest, ExtendRuntimeExtensionTest)
+{
+    gmx_mtop_t top;
+    t_inputrec ir;
+    readTprInput(tprFileHandle.tprName().c_str(), &top, &ir);
+
+    const int64_t originalNStep = ir.nsteps;
+
+    const int64_t extendByPs     = 100;
+    std::string   extendByString = std::to_string(extendByPs);
+
+    TestFileManager   fileManager;
+    std::string       outTprFilename = fileManager.getTemporaryFilePath("extended.tpr");
+    const char* const command[]      = {
+        "convert-tpr",          "-s",      tprFileHandle.tprName().c_str(), "-o",
+        outTprFilename.c_str(), "-extend", extendByString.c_str()
+    };
+    CommandLine cmdline(command);
+
+    gmx::test::CommandLineTestHelper::runModuleFactory(&gmx::ConvertTprInfo::create, &cmdline);
+
+    {
+        gmx_mtop_t top_after;
+        t_inputrec ir_after;
+        readTprInput(outTprFilename.c_str(), &top_after, &ir_after);
+
+        EXPECT_EQ(ir_after.nsteps, originalNStep + gmx::roundToInt64(extendByPs / ir.delta_t));
+    }
+
+
+    // Extending again (tests nsteps not zero initially
+
+    std::string anotherOutTprFilename = fileManager.getTemporaryFilePath("extended_again.tpr");
+
+    const char* const secondCommand[] = { "convert-tpr",
+                                          "-s",
+                                          outTprFilename.c_str(),
+                                          "-o",
+                                          anotherOutTprFilename.c_str(),
+                                          "-extend",
+                                          extendByString.c_str() };
+    CommandLine       secondCmdline(secondCommand);
+    gmx::test::CommandLineTestHelper::runModuleFactory(&gmx::ConvertTprInfo::create, &secondCmdline);
+
+
+    {
+        gmx_mtop_t top_after;
+        t_inputrec ir_after;
+        readTprInput(anotherOutTprFilename.c_str(), &top_after, &ir_after);
+
+        EXPECT_EQ(ir_after.nsteps, originalNStep + gmx::roundToInt64(2 * extendByPs / ir.delta_t));
+    }
+}
+
+TEST_F(ConvertTprTest, UntilRuntimeExtensionTest)
+{
+    gmx_mtop_t top;
+    t_inputrec ir;
+    readTprInput(tprFileHandle.tprName().c_str(), &top, &ir);
+
+    const int64_t originalNStep = ir.nsteps;
+
+    const int64_t untilPs       = 100;
+    std::string   untilPsString = std::to_string(untilPs);
+
+    TestFileManager   fileManager;
+    std::string       outTprFilename = fileManager.getTemporaryFilePath("extended.tpr");
+    const char* const command[]      = { "convert-tpr",
+                                    "-s",
+                                    tprFileHandle.tprName().c_str(),
+                                    "-o",
+                                    outTprFilename.c_str(),
+                                    "-until",
+                                    untilPsString.data() };
+    CommandLine       cmdline(command);
+
+    gmx::test::CommandLineTestHelper::runModuleFactory(&gmx::ConvertTprInfo::create, &cmdline);
+
+    {
+        gmx_mtop_t top_after;
+        t_inputrec ir_after;
+        readTprInput(outTprFilename.c_str(), &top_after, &ir_after);
+
+        EXPECT_EQ(ir_after.nsteps, originalNStep + gmx::roundToInt64(untilPs / ir.delta_t));
+    }
+}
+
+TEST_F(ConvertTprTest, nstepRuntimeExtensionTest)
+{
+    gmx_mtop_t top;
+    t_inputrec ir;
+    readTprInput(tprFileHandle.tprName().c_str(), &top, &ir);
+
+    const int64_t originalNStep = ir.nsteps;
+
+    const int64_t nsteps    = 102;
+    std::string   nstepsStr = std::to_string(nsteps);
+
+    TestFileManager   fileManager;
+    std::string       outTprFilename = fileManager.getTemporaryFilePath("extended.tpr");
+    const char* const command[]      = { "convert-tpr",
+                                    "-s",
+                                    tprFileHandle.tprName().c_str(),
+                                    "-o",
+                                    outTprFilename.c_str(),
+                                    "-nsteps",
+                                    nstepsStr.data() };
+    CommandLine       cmdline(command);
+
+    gmx::test::CommandLineTestHelper::runModuleFactory(&gmx::ConvertTprInfo::create, &cmdline);
+
+    {
+        gmx_mtop_t top_after;
+        t_inputrec ir_after;
+        readTprInput(outTprFilename.c_str(), &top_after, &ir_after);
+
+        EXPECT_EQ(ir_after.nsteps, originalNStep + nsteps);
+    }
+}
+
+} // namespace
+} // namespace test
+} // namespace gmx
index b0c890552b3dadd24fb46a505ee7c803a542fb07..25765ad1f390c6e36e6dad97822db73344b08e45 100644 (file)
@@ -44,7 +44,7 @@ Other options:
  -until  <time>             (0)
            Extend runtime until this ending time (ps)
  -nsteps <int>              (0)
-           Change the number of steps
+           Change the number of steps remaining to be made
  -[no]zeroq                 (no)
            Set the charges of a group (from the index) to zero
 ]]></String>