Fixes bug in gmx solvate with -shell that yielded 0 SOL.
authorDavid van der Spoel <spoel@xray.bmc.uu.se>
Fri, 10 Feb 2017 06:43:44 +0000 (07:43 +0100)
committerDavid van der Spoel <davidvanderspoel@gmail.com>
Mon, 27 Feb 2017 19:24:32 +0000 (20:24 +0100)
In the transition from genbox to the new solvate.cpp some
incorrect logic was introduced.

Added test for -shell.

Implemented test for conf files that check the title and
number of atoms. In principle this could be extended to
other structure files.

Fixes #2119

Change-Id: I88feef55f33c7076b9c27a9831ee04d890333b76

src/gromacs/gmxpreprocess/solvate.cpp
src/gromacs/gmxpreprocess/tests/refdata/SolvateTest_cs_box_Works.xml [new file with mode: 0644]
src/gromacs/gmxpreprocess/tests/refdata/SolvateTest_cs_cp_Works.xml [new file with mode: 0644]
src/gromacs/gmxpreprocess/tests/refdata/SolvateTest_cs_cp_p_Works.xml [new file with mode: 0644]
src/gromacs/gmxpreprocess/tests/refdata/SolvateTest_shell_Works.xml [new file with mode: 0644]
src/gromacs/gmxpreprocess/tests/solvate.cpp
src/testutils/CMakeLists.txt
src/testutils/conftest.cpp [new file with mode: 0644]
src/testutils/conftest.h [new file with mode: 0644]

index 2ee48d2fc70e2edf4bb57bd8299c4eb99e7791a2..e2c1d49ad9d28602826cfba70f393b783734536c 100644 (file)
@@ -491,50 +491,86 @@ static void removeSolventBoxOverlap(t_atoms *atoms, std::vector<RVec> *x,
 }
 
 /*! \brief
- * Removes solvent molecules that overlap with the solute, and optionally also
- * those that are outside a given shell radius from the solute.
+ * Remove all solvent molecules outside a give radius from the solute.
  *
- * \param[in,out] atoms      Solvent atoms.
- * \param[in,out] x          Solvent positions.
- * \param[in,out] v          Solvent velocities (can be empty).
- * \param[in,out] r          Solvent exclusion radii.
- * \param[in]     pbc        PBC information.
- * \param[in]     x_solute   Solute positions.
- * \param[in]     r_solute   Solute exclusion radii.
- * \param[in]     rshell     If >0, only keep solvent atoms within a shell of
- *     this size from the solute.
+ * \param[in,out] atoms     Solvent atoms.
+ * \param[in,out] x_solvent Solvent positions.
+ * \param[in,out] v_solvent Solvent velocities.
+ * \param[in,out] r         Atomic exclusion radii.
+ * \param[in]     pbc       PBC information.
+ * \param[in]     x_solute  Solute positions.
+ * \param[in]     rshell    The radius outside the solute molecule.
  */
