Add nblib tpr reading
authorJoe Jordan <ejjordan12@gmail.com>
Fri, 8 Oct 2021 14:22:08 +0000 (14:22 +0000)
committerJoe Jordan <ejjordan12@gmail.com>
Fri, 8 Oct 2021 14:22:08 +0000 (14:22 +0000)
api/nblib/CMakeLists.txt
api/nblib/gmxcalculatorcpu.cpp
api/nblib/gmxcalculatorcpu.h
api/nblib/tests/CMakeLists.txt
api/nblib/tests/tpr.cpp [new file with mode: 0644]
api/nblib/tpr.cpp [new file with mode: 0644]
api/nblib/tpr.h [new file with mode: 0644]

index 3619ddb510ee0e76199e1c22114e478be637308a..fcf9b2868fb7daa2e37c743d403a15cdf00ce6e5 100644 (file)
@@ -99,6 +99,7 @@ target_sources(nblib
         simulationstate.cpp
         topologyhelpers.cpp
         topology.cpp
+        tpr.cpp
         virials.cpp
         )
 
@@ -141,6 +142,7 @@ if(GMX_INSTALL_NBLIB_API)
             particletype.h
             simulationstate.h
             topology.h
+            tpr.h
             vector.h
             DESTINATION include/nblib)
 endif()
index a28a3760657bd29e0f63551a5757a0f60d15fa2d..a976e5d370c9faee432d634376642f5660c05e44 100644 (file)
@@ -57,6 +57,7 @@
 #include "nblib/pbc.hpp"
 #include "nblib/systemdescription.h"
 #include "nblib/topology.h"
+#include "nblib/tpr.h"
 #include "nblib/virials.h"
 
 namespace nblib
@@ -336,4 +337,16 @@ std::unique_ptr<GmxNBForceCalculatorCpu> setupGmxForceCalculatorCpu(const Topolo
                                                      options);
 }
 
+std::unique_ptr<GmxNBForceCalculatorCpu> setupGmxForceCalculatorCpu(TprReader& tprReader,
+                                                                    const NBKernelOptions& options)
+{
+    return std::make_unique<GmxNBForceCalculatorCpu>(tprReader.particleTypeIdOfAllParticles_,
+                                                     tprReader.nonbondedParameters_,
+                                                     tprReader.charges_,
+                                                     tprReader.particleInteractionFlags_,
+                                                     tprReader.exclusionListRanges_,
+                                                     tprReader.exclusionListElements_,
+                                                     options);
+}
+
 } // namespace nblib
