Support trajectory analysis with any atom subset
authorTeemu Murtola <teemu.murtola@gmail.com>
Sat, 28 Nov 2015 19:47:46 +0000 (21:47 +0200)
committerTeemu Murtola <teemu.murtola@gmail.com>
Sun, 10 Jan 2016 16:19:44 +0000 (18:19 +0200)
Make tools written for the new C++ analysis framework support analyzing
trajectories that contain an arbitrary subset of atoms.

Part of #1861.

Change-Id: I4a658e953f6f4e3d2ec1151c9c4d405c2e888780

src/gromacs/trajectoryanalysis/cmdlinerunner.cpp
src/gromacs/trajectoryanalysis/modules/sasa.cpp
src/gromacs/trajectoryanalysis/runnercommon.cpp
src/gromacs/trajectoryanalysis/runnercommon.h
src/gromacs/trajectoryanalysis/tests/cmdlinerunner.cpp
src/gromacs/trajectoryanalysis/tests/refdata/TrajectoryAnalysisCommandLineRunnerTest_WritesHelp.xml
src/gromacs/trajectoryanalysis/tests/simple-subset.gro [new file with mode: 0644]

index bd7d76c9c1fe6471d32e7ae015446f11753753be..71992a21f2878372a804138a500a4787716e9ade 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2012,2013,2014,2015,2016, 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.
@@ -124,6 +124,7 @@ int RunnerModule::run()
 
     // Load first frame.
     common_.initFirstFrame();
+    common_.initFrameIndexGroup();
     module_->initAfterFirstFrame(settings_, common_.frame());
 
     t_pbc  pbc;