-static void removeSoluteOverlap(t_atoms *atoms, std::vector<RVec> *x,
-                                std::vector<RVec> *v, std::vector<real> *r,
-                                const t_pbc &pbc,
-                                const std::vector<RVec> &x_solute,
-                                const std::vector<real> &r_solute,
-                                real rshell)
+static void removeSolventOutsideShell(t_atoms                 *atoms,
+                                      std::vector<RVec>       *x_solvent,
+                                      std::vector<RVec>       *v_solvent,
+                                      std::vector<real>       *r,
+                                      const t_pbc             &pbc,
+                                      const std::vector<RVec> &x_solute,
+                                      real                     rshell)
 {
-    const real                          maxRadius1
-        = *std::max_element(r->begin(), r->end());
-    const real                          maxRadius2
-        = *std::max_element(r_solute.begin(), r_solute.end());
-
     gmx::AtomsRemover                   remover(*atoms);
-    // If rshell is >0, the neighborhood search looks at all pairs
-    // within rshell, and unmarks those that are within the cutoff.
-    // This line marks everything, so that solvent outside rshell remains
-    // marked after the loop.
-    // Without rshell, the neighborhood search only marks the overlapping
-    // solvent atoms, and all others are left alone.
-    if (rshell > 0.0)
+    gmx::AnalysisNeighborhood           nb;
+    nb.setCutoff(rshell);
+    gmx::AnalysisNeighborhoodPositions  posSolute(x_solute);
+    gmx::AnalysisNeighborhoodSearch     search     = nb.initSearch(&pbc, posSolute);
+    gmx::AnalysisNeighborhoodPositions  pos(*x_solvent);
+    gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(pos);
+    gmx::AnalysisNeighborhoodPair       pair;
+
+    // Remove everything
+    remover.markAll();
+    // Now put back those within the shell without checking for overlap
+    while (pairSearch.findNextPair(&pair))
     {
-        remover.markAll();
+        remover.markResidue(*atoms, pair.testIndex(), false);
+        pairSearch.skipRemainingPairsForTestPosition();
     }
+    remover.removeMarkedElements(x_solvent);
+    if (!v_solvent->empty())
+    {
+        remover.removeMarkedElements(v_solvent);
+    }
+    remover.removeMarkedElements(r);
+    const int originalAtomCount = atoms->nr;
+    remover.removeMarkedAtoms(atoms);
+    fprintf(stderr, "Removed %d solvent atoms more than %f nm from solute.\n",
+            originalAtomCount - atoms->nr, rshell);
+}
+
+/*! \brief
+ * Removes solvent molecules that overlap with the solute.
+ *
+ * \param[in,out] atoms    Solvent atoms.
+ * \param[in,out] x        Solvent positions.
+ * \param[in,out] v        Solvent velocities (can be empty).
+ * \param[in,out] r        Solvent exclusion radii.
+ * \param[in]     pbc      PBC information.
+ * \param[in]     x_solute Solute positions.
+ * \param[in]     r_solute Solute exclusion radii.
+ */
+static void removeSolventOverlappingWithSolute(t_atoms                 *atoms,
+                                               std::vector<RVec>       *x,
+                                               std::vector<RVec>       *v,
+                                               std::vector<real>       *r,
+                                               const t_pbc             &pbc,
+                                               const std::vector<RVec> &x_solute,
+                                               const std::vector<real> &r_solute)
+{
+    gmx::AtomsRemover             remover(*atoms);
+    const real                    maxRadius1
+        = *std::max_element(r->begin(), r->end());
+    const real                    maxRadius2
+        = *std::max_element(r_solute.begin(), r_solute.end());
 
+    // Now check for overlap.
     gmx::AnalysisNeighborhood           nb;
-    nb.setCutoff(std::max(maxRadius1 + maxRadius2, rshell));
+    gmx::AnalysisNeighborhoodPair       pair;
+    nb.setCutoff(maxRadius1 + maxRadius2);
     gmx::AnalysisNeighborhoodPositions  posSolute(x_solute);
     gmx::AnalysisNeighborhoodSearch     search     = nb.initSearch(&pbc, posSolute);
     gmx::AnalysisNeighborhoodPositions  pos(*x);
     gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(pos);
-    gmx::AnalysisNeighborhoodPair       pair;
     while (pairSearch.findNextPair(&pair))
     {
         if (remover.isMarked(pair.testIndex()))
@@ -647,8 +683,14 @@ static void add_solv(const char *fn, t_topology *top,
     }
     if (top->atoms.nr > 0)
     {
-        removeSoluteOverlap(atoms_solvt, &x_solvt, &v_solvt, &exclusionDistances_solvt, pbc,
-                            *x, exclusionDistances, rshell);
+        if (rshell > 0.0)
+        {
+            removeSolventOutsideShell(atoms_solvt, &x_solvt,  &v_solvt,
+                                      &exclusionDistances_solvt, pbc, *x, rshell);
+        }
+        removeSolventOverlappingWithSolute(atoms_solvt, &x_solvt, &v_solvt,
+                                           &exclusionDistances_solvt, pbc, *x,
+                                           exclusionDistances);
     }
 
     if (max_sol > 0 && atoms_solvt->nres > max_sol)
diff --git a/src/gromacs/gmxpreprocess/tests/refdata/SolvateTest_cs_box_Works.xml b/src/gromacs/gmxpreprocess/tests/refdata/SolvateTest_cs_box_Works.xml
new file mode 100644 (file)
index 0000000..8bd8870
--- /dev/null
@@ -0,0 +1,12 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <OutputFiles Name="Files">
+    <File Name="-o">
+      <GroFile Name="Header">
+        <String Name="Title">Generated by gmx solvate</String>
+        <Int Name="Number of atoms">141</Int>
+      </GroFile>
+    </File>
+  </OutputFiles>
+</ReferenceData>
diff --git a/src/gromacs/gmxpreprocess/tests/refdata/SolvateTest_cs_cp_Works.xml b/src/gromacs/gmxpreprocess/tests/refdata/SolvateTest_cs_cp_Works.xml
new file mode 100644 (file)
index 0000000..b38a36e
--- /dev/null
@@ -0,0 +1,12 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <OutputFiles Name="Files">
+    <File Name="-o">
+      <GroFile Name="Header">
+        <String Name="Title">Test system for solvate/insert-molecules</String>
+        <Int Name="Number of atoms">2664</Int>
+      </GroFile>
+    </File>
+  </OutputFiles>
+</ReferenceData>
diff --git a/src/gromacs/gmxpreprocess/tests/refdata/SolvateTest_cs_cp_p_Works.xml b/src/gromacs/gmxpreprocess/tests/refdata/SolvateTest_cs_cp_p_Works.xml
new file mode 100644 (file)
index 0000000..b38a36e
--- /dev/null
@@ -0,0 +1,12 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <OutputFiles Name="Files">
+    <File Name="-o">
+      <GroFile Name="Header">
+        <String Name="Title">Test system for solvate/insert-molecules</String>
+        <Int Name="Number of atoms">2664</Int>
+      </GroFile>
+    </File>
+  </OutputFiles>
+</ReferenceData>
diff --git a/src/gromacs/gmxpreprocess/tests/refdata/SolvateTest_shell_Works.xml b/src/gromacs/gmxpreprocess/tests/refdata/SolvateTest_shell_Works.xml
new file mode 100644 (file)
index 0000000..f3cdafb
--- /dev/null
@@ -0,0 +1,12 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <OutputFiles Name="Files">
+    <File Name="-o">
+      <GroFile Name="Header">
+        <String Name="Title">Test system for solvate/insert-molecules</String>
+        <Int Name="Number of atoms">762</Int>
+      </GroFile>
+    </File>
+  </OutputFiles>
+</ReferenceData>
index 445e184f75f293e6f5c399bdeb84cdc167112b21..f1858fc9b538b1b2f596c46eafcb03f1f8da0746 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2017, 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/gmxpreprocess/solvate.h"
 
 #include "gromacs/utility/futil.h"
+#include "gromacs/utility/textreader.h"
 
 #include "testutils/cmdlinetest.h"
+#include "testutils/conftest.h"
+#include "testutils/refdata.h"
 #include "testutils/testfilemanager.h"
 #include "testutils/textblockmatchers.h"
 
@@ -54,14 +57,14 @@ namespace
 {
 
 using gmx::test::CommandLine;
-using gmx::test::NoTextMatch;
+using gmx::test::ConfMatch;
 
 class SolvateTest : public gmx::test::CommandLineTestBase
 {
     public:
         SolvateTest()
         {
-            setOutputFile("-o", "out.gro", NoTextMatch());
+            setOutputFile("-o", "out.gro", ConfMatch());
         }
 
         void runTest(const CommandLine &args)
@@ -70,6 +73,7 @@ class SolvateTest : public gmx::test::CommandLineTestBase
             cmdline.merge(args);
 
             ASSERT_EQ(0, gmx_solvate(cmdline.argc(), cmdline.argv()));
+            checkOutputFiles();
         }
 };
 
@@ -109,4 +113,16 @@ TEST_F(SolvateTest, cs_cp_p_Works)
     runTest(CommandLine(cmdline));
 }
 
+TEST_F(SolvateTest, shell_Works)
+{
+    // use default solvent box (-cs without argument)
+    const char *const cmdline[] = {
+        "solvate", "-cs"
+    };
+    setInputFile("-cp", "spc-and-methanol.gro");
+    commandLine().addOption("-shell", 1.0);
+
+    runTest(CommandLine(cmdline));
+}
+
 } // namespace