index 2fec548256a8b7fce9d7050efc0e9490d66bf84d..efc079896ef7d7baccbcfa5f9b39328b24a146c0 100644 (file)
@@ -63,6 +63,7 @@ namespace nblib
 {
 struct NBKernelOptions;
 class Topology;
+struct TprReader;
 
 class GmxNBForceCalculatorCpu final
 {
@@ -107,6 +108,10 @@ private:
 std::unique_ptr<GmxNBForceCalculatorCpu> setupGmxForceCalculatorCpu(const Topology&        topology,
                                                                     const NBKernelOptions& options);
 
+//! Sets up and returns a GmxForceCalculatorCpu based on a TPR file as input
+std::unique_ptr<GmxNBForceCalculatorCpu> setupGmxForceCalculatorCpu(TprReader& tprReader,
+                                                                    const NBKernelOptions& options);
+
 } // namespace nblib
 
 #endif // NBLIB_GMXCALCULATORCPU_H
index 6346ff29fb326058ea824e90696d6e4574b8e111..ef2ace7d19a30b4bd0a653694a8d6648cef10fd3 100644 (file)
@@ -61,6 +61,7 @@ gmx_add_gtest_executable(
         molecules.cpp
         nbnxmsetup.cpp
         topology.cpp
+        tpr.cpp
         virials.cpp
     )
 target_link_libraries(${exename} PRIVATE mdrun_test_infrastructure)
diff --git a/api/nblib/tests/tpr.cpp b/api/nblib/tests/tpr.cpp
new file mode 100644 (file)
index 0000000..0c39288
--- /dev/null
@@ -0,0 +1,146 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 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.
+ *
+ * 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
+ * This implements TPR reading tests
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ */
+#include <gtest/gtest.h>
+
+#include "nblib/tpr.h"
+
+#include "gromacs/tools/convert_tpr.h"
+
+#include "nblib/gmxcalculatorcpu.h"
+#include "nblib/tests/testsystems.h"
+#include "nblib/tests/testhelpers.h"
+
+#include "programs/mdrun/tests/moduletest.h"
+
+#include "testutils/cmdlinetest.h"
+#include "testutils/simulationdatabase.h"
+
+namespace nblib
+{
+namespace test
+{
+
+//! Test fixture for running the utils needed for reading TPR files from simulation database
+class TprReaderTest : public gmx::test::MdrunTestFixture
+{
+public:
+    //! Returns a TprReader object populated with the given files in the simulation database
+    TprReader makeTPRfromSimulationDatabase(const std::string& filename)
+    {
+        gmx::test::CommandLine caller;
+        runner_.useTopGroAndNdxFromDatabase(filename);
+        auto mdpFieldValues = gmx::test::prepareMdpFieldValues(filename.c_str(), "md", "no", "no");
+        mdpFieldValues["init-lambda-state"] = "0";
+        mdpFieldValues["constraints"]       = "none"; // constraints not supported in NBLIB
+        runner_.useStringAsMdpFile(gmx::test::prepareMdpFileContents(mdpFieldValues));
+        runner_.callGrompp(caller);
+
+        return TprReader(runner_.tprFileName_);
+    }
+};
+
+TEST_F(TprReaderTest, SimDBTprIsCreated)
+{
+    EXPECT_NO_THROW(makeTPRfromSimulationDatabase("argon12"));
+}
+
+TEST_F(TprReaderTest, Spc2Reads)
+{
+    TprReader tprReader = makeTPRfromSimulationDatabase("spc2");
+
+    EXPECT_NO_THROW(tprReader.coordinates_.size());
+}
+
+TEST_F(TprReaderTest, ArgonImportedDataIsCorrect)
+{
+    TprReader      tprReader = makeTPRfromSimulationDatabase("argon12");
+    RefDataChecker dataImportTest(5e-5);
+
+    // check number of coordinates in the system
+    EXPECT_EQ(tprReader.coordinates_.size(), 12);
+
+    // check coordinates
+    dataImportTest.testArrays<Vec3>(tprReader.coordinates_, "coordinates");
+
+    // check number of velocities in the system
+    EXPECT_EQ(tprReader.velocities_.size(), 12);
+
+    // check velocities
+    dataImportTest.testArrays<Vec3>(tprReader.velocities_, "velocities");
+
+    // check box
+    EXPECT_TRUE(tprReader.getBox() == Box(6.05449));
+
+    // check non-bonded params
+    dataImportTest.testArrays<real>(tprReader.nonbondedParameters_, "nbparams");
+
+    // check exclusion elements
+    dataImportTest.testArrays<int>(tprReader.exclusionListElements_, "exclusion elements");
+
+    // check exclusion ranges
+    dataImportTest.testArrays<int>(tprReader.exclusionListRanges_, "exclusion ranges");
+}
+
+TEST_F(TprReaderTest, FCfromTprDataWorks)
+{
+    auto options        = NBKernelOptions();
+    options.nbnxmSimd   = SimdKernels::SimdNo;
+    options.coulombType = CoulombType::Cutoff;
+
+    TprReader tprReader = makeTPRfromSimulationDatabase("argon12");
+
+    auto forceCalculatorTpr = setupGmxForceCalculatorCpu(tprReader, options);
+
+    std::vector<Vec3> forceOutputTpr(tprReader.coordinates_.size(), Vec3{ 0, 0, 0 });
+
+    forceCalculatorTpr->updatePairlist(tprReader.coordinates_, tprReader.getBox());
+    forceCalculatorTpr->compute(tprReader.coordinates_, tprReader.getBox(), forceOutputTpr);
+
+    RefDataChecker forcesOutputTest;
+    forcesOutputTest.testArrays<Vec3>(forceOutputTpr, "Argon forces");
+}
+
+} // namespace test
+
+} // namespace nblib
diff --git a/api/nblib/tpr.cpp b/api/nblib/tpr.cpp
new file mode 100644 (file)
index 0000000..f534dcd
--- /dev/null
@@ -0,0 +1,145 @@
+/*
+ * 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.
+ */
+/*! \inpublicapi \file
+ * \brief
+ * Implements nblib tpr reading
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ */
+
+#include "nblib/tpr.h"
+
+#include "gromacs/fileio/tpxio.h"
+#include "gromacs/gmxlib/network.h"
+#include "gromacs/listed_forces/listed_forces.h"
+#include "gromacs/mdtypes/forcerec.h"
+#include "gromacs/mdlib/mdatoms.h"
+#include "gromacs/mdlib/forcerec.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/iforceprovider.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/mdatom.h"
+#include "gromacs/mdtypes/simulation_workload.h"
+#include "gromacs/mdtypes/state.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/logger.h"
+#include "nblib/box.h"
+
+namespace nblib
+{
+
+TprReader::TprReader(std::string filename)
+{
+    t_inputrec              inputRecord;
+    t_state                 globalState;
+    gmx_mtop_t              molecularTopology;
+    gmx::SimulationWorkload simulationWorkload;
+
+    // If the file does not exist, this function will throw
+    PartialDeserializedTprFile partialDeserializedTpr =
+            read_tpx_state(filename.c_str(), &inputRecord, &globalState, &molecularTopology);
+
+    // init commrec
+    MPI_Comm      simulationCommunicator = MPI_COMM_WORLD;
+    CommrecHandle crHandle               = init_commrec(simulationCommunicator);
+    t_commrec*    commrec                = crHandle.get();
+    assert((commrec != nullptr) && "Must have valid commrec");
+
+    // init forcerec
+    t_forcerec          forceRecord;
+    gmx::ForceProviders forceProviders;
+    forceRecord.forceProviders = &forceProviders;
+    init_forcerec(nullptr,
+                  gmx::MDLogger(),
+                  simulationWorkload,
+                  &forceRecord,
+                  inputRecord,
+                  molecularTopology,
+                  commrec,
+                  globalState.box,
+                  nullptr,
+                  nullptr,
+                  gmx::ArrayRef<const std::string>{},
+                  -1);
+
+    nonbondedParameters_ = makeNonBondedParameterLists(molecularTopology.ffparams.atnr,
+                                                       molecularTopology.ffparams.iparams,
+                                                       forceRecord.haveBuckingham);
+
+    gmx_localtop_t localtop(molecularTopology.ffparams);
+    gmx_mtop_generate_local_top(molecularTopology, &localtop, false);
+    exclusionListElements_ = std::vector<int>(localtop.excls.elementsView().begin(),
+                                              localtop.excls.elementsView().end());
+    exclusionListRanges_   = std::vector<int>(localtop.excls.listRangesView().begin(),
+                                            localtop.excls.listRangesView().end());
+
+    int                           ntopatoms = molecularTopology.natoms;
+    std::unique_ptr<gmx::MDAtoms> mdAtoms =
+            gmx::makeMDAtoms(nullptr, molecularTopology, inputRecord, false);
+    atoms2md(molecularTopology, inputRecord, -1, {}, ntopatoms, mdAtoms.get());
+    update_mdatoms(mdAtoms->mdatoms(), inputRecord.fepvals->init_lambda);
+    auto numParticles = mdAtoms->mdatoms()->nr;
+    charges_.resize(numParticles);
+    particleTypeIdOfAllParticles_.resize(numParticles);
+    inverseMasses_.resize(numParticles);
+    for (int i = 0; i < numParticles; i++)
+    {
+        charges_[i]                      = mdAtoms->mdatoms()->chargeA[i];
+        particleTypeIdOfAllParticles_[i] = mdAtoms->mdatoms()->typeA[i];
+        inverseMasses_[i]                = mdAtoms->mdatoms()->invmass[i];
+    }
+    particleInteractionFlags_ = forceRecord.atomInfo;
+
+    if (TRICLINIC(globalState.box))
+    {
+        throw InputException("Only rectangular unit-cells are supported here");
+    }
+    boxX_ = globalState.box[XX][XX];
+    boxY_ = globalState.box[YY][YY];
+    boxZ_ = globalState.box[ZZ][ZZ];
+    coordinates_.assign(globalState.x.begin(), globalState.x.end());
+    velocities_.assign(globalState.v.begin(), globalState.v.end());
+}
+
+Box TprReader::getBox() const
+{
+    return Box(boxX_, boxY_, boxZ_);
+}
+
+} // namespace nblib
diff --git a/api/nblib/tpr.h b/api/nblib/tpr.h
new file mode 100644 (file)
index 0000000..2ed3f3c
--- /dev/null
@@ -0,0 +1,107 @@
+/*
+ * 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.
+ */
+/*! \inpublicapi \file
+ * \brief
+ * Implements nblib tpr reading
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ */
+#ifndef NBLIB_TPR_H
+#define NBLIB_TPR_H
+
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "nblib/basicdefinitions.h"
+#include "nblib/box.h"
+#include "nblib/vector.h"
+#include "nblib/topology.h"
+
+namespace nblib
+{
+template<typename T>
+struct ExclusionLists;
+class Box;
+
+/*! \brief Reads in data from tpr which can be used in construction of force calculators
+ *
+ *  This object reads data such as topology and coordinates from a tpr file, with the
+ *  intent of using that data for construction of nblib force calculator objects. To
+ *  allow for maximum flexibility in the calling code, it is the responsibility of the
+ *  caller to, e.g., make copies if they want to preserve initial state information or
+ *  use the data to construct multiple force calculator objects.
+ */
+struct TprReader
+{
+public:
+    TprReader(std::string filename);
+
+    //! Particle info where all particles are marked to have Van der Waals interactions
+    std::vector<int64_t> particleInteractionFlags_;
+    //! particle type id of all particles
+    std::vector<int> particleTypeIdOfAllParticles_;
+    //! Storage for parameters for short range interactions.
+    std::vector<real> nonbondedParameters_;
+    //! electrostatic charges
+    std::vector<real> charges_;
+    //! inverse masses
+    std::vector<real> inverseMasses_;
+    //! stores information about particles pairs to be excluded from the non-bonded force
+    //! calculation exclusion list ranges
+    std::vector<int> exclusionListRanges_;
+    //! exclusion list elements
+    std::vector<int> exclusionListElements_;
+    //! coordinates
+    std::vector<Vec3> coordinates_;
+    //! velocities
+    std::vector<Vec3> velocities_;
+    //! listed forces data
+    ListedInteractionData listedInteractionData_;
+
+    //! bounding box of particle coordinates
+    [[nodiscard]] Box getBox() const;
+
+private:
+    real boxX_;
+    real boxY_;
+    real boxZ_;
+};
+
+} // namespace nblib
+#endif // NBLIB_TPR_H