index 801b2951ea4771ffeb70b8b0075d8374b1a21997..1e30e4ff5184df5e3baab4ef83b73fa4ff49b2fe 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2006, The GROMACS development team.
- * Copyright (c) 2008,2009,2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2008,2009,2010,2011,2012,2013,2014,2015,2016, 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.
@@ -954,6 +954,10 @@ Sasa::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 
     if (bConnolly)
     {
+        if (fr.natoms != top_->atoms.nr)
+        {
+            GMX_THROW(InconsistentInputError("Connolly plot (-q) is only supported for trajectories that contain all the atoms"));
+        }
         // This is somewhat nasty, as it modifies the atoms and symtab
         // structures.  But since it is only used in the first frame, and no
         // one else uses the topology after initialization, it may just work
index e757487e125586eb0dd09e9f8a60501257295b9f..c26a50044711180cdb17af17520e4d096f7c4997 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2012,2013,2014,2015,2016, 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.
@@ -45,6 +45,9 @@
 
 #include <string.h>
 
+#include <algorithm>
+#include <string>
+
 #include "gromacs/fileio/confio.h"
 #include "gromacs/fileio/oenv.h"
 #include "gromacs/fileio/timecontrol.h"
@@ -54,7 +57,9 @@
 #include "gromacs/options/filenameoption.h"
 #include "gromacs/options/ioptionscontainer.h"
 #include "gromacs/pbcutil/rmpbc.h"
+#include "gromacs/selection/selection.h"
 #include "gromacs/selection/selectioncollection.h"
+#include "gromacs/selection/selectionoption.h"
 #include "gromacs/selection/selectionoptionbehavior.h"
 #include "gromacs/topology/topology.h"
 #include "gromacs/trajectory/trajectoryframe.h"
@@ -81,6 +86,7 @@ class TrajectoryAnalysisRunnerCommon::Impl : public ITopologyProvider
 
         void initTopology(bool required);
         void initFirstFrame();
+        void initFrameIndexGroup();
         void finishTrajectory();
 
         // From ITopologyProvider
@@ -93,6 +99,10 @@ class TrajectoryAnalysisRunnerCommon::Impl : public ITopologyProvider
         {
             if (!topInfo_.hasTopology())
             {
+                if (trajectoryGroup_.isValid())
+                {
+                    GMX_THROW(InconsistentInputError("-fgroup is only supported when -s is also specified"));
+                }
                 // Read the first frame if we don't know the maximum number of
                 // atoms otherwise.
                 initFirstFrame();
@@ -108,6 +118,7 @@ class TrajectoryAnalysisRunnerCommon::Impl : public ITopologyProvider
         std::string                 trjfile_;
         //! Name of the topology file (empty if no topology provided).
         std::string                 topfile_;
+        Selection                   trajectoryGroup_;
         double                      startTime_;
         double                      endTime_;
         double                      deltaTime_;
@@ -137,15 +148,16 @@ TrajectoryAnalysisRunnerCommon::Impl::Impl(TrajectoryAnalysisSettings *settings)
 TrajectoryAnalysisRunnerCommon::Impl::~Impl()
 {
     finishTrajectory();
-    if (fr)
+    if (fr != nullptr)
     {
         // There doesn't seem to be a function for freeing frame data
         sfree(fr->x);
         sfree(fr->v);
         sfree(fr->f);
+        sfree(fr->index);
         sfree(fr);
     }
-    if (oenv_ != NULL)
+    if (oenv_ != nullptr)
     {
         output_env_done(oenv_);
     }
@@ -247,6 +259,30 @@ TrajectoryAnalysisRunnerCommon::Impl::initFirstFrame()
     }
 }
 
+void
+TrajectoryAnalysisRunnerCommon::Impl::initFrameIndexGroup()
+{
+    if (!trajectoryGroup_.isValid())
+    {
+        return;
+    }
+    GMX_RELEASE_ASSERT(bTrajOpen_,
+                       "Trajectory index only makes sense with a real trajectory");
+    if (trajectoryGroup_.atomCount() != fr->natoms)
+    {
+        const std::string message = formatString(
+                    "Selection specified with -fgroup has %d atoms, but "
+                    "the trajectory (-f) has %d atoms.",
+                    trajectoryGroup_.atomCount(), fr->natoms);
+        GMX_THROW(InconsistentInputError(message));
+    }
+    fr->bIndex = TRUE;
+    snew(fr->index, trajectoryGroup_.atomCount());
+    std::copy(trajectoryGroup_.atomIndices().begin(),
+              trajectoryGroup_.atomIndices().end(),
+              fr->index);
+}
+
 void
 TrajectoryAnalysisRunnerCommon::Impl::finishTrajectory()
 {
@@ -322,6 +358,12 @@ TrajectoryAnalysisRunnerCommon::initOptions(IOptionsContainer *options,
     timeUnitBehavior->addTimeUnitOption(options, "tu");
     timeUnitBehavior->setTimeUnitStore(&impl_->settings_.impl_->timeUnit);
 
+    options->addOption(SelectionOption("fgroup")
+                           .store(&impl_->trajectoryGroup_)
+                           .onlySortedAtoms().onlyStatic()
+                           .description("Atoms stored in the trajectory file "
+                                        "(if not set, assume first N atoms)"));
+
     // Add plot options.
     settings.impl_->plotSettings.initOptions(options);
 
@@ -342,13 +384,18 @@ TrajectoryAnalysisRunnerCommon::initOptions(IOptionsContainer *options,
 void
 TrajectoryAnalysisRunnerCommon::optionsFinished()
 {
-    impl_->settings_.impl_->plotSettings.setTimeUnit(impl_->settings_.timeUnit());
-
     if (impl_->trjfile_.empty() && impl_->topfile_.empty())
     {
         GMX_THROW(InconsistentInputError("No trajectory or topology provided, nothing to do!"));
     }
 
+    if (impl_->trajectoryGroup_.isValid() && impl_->trjfile_.empty())
+    {
+        GMX_THROW(InconsistentInputError("-fgroup only makes sense together with a trajectory (-f)"));
+    }
+
+    impl_->settings_.impl_->plotSettings.setTimeUnit(impl_->settings_.timeUnit());
+
     if (impl_->bStartTimeSet_)
     {
         setTimeValue(TBEGIN, impl_->startTime_);
@@ -380,6 +427,13 @@ TrajectoryAnalysisRunnerCommon::initFirstFrame()
 }
 
 
+void
+TrajectoryAnalysisRunnerCommon::initFrameIndexGroup()
+{
+    impl_->initFrameIndexGroup();
+}
+
+
 bool
 TrajectoryAnalysisRunnerCommon::readNextFrame()
 {
index eb7001705f2b44cfbc92ba9a40e4d313224d2be8..20d703396950bbe9b74c38970918e51fd9f5f969 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2012,2013,2014,2015,2016, 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.
@@ -97,6 +97,12 @@ class TrajectoryAnalysisRunnerCommon
          * After this call, frame() returns the first frame.
          */
         void initFirstFrame();
+        /*! \brief
+         * Initializes the index in frame() that specifies the atoms contained.
+         *
+         * Can be called after selections have been compiled.
+         */
+        void initFrameIndexGroup();
         /*! \brief
          * Reads the next frame from the trajectory.
          *
index c78d53c3acef81d07b32c58902ef7bf69ffd4e60..221651b7f488550a7a4f138861566a7478a6fb14 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2015, by the GROMACS development team, led by
+ * Copyright (c) 2015,2016, 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.
 #include "gromacs/trajectory/trajectoryframe.h"
 #include "gromacs/trajectoryanalysis/analysismodule.h"
 #include "gromacs/trajectoryanalysis/analysissettings.h"
+#include "gromacs/utility/exceptions.h"
 
 #include "testutils/cmdlinetest.h"
+#include "testutils/testasserts.h"
 
 namespace
 {
@@ -75,6 +77,30 @@ class MockModule : public gmx::TrajectoryAnalysisModule
 
 using gmx::test::CommandLine;
 
+class TrajectoryAnalysisCommandLineRunnerTest : public gmx::test::CommandLineTestBase
+{
+    public:
+        TrajectoryAnalysisCommandLineRunnerTest()
+            : mockModule_(new MockModule())
+        {
+        }
+
+        gmx::ICommandLineOptionsModulePointer createRunner()
+        {
+            return gmx::TrajectoryAnalysisCommandLineRunner::createModule(
+                    std::move(mockModule_));
+        }
+
+        void runTest(const CommandLine &args)
+        {
+            CommandLine &cmdline = commandLine();
+            cmdline.merge(args);
+            ASSERT_EQ(0, gmx::test::CommandLineTestHelper::runModuleDirect(createRunner(), &cmdline));
+        }
+
+        std::unique_ptr<MockModule> mockModule_;
+};
+
 //! Initializes options for help testing.
 void initOptions(gmx::IOptionsContainer *options, gmx::TrajectoryAnalysisSettings *settings)
 {
@@ -87,22 +113,63 @@ void initOptions(gmx::IOptionsContainer *options, gmx::TrajectoryAnalysisSetting
     options->addOption(gmx::BooleanOption("test").description("Test option"));
 }
 
-//! Test fixture for TrajectoryAnalysisCommandLineRunner tests.
-typedef gmx::test::CommandLineTestBase TrajectoryAnalysisCommandLineRunnerTest;
-
 TEST_F(TrajectoryAnalysisCommandLineRunnerTest, WritesHelp)
 {
-    std::unique_ptr<MockModule>                    mockModule(new MockModule());
     using ::testing::_;
     using ::testing::Invoke;
-    EXPECT_CALL(*mockModule, initOptions(_, _)).WillOnce(Invoke(&initOptions));
+    EXPECT_CALL(*mockModule_, initOptions(_, _)).WillOnce(Invoke(&initOptions));
 
-    gmx::ICommandLineOptionsModulePointer          runner(
-            gmx::TrajectoryAnalysisCommandLineRunner::createModule(std::move(mockModule)));
     const std::unique_ptr<gmx::ICommandLineModule> module(
             gmx::ICommandLineOptionsModule::createModule("mod", "Description",
-                                                         std::move(runner)));
+                                                         createRunner()));
     testWriteHelp(module.get());
 }
 
+TEST_F(TrajectoryAnalysisCommandLineRunnerTest, RunsWithSubsetTrajectory)
+{
+    const char *const cmdline[] = {
+        "-fgroup", "atomnr 4 5 6 10 to 14"
+    };
+
+    using ::testing::_;
+    EXPECT_CALL(*mockModule_, initOptions(_, _));
+    EXPECT_CALL(*mockModule_, initAnalysis(_, _));
+    EXPECT_CALL(*mockModule_, analyzeFrame(0, _, _, _));
+    EXPECT_CALL(*mockModule_, analyzeFrame(1, _, _, _));
+    EXPECT_CALL(*mockModule_, finishAnalysis(2));
+    EXPECT_CALL(*mockModule_, writeOutput());
+
+    setInputFile("-s", "simple.gro");
+    setInputFile("-f", "simple-subset.gro");
+    EXPECT_NO_THROW_GMX(runTest(CommandLine(cmdline)));
+}
+
+TEST_F(TrajectoryAnalysisCommandLineRunnerTest, DetectsIncorrectTrajectorySubset)
+{
+    const char *const cmdline[] = {
+        "-fgroup", "atomnr 3 to 6 10 to 14"
+    };
+
+    using ::testing::_;
+    EXPECT_CALL(*mockModule_, initOptions(_, _));
+    EXPECT_CALL(*mockModule_, initAnalysis(_, _));
+
+    setInputFile("-s", "simple.gro");
+    setInputFile("-f", "simple-subset.gro");
+    EXPECT_THROW_GMX(runTest(CommandLine(cmdline)), gmx::InconsistentInputError);
+}
+
+TEST_F(TrajectoryAnalysisCommandLineRunnerTest, FailsWithTrajectorySubsetWithoutTrajectory)
+{
+    const char *const cmdline[] = {
+        "-fgroup", "atomnr 3 to 6 10 to 14"
+    };
+
+    using ::testing::_;
+    EXPECT_CALL(*mockModule_, initOptions(_, _));
+
+    setInputFile("-s", "simple.gro");
+    EXPECT_THROW_GMX(runTest(CommandLine(cmdline)), gmx::InconsistentInputError);
+}
+
 } // namespace
index af1c1c831ef3667071ed0dd58849ec7f3f396bf7..7af475e6c772c26f0df6fdaff15a0676386705d4 100644 (file)
@@ -5,8 +5,9 @@
 SYNOPSIS
 
 test mod [-f [<.xtc/.trr/...>]] [-s [<.tpr/.gro/...>]] [-n [<.ndx>]]
-         [-b <time>] [-e <time>] [-dt <time>] [-tu <enum>] [-xvg <enum>]
-         [-[no]rmpbc] [-[no]pbc] [-sf <file>] [-selrpos <enum>] [-[no]test]
+         [-b <time>] [-e <time>] [-dt <time>] [-tu <enum>]
+         [-fgroup <selection>] [-xvg <enum>] [-[no]rmpbc] [-[no]pbc]
+         [-sf <file>] [-selrpos <enum>] [-[no]test]
 
 DESCRIPTION
 
@@ -34,6 +35,9 @@ Other options:
            Only use frame if t MOD dt == first time (ps)
  -tu     <enum>             (ps)
            Unit for time values: fs, ps, ns, us, ms, s
+ -fgroup <selection>
+           Atoms stored in the trajectory file (if not set, assume first N
+           atoms)
  -xvg    <enum>             (xmgrace)
            Plot formatting: none, xmgrace, xmgr
  -[no]rmpbc                 (yes)
diff --git a/src/gromacs/trajectoryanalysis/tests/simple-subset.gro b/src/gromacs/trajectoryanalysis/tests/simple-subset.gro
new file mode 100644 (file)
index 0000000..4446b20
--- /dev/null
@@ -0,0 +1,22 @@
+Test system
+  8
+    3RB      CB    4   1.000   4.000   0.000
+    3RB      S1    5   2.000   1.000   0.000
+    3RB      S2    6   2.000   2.000   0.000
+    5RC      CB   10   3.000   2.000   0.000
+    5RC      S1   11   3.000   3.000   0.000
+    5RC      S2   12   3.000   4.000   0.000
+    1RD      CB   13   4.000   1.000   0.000
+    1RD      S1   14   4.000   2.000   0.000
+  10.00000  10.00000  10.00000
+Test system
+  8
+    3RB      CB    4   1.000   2.000   0.000
+    3RB      S1    5   2.000   3.000   0.000
+    3RB      S2    6   2.000   4.000   0.000
+    5RC      CB   10   3.000   3.000   0.000
+    5RC      S1   11   3.000   4.000   0.000
+    5RC      S2   12   3.000   1.000   0.000
+    1RD      CB   13   4.000   2.000   0.000
+    1RD      S1   14   4.000   3.000   0.000
+  10.00000  10.00000  10.00000