index 2bd34f4b64dc1b13c057bc8cbc2ed17eabf7fbc6..da1fc6a6603530da8b06e35fd1c8540abc21220d 100644 (file)
@@ -1,7 +1,7 @@
 #
 # This file is part of the GROMACS molecular simulation package.
 #
-# Copyright (c) 2011,2012,2013,2014,2015,2016, by the GROMACS development team, led by
+# Copyright (c) 2011,2012,2013,2014,2015,2016,2017, 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.
@@ -35,6 +35,7 @@
 include_directories(BEFORE SYSTEM ${GMOCK_INCLUDE_DIRS})
 set(TESTUTILS_SOURCES
     cmdlinetest.cpp
+    conftest.cpp
     integrationtests.cpp
     interactivetest.cpp
     mpi-printer.cpp
diff --git a/src/testutils/conftest.cpp b/src/testutils/conftest.cpp
new file mode 100644 (file)
index 0000000..dc591ee
--- /dev/null
@@ -0,0 +1,106 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2017, 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
+ * Implements routine to check the content of conf files.
+ *
+ * \author David van der Spoel <david.vanderspoel@icm.uu.se>
+ * \ingroup module_testutils
+ */
+#include "gmxpre.h"
+
+#include "conftest.h"
+
+#include <cstdio>
+#include <cstdlib>
+
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/stringutil.h"
+#include "gromacs/utility/textstream.h"
+
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+#include "testutils/textblockmatchers.h"
+
+namespace gmx
+{
+
+namespace test
+{
+
+namespace
+{
+
+class ConfMatcher : public ITextBlockMatcher
+{
+    public:
+        explicit ConfMatcher(const ConfMatchSettings &settings) : settings_(settings)
+        {
+        }
+
+        virtual void checkStream(TextInputStream      *stream,
+                                 TestReferenceChecker *checker)
+        {
+            checkConfFile(stream, checker, settings_);
+        }
+    private:
+        ConfMatchSettings settings_;
+};
+
+}       // namespace
+
+void checkConfFile(TextInputStream         *input,
+                   TestReferenceChecker    *checker,
+                   const ConfMatchSettings  &)
+{
+
+    TestReferenceChecker groChecker(checker->checkCompound("GroFile", "Header"));
+    // Just check the first two lines of the output file
+    std::string          line;
+    EXPECT_TRUE(input->readLine(&line));
+    line = stripSuffixIfPresent(line, "\n");
+    groChecker.checkString(line, "Title");
+    EXPECT_TRUE(input->readLine(&line));
+    line = stripSuffixIfPresent(line, "\n");
+    groChecker.checkInteger(std::atoi(line.c_str()), "Number of atoms");
+}
+
+TextBlockMatcherPointer ConfMatch::createMatcher() const
+{
+    return TextBlockMatcherPointer(new ConfMatcher(settings_));
+}
+
+} // namespace test
+} // namespace gmx
diff --git a/src/testutils/conftest.h b/src/testutils/conftest.h
new file mode 100644 (file)
index 0000000..1ef81de
--- /dev/null
@@ -0,0 +1,114 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2017, 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.
+ */
+/*! \libinternal \file
+ * \brief
+ * Declares function to add the content of a conf file to a checker.
+ *
+ * \author David van der Spoel <david.vanderspoel@icm.uu.se>
+ * \inlibraryapi
+ * \ingroup module_testutils
+ */
+#ifndef GMX_TESTUTILS_CONFTEST_H
+#define GMX_TESTUTILS_CONFTEST_H
+
+#include <string>
+
+#include "testutils/testasserts.h"
+#include "testutils/textblockmatchers.h"
+
+namespace gmx
+{
+
+class TextInputStream;
+
+namespace test
+{
+
+class TestReferenceChecker;
+
+struct ConfMatchSettings
+{
+    ConfMatchSettings() : tolerance(defaultRealTolerance())
+    {
+    }
+
+    FloatingPointTolerance  tolerance;
+};
+
+/*! \brief
+ * Adds content of a gro file to TestReferenceChecker object.
+ *
+ * \param[in] input       Stream that provides the gro content.
+ * \param[in,out] checker Checker to use.
+ * \param[in] settings    Settings to use for matching.
+ *
+ * Parses a gro file from the input stream, and checks the contents against
+ * reference data (only first two lines for now).
+ *
+ * \see ConfMatch
+ */
+void checkConfFile(TextInputStream         *input,
+                   TestReferenceChecker    *checker,
+                   const ConfMatchSettings &settings);
+
+/*! \libinternal \brief
+ * Match the contents as an gro file.
+ *
+ * \see checkGroFile()
+ *
+ * \inlibraryapi
+ * \ingroup module_testutils
+ */
+class ConfMatch : public ITextBlockMatcherSettings
+{
+    public:
+        //! Sets the tolerance for matching floating point values.
+        ConfMatch &tolerance(const FloatingPointTolerance &tolerance)
+        {
+            settings_.tolerance = tolerance;
+            return *this;
+        }
+
+        virtual TextBlockMatcherPointer createMatcher() const;
+
+    private:
+        ConfMatchSettings  settings_;
+};
+
+} // namespace test
+
+} // namespace gmx
+
+#endif