Introduce nblib nonbonded force calculator impl class
authorJoe Jordan <ejjordan12@gmail.com>
Thu, 7 Oct 2021 07:35:17 +0000 (07:35 +0000)
committerMark Abraham <mark.j.abraham@gmail.com>
Thu, 7 Oct 2021 07:35:17 +0000 (07:35 +0000)
15 files changed:
api/nblib/CMakeLists.txt
api/nblib/gmxbackenddata.h [moved from api/nblib/gmxcalculator.h with 57% similarity]
api/nblib/gmxcalculator.cpp [deleted file]
api/nblib/gmxcalculatorcpu.cpp [new file with mode: 0644]
api/nblib/gmxcalculatorcpu.h [new file with mode: 0644]
api/nblib/nblib.h
api/nblib/nbnxmsetuphelpers.cpp
api/nblib/nbnxmsetuphelpers.h
api/nblib/samples/argon-forces-integration.cpp
api/nblib/samples/methane-water-integration.cpp
api/nblib/systemdescription.h [new file with mode: 0644]
api/nblib/tests/gmxcalculator.cpp
api/nblib/tests/nbkernelsystem.cpp
api/nblib/tests/refdata/NBlibTest_ArgonEnergiesAreCorrect.xml [new file with mode: 0644]
api/nblib/tests/refdata/NBlibTest_ArgonVirialsAreCorrect.xml [new file with mode: 0644]

index 31eca82e77a6ffb630b1c8f03cc5927e18f673fd..3619ddb510ee0e76199e1c22114e478be637308a 100644 (file)
@@ -52,14 +52,14 @@ add_custom_target(nblib-tests
 # Ensure that "make tests" builds all nblib tests so the top-level
 # "make check" will work.
 if (BUILD_TESTING)
-       add_dependencies(tests nblib-tests)
+        add_dependencies(tests nblib-tests)
 
-       # this allows all nblib tests to be run with "make check-nblib"
-       add_custom_target(check-nblib
-               COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -R NbLib
-               COMMENT "Running nblib tests"
-               USES_TERMINAL VERBATIM)
-       add_dependencies(check-nblib nblib-tests)
+        # this allows all nblib tests to be run with "make check-nblib"
+        add_custom_target(check-nblib
+                COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -R NbLib
+                COMMENT "Running nblib tests"
+                USES_TERMINAL VERBATIM)
+        add_dependencies(check-nblib nblib-tests)
 endif()
 
 set(NBLIB_MAJOR 0)
@@ -89,8 +89,7 @@ endif()
 target_sources(nblib
         PRIVATE
         box.cpp
-        gmxcalculator.cpp
-        gmxsetup.cpp
+        gmxcalculatorcpu.cpp
         integrator.cpp
         interactions.cpp
         molecules.cpp
@@ -131,7 +130,7 @@ if(GMX_INSTALL_NBLIB_API)
             basicdefinitions.h
             box.h
             exception.h
-            gmxcalculator.h
+            gmxcalculatorcpu.h
             integrator.h
             interactions.h
             molecules.h
similarity index 57%
rename from api/nblib/gmxcalculator.h
rename to api/nblib/gmxbackenddata.h
index 71a56ebb08adbe3d9adb66a28e03e555018a4237..6794bb47debd747241e8fc99336e816c9e238330 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2020,2021, by the GROMACS development team, led by
+ * 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.
  * 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
+/*! \inpublicapi \file
  * \brief
- * Implements a force calculator based on GROMACS data structures.
- *
- * Intended for internal use inside the ForceCalculator.
+ * Implements nblib simulation box
  *
  * \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_GMXCALCULATOR_H
-#define NBLIB_GMXCALCULATOR_H
-
-#include <memory>
-
-#include "nblib/vector.h"
-
-struct nonbonded_verlet_t;
-struct t_forcerec;
-struct t_nrnb;
-struct interaction_const_t;
-struct gmx_enerdata_t;
-
-namespace gmx
-{
-template<typename T>
-class ArrayRef;
-class StepWorkload;
-} // namespace gmx
+#ifndef NBLIB_GMXBACKENDDATA_H
+#define NBLIB_GMXBACKENDDATA_H
+
+#include "gromacs/gmxlib/nrnb.h"
+#include "gromacs/mdtypes/enerdata.h"
+#include "gromacs/mdtypes/forcerec.h"
+#include "gromacs/mdtypes/interaction_const.h"
+#include "gromacs/mdtypes/simulation_workload.h"
+#include "gromacs/nbnxm/nbnxm.h"
+#include "gromacs/timing/wallcycle.h"
+#include "gromacs/utility/listoflists.h"
+#include "gromacs/utility/range.h"
+#include "nblib/kerneloptions.h"
+#include "nblib/nbnxmsetuphelpers.h"
 
 namespace nblib
 {
-class Box;
-class NbvSetupUtil;
-class SimulationState;
-struct NBKernelOptions;
 
 /*! \brief GROMACS non-bonded force calculation backend
  *
@@ -82,48 +70,59 @@ struct NBKernelOptions;
  * interactions (PP) on a single node CPU.
  *
  */
-
-class GmxForceCalculator final
+class GmxBackendData
 {
 public:
-    GmxForceCalculator();
+    GmxBackendData() = default;
+    GmxBackendData(const NBKernelOptions& options,
+                   int                    numEnergyGroups,
+                   gmx::ArrayRef<int>     exclusionRanges,
+                   gmx::ArrayRef<int>     exclusionElements)
+    {
+        // Set hardware params from the execution context
+        setGmxNonBondedNThreads(options.numOpenMPThreads);
+
+        // Set interaction constants struct
+        interactionConst_ = createInteractionConst(options);
 
-    ~GmxForceCalculator();
+        // Set up StepWorkload data
+        stepWork_ = createStepWorkload();
 
-    //! Compute forces and return
-    void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput, gmx::ArrayRef<gmx::RVec> forceOutput);
+        // Set up gmx_enerdata_t (holds energy information)
+        enerd_ = gmx_enerdata_t{ numEnergyGroups, 0 };
 
-    //! Puts particles on a grid based on bounds specified by the box (for every NS step)
-    void setParticlesOnGrid(gmx::ArrayRef<const int64_t>   particleInfoAllVdw,
-                            gmx::ArrayRef<const gmx::RVec> coordinates,
-                            const Box&                     box);
+        // Construct pair lists
+        std::vector<int> exclusionRanges_(exclusionRanges.begin(), exclusionRanges.end());
+        std::vector<int> exclusionElements_(exclusionElements.begin(), exclusionElements.end());
+        exclusions_ = gmx::ListOfLists<int>(std::move(exclusionRanges_), std::move(exclusionElements_));
+    }
 
-private:
-    //! Friend to allow setting up private members in this class
-    friend class NbvSetupUtil;
+    //! exclusions in gmx format
+    gmx::ListOfLists<int> exclusions_;
 
     //! Non-Bonded Verlet object for force calculation
     std::unique_ptr<nonbonded_verlet_t> nbv_;
 
-    //! Only nbfp and shift_vec are used
-    std::unique_ptr<t_forcerec> forcerec_;
+    //! Only shift_vec is used
+    t_forcerec forcerec_;
 
     //! Parameters for various interactions in the system
-    std::unique_ptr<interaction_const_t> interactionConst_;
+    interaction_const_t interactionConst_;
 
     //! Tasks to perform in an MD Step
-    std::unique_ptr<gmx::StepWorkload> stepWork_;
+    gmx::StepWorkload stepWork_;
+
+    gmx::SimulationWorkload simulationWork_;
 
     //! Energies of different interaction types; currently only needed as an argument for dispatchNonbondedKernel
-    std::unique_ptr<gmx_enerdata_t> enerd_;
+    gmx_enerdata_t enerd_{ 1, 0 };
 
     //! Non-bonded flop counter; currently only needed as an argument for dispatchNonbondedKernel
-    std::unique_ptr<t_nrnb> nrnb_;
+    t_nrnb nrnb_;
 
-    //! Legacy matrix for box
-    matrix box_{ { 0 } };
+    //! Keep track of whether updatePairlist has been called at least once
+    bool updatePairlistCalled{ false };
 };
 
 } // namespace nblib
-
-#endif // NBLIB_GMXCALCULATOR_H
+#endif // NBLIB_GMXBACKENDDATA_H
diff --git a/api/nblib/gmxcalculator.cpp b/api/nblib/gmxcalculator.cpp
deleted file mode 100644 (file)
index 81e740f..0000000
+++ /dev/null
@@ -1,119 +0,0 @@
-/*
- * 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 Implements a force calculator based on GROMACS data structures.
- *
- * \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/gmxcalculator.h"
-#include "gromacs/ewald/ewald_utils.h"
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/mdlib/rf_util.h"
-#include "gromacs/mdtypes/enerdata.h"
-#include "gromacs/mdtypes/forcerec.h"
-#include "gromacs/mdtypes/interaction_const.h"
-#include "gromacs/mdtypes/simulation_workload.h"
-#include "gromacs/nbnxm/nbnxm.h"
-#include "gromacs/utility/range.h"
-#include "nblib/exception.h"
-#include "nblib/simulationstate.h"
-
-namespace nblib
-{
-
-GmxForceCalculator::GmxForceCalculator()
-{
-    enerd_            = std::make_unique<gmx_enerdata_t>(1, 0);
-    forcerec_         = std::make_unique<t_forcerec>();
-    interactionConst_ = std::make_unique<interaction_const_t>();
-    stepWork_         = std::make_unique<gmx::StepWorkload>();
-    nrnb_             = std::make_unique<t_nrnb>();
-}
-
-GmxForceCalculator::~GmxForceCalculator() = default;
-
-void GmxForceCalculator::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
-                                 gmx::ArrayRef<gmx::RVec>       forceOutput)
-{
-    // update the coordinates in the backend
-    nbv_->convertCoordinates(gmx::AtomLocality::Local, coordinateInput);
-
-    nbv_->dispatchNonbondedKernel(
-            gmx::InteractionLocality::Local,
-            *interactionConst_,
-            *stepWork_,
-            enbvClearFYes,
-            forcerec_->shift_vec,
-            enerd_->grpp.energyGroupPairTerms[forcerec_->haveBuckingham ? NonBondedEnergyTerms::BuckinghamSR
-                                                                        : NonBondedEnergyTerms::LJSR],
-            enerd_->grpp.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR],
-            nrnb_.get());
-
-    nbv_->atomdata_add_nbat_f_to_f(gmx::AtomLocality::All, forceOutput);
-}
-
-void GmxForceCalculator::setParticlesOnGrid(gmx::ArrayRef<const int64_t>   particleInfoAllVdw,
-                                            gmx::ArrayRef<const gmx::RVec> coordinates,
-                                            const Box&                     box)
-{
-    auto legacyBox = box.legacyMatrix();
-
-    if (TRICLINIC(legacyBox))
-    {
-        throw InputException("Only rectangular unit-cells are supported here");
-    }
-    const rvec lowerCorner = { 0, 0, 0 };
-    const rvec upperCorner = { legacyBox[dimX][dimX], legacyBox[dimY][dimY], legacyBox[dimZ][dimZ] };
-
-    const real particleDensity = coordinates.size() / det(legacyBox);
-
-    nbnxn_put_on_grid(nbv_.get(),
-                      legacyBox,
-                      0,
-                      lowerCorner,
-                      upperCorner,
-                      nullptr,
-                      { 0, int(coordinates.size()) },
-                      particleDensity,
-                      particleInfoAllVdw,
-                      coordinates,
-                      0,
-                      nullptr);
-}
-
-} // namespace nblib
diff --git a/api/nblib/gmxcalculatorcpu.cpp b/api/nblib/gmxcalculatorcpu.cpp
new file mode 100644 (file)
index 0000000..a28a376
--- /dev/null
@@ -0,0 +1,339 @@
+/*
+ * 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 Implements a force calculator based on GROMACS data structures.
+ *
+ * \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 "gromacs/ewald/ewald_utils.h"
+#include "gromacs/mdtypes/enerdata.h"
+#include "gromacs/mdtypes/interaction_const.h"
+#include "gromacs/nbnxm/atomdata.h"
+#include "gromacs/nbnxm/nbnxm.h"
+#include "gromacs/nbnxm/pairlistset.h"
+#include "gromacs/nbnxm/pairlistsets.h"
+#include "gromacs/nbnxm/pairsearch.h"
+#include "gromacs/utility/listoflists.h"
+#include "gromacs/utility/range.h"
+#include "nblib/exception.h"
+#include "nblib/gmxbackenddata.h"
+#include "nblib/gmxcalculatorcpu.h"
+#include "nblib/nbnxmsetuphelpers.h"
+#include "nblib/pbc.hpp"
+#include "nblib/systemdescription.h"
+#include "nblib/topology.h"
+#include "nblib/virials.h"
+
+namespace nblib
+{
+
+class GmxNBForceCalculatorCpu::CpuImpl final
+{
+public:
+    CpuImpl(gmx::ArrayRef<int>     particleTypeIdOfAllParticles,
+            gmx::ArrayRef<real>    nonBondedParams,
+            gmx::ArrayRef<real>    charges,
+            gmx::ArrayRef<int64_t> particleInteractionFlags,
+            gmx::ArrayRef<int>     exclusionRanges,
+            gmx::ArrayRef<int>     exclusionElements,
+            const NBKernelOptions& options);
+
+    //! calculates a new pair list based on new coordinates (for every NS step)
+    void updatePairlist(gmx::ArrayRef<const gmx::RVec> coordinates, const Box& box);
+
+    //! Compute forces and return
+    void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
+                 const Box&                     box,
+                 gmx::ArrayRef<gmx::RVec>       forceOutput);
+
+    //! Compute forces and virial tensor
+    void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
+                 const Box&                     box,
+                 gmx::ArrayRef<gmx::RVec>       forceOutput,
+                 gmx::ArrayRef<real>            virialOutput);
+
+    //! Compute forces, virial tensor and potential energies
+    void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
+                 const Box&                     box,
+                 gmx::ArrayRef<gmx::RVec>       forceOutput,
+                 gmx::ArrayRef<real>            virialOutput,
+                 gmx::ArrayRef<real>            energyOutput);
+
+private:
+    //! \brief client-side provided system description data
+    SystemDescription system_;
+
+    //! \brief Gmx backend objects, employed for calculating the forces
+    GmxBackendData backend_;
+};
+
+GmxNBForceCalculatorCpu::CpuImpl::CpuImpl(gmx::ArrayRef<int>     particleTypeIdOfAllParticles,
+                                          gmx::ArrayRef<real>    nonBondedParams,
+                                          gmx::ArrayRef<real>    charges,
+                                          gmx::ArrayRef<int64_t> particleInteractionFlags,
+                                          gmx::ArrayRef<int>     exclusionRanges,
+                                          gmx::ArrayRef<int>     exclusionElements,
+                                          const NBKernelOptions& options) :
+    system_(SystemDescription(particleTypeIdOfAllParticles, nonBondedParams, charges, particleInteractionFlags)),
+    backend_(GmxBackendData(options, findNumEnergyGroups(particleInteractionFlags), exclusionRanges, exclusionElements))
+{
+    // Set up non-bonded verlet in the backend
+    backend_.nbv_ = createNbnxmCPU(system_.numParticleTypes_,
+                                   options,
+                                   findNumEnergyGroups(particleInteractionFlags),
+                                   system_.nonBondedParams_);
+}
+
+void GmxNBForceCalculatorCpu::CpuImpl::updatePairlist(gmx::ArrayRef<const gmx::RVec> coordinates,
+                                                      const Box&                     box)
+{
+    if (coordinates.size() != system_.numParticles_)
+    {
+        throw InputException(
+                "Coordinate array containing different number of entries than particles in the "
+                "system");
+    }
+
+    const auto* legacyBox = box.legacyMatrix();
+    system_.box_          = box;
+    updateForcerec(&backend_.forcerec_, box.legacyMatrix());
+    if (TRICLINIC(legacyBox))
+    {
+        throw InputException("Only rectangular unit-cells are supported here");
+    }
+
+    const rvec lowerCorner = { 0, 0, 0 };
+    const rvec upperCorner = { legacyBox[dimX][dimX], legacyBox[dimY][dimY], legacyBox[dimZ][dimZ] };
+
+    const real particleDensity = static_cast<real>(coordinates.size()) / det(legacyBox);
+
+    // Put particles on a grid based on bounds specified by the box
+    nbnxn_put_on_grid(backend_.nbv_.get(),
+                      legacyBox,
+                      0,
+                      lowerCorner,
+                      upperCorner,
+                      nullptr,
+                      { 0, int(coordinates.size()) },
+                      particleDensity,
+                      system_.particleInfo_,
+                      coordinates,
+                      0,
+                      nullptr);
+
+    backend_.nbv_->constructPairlist(
+            gmx::InteractionLocality::Local, backend_.exclusions_, 0, &backend_.nrnb_);
+
+    // Set Particle Types and Charges and VdW params
+    backend_.nbv_->setAtomProperties(
+            system_.particleTypeIdOfAllParticles_, system_.charges_, system_.particleInfo_);
+    backend_.updatePairlistCalled = true;
+}
+
+void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
+                                               const Box&                     box,
+                                               gmx::ArrayRef<gmx::RVec>       forceOutput,
+                                               gmx::ArrayRef<real>            virialOutput,
+                                               gmx::ArrayRef<real>            energyOutput)
+{
+    if (coordinateInput.size() != forceOutput.size())
+    {
+        throw InputException("coordinate array and force buffer size mismatch");
+    }
+
+    if (!backend_.updatePairlistCalled)
+    {
+        throw InputException("compute called without updating pairlist at least once");
+    }
+
+    // update the box if changed
+    if (!(system_.box_ == box))
+    {
+        system_.box_ = box;
+        updateForcerec(&backend_.forcerec_, box.legacyMatrix());
+    }
+
+    bool computeVirial               = !virialOutput.empty();
+    bool computeEnergies             = !energyOutput.empty();
+    backend_.stepWork_.computeVirial = computeVirial;
+    backend_.stepWork_.computeEnergy = computeEnergies;
+
+    // update the coordinates in the backend
+    backend_.nbv_->convertCoordinates(gmx::AtomLocality::Local, coordinateInput);
+
+    backend_.nbv_->dispatchNonbondedKernel(
+            gmx::InteractionLocality::Local,
+            backend_.interactionConst_,
+            backend_.stepWork_,
+            enbvClearFYes,
+            backend_.forcerec_.shift_vec,
+            backend_.enerd_.grpp.energyGroupPairTerms[backend_.forcerec_.haveBuckingham ? NonBondedEnergyTerms::BuckinghamSR
+                                                                                        : NonBondedEnergyTerms::LJSR],
+            backend_.enerd_.grpp.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR],
+            &backend_.nrnb_);
+
+    backend_.nbv_->atomdata_add_nbat_f_to_f(gmx::AtomLocality::All, forceOutput);
+
+    if (computeVirial)
+    {
+        // calculate shift forces and turn into an array ref
+        std::vector<Vec3> shiftForcesVector(gmx::c_numShiftVectors, Vec3(0.0, 0.0, 0.0));
+        nbnxn_atomdata_add_nbat_fshift_to_fshift(*backend_.nbv_->nbat, shiftForcesVector);
+        auto shiftForcesRef = constArrayRefFromArray(shiftForcesVector.data(), shiftForcesVector.size());
+
+        std::vector<Vec3> shiftVectorsArray(gmx::c_numShiftVectors);
+
+        // copy shift vectors from ForceRec
+        std::copy(backend_.forcerec_.shift_vec.begin(),
+                  backend_.forcerec_.shift_vec.end(),
+                  shiftVectorsArray.begin());
+
+        computeVirialTensor(
+                coordinateInput, forceOutput, shiftVectorsArray, shiftForcesRef, box, virialOutput);
+    }
+
+    // extract term energies (per interaction type)
+    if (computeEnergies)
+    {
+        int nGroupPairs = backend_.enerd_.grpp.nener;
+        if (int(energyOutput.size()) != int(NonBondedEnergyTerms::Count) * nGroupPairs)
+        {
+            throw InputException("Array size for energy output is wrong\n");
+        }
+
+        for (int eg = 0; eg < int(NonBondedEnergyTerms::Count); ++eg)
+        {
+            std::copy(begin(backend_.enerd_.grpp.energyGroupPairTerms[eg]),
+                      end(backend_.enerd_.grpp.energyGroupPairTerms[eg]),
+                      energyOutput.begin() + eg * nGroupPairs);
+        }
+    }
+}
+
+void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
+                                               const Box&                     box,
+                                               gmx::ArrayRef<gmx::RVec>       forceOutput)
+{
+    // compute forces and fill in force buffer
+    compute(coordinateInput, box, forceOutput, gmx::ArrayRef<real>{}, gmx::ArrayRef<real>{});
+}
+
+void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
+                                               const Box&                     box,
+                                               gmx::ArrayRef<gmx::RVec>       forceOutput,
+                                               gmx::ArrayRef<real>            virialOutput)
+{
+    // compute forces and fill in force buffer
+    compute(coordinateInput, box, forceOutput, virialOutput, gmx::ArrayRef<real>{});
+}
+
+GmxNBForceCalculatorCpu::GmxNBForceCalculatorCpu(gmx::ArrayRef<int>  particleTypeIdOfAllParticles,
+                                                 gmx::ArrayRef<real> nonBondedParams,
+                                                 gmx::ArrayRef<real> charges,
+                                                 gmx::ArrayRef<int64_t> particleInteractionFlags,
+                                                 gmx::ArrayRef<int>     exclusionRanges,
+                                                 gmx::ArrayRef<int>     exclusionElements,
+                                                 const NBKernelOptions& options)
+{
+    if (options.useGpu)
+    {
+        throw InputException("Use GmxNBForceCalculatorGpu for GPU support");
+    }
+
+    impl_ = std::make_unique<CpuImpl>(particleTypeIdOfAllParticles,
+                                      nonBondedParams,
+                                      charges,
+                                      particleInteractionFlags,
+                                      exclusionRanges,
+                                      exclusionElements,
+                                      options);
+}
+
+GmxNBForceCalculatorCpu::~GmxNBForceCalculatorCpu() = default;
+
+//! calculates a new pair list based on new coordinates (for every NS step)
+void GmxNBForceCalculatorCpu::updatePairlist(gmx::ArrayRef<const gmx::RVec> coordinates, const Box& box)
+{
+    impl_->updatePairlist(coordinates, box);
+}
+
+//! Compute forces and return
+void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
+                                      const Box&                     box,
+                                      gmx::ArrayRef<gmx::RVec>       forceOutput)
+{
+    impl_->compute(coordinateInput, box, forceOutput);
+}
+
+//! Compute forces and virial tensor
+void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
+                                      const Box&                     box,
+                                      gmx::ArrayRef<gmx::RVec>       forceOutput,
+                                      gmx::ArrayRef<real>            virialOutput)
+{
+    impl_->compute(coordinateInput, box, forceOutput, virialOutput);
+}
+
+//! Compute forces, virial tensor and potential energies
+void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
+                                      const Box&                     box,
+                                      gmx::ArrayRef<gmx::RVec>       forceOutput,
+                                      gmx::ArrayRef<real>            virialOutput,
+                                      gmx::ArrayRef<real>            energyOutput)
+{
+    impl_->compute(coordinateInput, box, forceOutput, virialOutput, energyOutput);
+}
+
+std::unique_ptr<GmxNBForceCalculatorCpu> setupGmxForceCalculatorCpu(const Topology&        topology,
+                                                                    const NBKernelOptions& options)
+{
+    std::vector<real> nonBondedParameters = createNonBondedParameters(
+            topology.getParticleTypes(), topology.getNonBondedInteractionMap());
+
+    std::vector<int64_t> particleInteractionFlags = createParticleInfoAllVdw(topology.numParticles());
+
+    return std::make_unique<GmxNBForceCalculatorCpu>(topology.getParticleTypeIdOfAllParticles(),
+                                                     nonBondedParameters,
+                                                     topology.getCharges(),
+                                                     particleInteractionFlags,
+                                                     topology.exclusionLists().ListRanges,
+                                                     topology.exclusionLists().ListElements,
+                                                     options);
+}
+
+} // namespace nblib
diff --git a/api/nblib/gmxcalculatorcpu.h b/api/nblib/gmxcalculatorcpu.h
new file mode 100644 (file)
index 0000000..2fec548
--- /dev/null
@@ -0,0 +1,112 @@
+/*
+ * 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.
+ */
+/*! \libinternal \file
+ * \brief
+ * Implements a force calculator based on GROMACS data structures.
+ *
+ * Intended for internal use inside the ForceCalculator.
+ *
+ * \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_GMXCALCULATORCPU_H
+#define NBLIB_GMXCALCULATORCPU_H
+
+#include <memory>
+#include <vector>
+
+#include "nblib/box.h"
+#include "nblib/vector.h"
+
+namespace gmx
+{
+template<typename T>
+class ArrayRef;
+} // namespace gmx
+
+namespace nblib
+{
+struct NBKernelOptions;
+class Topology;
+
+class GmxNBForceCalculatorCpu final
+{
+public:
+    GmxNBForceCalculatorCpu(gmx::ArrayRef<int>     particleTypeIdOfAllParticles,
+                            gmx::ArrayRef<real>    nonBondedParams,
+                            gmx::ArrayRef<real>    charges,
+                            gmx::ArrayRef<int64_t> particleInteractionFlags,
+                            gmx::ArrayRef<int>     exclusionRanges,
+                            gmx::ArrayRef<int>     exclusionElements,
+                            const NBKernelOptions& options);
+
+    ~GmxNBForceCalculatorCpu();
+
+    //! calculates a new pair list based on new coordinates (for every NS step)
+    void updatePairlist(gmx::ArrayRef<const gmx::RVec> coordinates, const Box& box);
+
+    //! Compute forces and return
+    void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
+                 const Box&                     box,
+                 gmx::ArrayRef<gmx::RVec>       forceOutput);
+
+    //! Compute forces and virial tensor
+    void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
+                 const Box&                     box,
+                 gmx::ArrayRef<gmx::RVec>       forceOutput,
+                 gmx::ArrayRef<real>            virialOutput);
+
+    //! Compute forces, virial tensor and potential energies
+    void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
+                 const Box&                     box,
+                 gmx::ArrayRef<gmx::RVec>       forceOutput,
+                 gmx::ArrayRef<real>            virialOutput,
+                 gmx::ArrayRef<real>            energyOutput);
+
+private:
+    //! Private implementation
+    class CpuImpl;
+    std::unique_ptr<CpuImpl> impl_;
+};
+
+std::unique_ptr<GmxNBForceCalculatorCpu> setupGmxForceCalculatorCpu(const Topology&        topology,
+                                                                    const NBKernelOptions& options);
+
+} // namespace nblib
+
+#endif // NBLIB_GMXCALCULATORCPU_H
index 636c92bc47bb41ba913a2639d25144be7380f528..55f3620c08daab3c3fba5e468f64be244560bca6 100644 (file)
@@ -47,8 +47,7 @@
 
 #include "nblib/basicdefinitions.h"
 #include "nblib/box.h"
-#include "nblib/gmxcalculator.h"
-#include "nblib/gmxsetup.h"
+#include "nblib/gmxcalculatorcpu.h"
 #include "nblib/integrator.h"
 #include "nblib/interactions.h"
 #include "nblib/kerneloptions.h"
@@ -56,6 +55,7 @@
 #include "nblib/listed_forces/definitions.h"
 #include "nblib/listed_forces/calculator.h"
 #include "nblib/molecules.h"
+#include "nblib/nbnxmsetuphelpers.h"
 #include "nblib/particlesequencer.h"
 #include "nblib/particletype.h"
 #include "nblib/simulationstate.h"
index 027a0fdeb1d4018e55a498ea9262b3fac4afeb2f..160612f8e8228fd2d170d962b7924878074c2b25 100644 (file)
@@ -167,7 +167,7 @@ std::vector<real> createNonBondedParameters(const std::vector<ParticleType>& par
     return nonbondedParameters;
 }
 
-gmx::StepWorkload createStepWorkload([[maybe_unused]] const NBKernelOptions& options)
+gmx::StepWorkload createStepWorkload()
 {
     gmx::StepWorkload stepWorkload;
     stepWorkload.computeForces          = true;
index 9c3ee71f4697b764f7bef79cf70823c018491b75..2f22bd825b80699aca222bce185f8bdb4b94c6b0 100644 (file)
@@ -99,7 +99,7 @@ std::vector<real> createNonBondedParameters(const std::vector<ParticleType>& par
                                             const NonBondedInteractionMap& nonBondedInteractionMap);
 
 //! Create a step work object
-gmx::StepWorkload createStepWorkload(const NBKernelOptions& options);
+gmx::StepWorkload createStepWorkload();
 
 //! Computes the Ewald splitting coefficient for Coulomb
 real ewaldCoeff(real ewald_rtol, real pairlistCutoff);
index d979c4944a399b787bfe2779ac4cf3e873a36415..3daa462ba95a74b15772ee3413ed0873d6a3e96d 100644 (file)
@@ -107,15 +107,16 @@ int main()
     // Some performance flags can be set a run time
     options.nbnxmSimd = nblib::SimdKernels::SimdNo;
     // The force calculator contains all the data needed to compute forces.
-    auto forceCalculator = nblib::setupGmxForceCalculator(
-            simState.topology(), simState.coordinates(), simState.box(), options);
+    auto forceCalculator = nblib::setupGmxForceCalculatorCpu(simState.topology(), options);
+    // build the pairlist
+    forceCalculator->updatePairlist(simState.coordinates(), simState.box());
     // Integration requires masses, positions, and forces
     nblib::LeapFrog integrator(simState.topology(), simState.box());
     // Print some diagnostic info
     printf("initial forces on particle 0: x %4f y %4f z %4f\n", forces[0][0], forces[0][1], forces[0][2]);
     // The forces are computed for the user
     gmx::ArrayRef<nblib::Vec3> userForces(simState.forces());
-    forceCalculator->compute(simState.coordinates(), userForces);
+    forceCalculator->compute(simState.coordinates(), simState.box(), userForces);
     // Print some diagnostic info
     printf("  final forces on particle 0: x %4f y %4f z %4f\n",
            userForces[0][0],
index 4dbf3b46d702ab80646d4327324a42d689cb4879..06bc92b6680341e5050e4fe11b7e3da2cfd3b7e5 100644 (file)
@@ -174,8 +174,10 @@ int main()
     SimulationState simulationState(coordinates, velocities, forces, box, topology);
 
     // The non-bonded force calculator contains all the data needed to compute forces
-    auto forceCalculator = nblib::setupGmxForceCalculator(
-            simulationState.topology(), simulationState.coordinates(), simulationState.box(), options);
+    auto forceCalculator = setupGmxForceCalculatorCpu(simulationState.topology(), options);
+
+    // build the pair list
+    forceCalculator->updatePairlist(simulationState.coordinates(), simulationState.box());
 
     // The listed force calculator is also initialized with the required arguments
     ListedForceCalculator listedForceCalculator(
@@ -183,7 +185,7 @@ int main()
 
     // Integrator is initialized with an array of inverse masses (constructed from topology) and
     // the bounding box
-    LeapFrog integrator(topology, box);
+    LeapFrog integrator(simulationState.topology(), simulationState.box());
 
     // Print some diagnostic info
     printf("initial position of particle 0: x %4f y %4f z %4f\n",
@@ -198,7 +200,8 @@ int main()
     {
         zeroCartesianArray(simulationState.forces());
 
-        forceCalculator->compute(simulationState.coordinates(), simulationState.forces());
+        forceCalculator->compute(
+                simulationState.coordinates(), simulationState.box(), simulationState.forces());
 
         listedForceCalculator.compute(simulationState.coordinates(), simulationState.forces());
 
diff --git a/api/nblib/systemdescription.h b/api/nblib/systemdescription.h
new file mode 100644 (file)
index 0000000..0d2a3e0
--- /dev/null
@@ -0,0 +1,130 @@
+/*
+ * 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 simulation box
+ *
+ * \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_SYSTEMDESCRIPTION_H
+#define NBLIB_SYSTEMDESCRIPTION_H
+
+#include "gromacs/ewald/ewald_utils.h"
+#include "gromacs/gmxlib/nrnb.h"
+#include "gromacs/gpu_utils/device_stream_manager.h"
+#include "gromacs/mdlib/rf_util.h"
+#include "gromacs/mdtypes/enerdata.h"
+#include "gromacs/mdtypes/forcerec.h"
+#include "gromacs/mdtypes/interaction_const.h"
+#include "gromacs/mdtypes/simulation_workload.h"
+#include "gromacs/nbnxm/gpu_data_mgmt.h"
+#include "gromacs/nbnxm/nbnxm_gpu.h"
+#include "gromacs/nbnxm/atomdata.h"
+#include "gromacs/nbnxm/nbnxm.h"
+#include "gromacs/nbnxm/pairlistset.h"
+#include "gromacs/nbnxm/pairlistsets.h"
+#include "gromacs/nbnxm/pairsearch.h"
+#include "gromacs/timing/wallcycle.h"
+#include "gromacs/utility/listoflists.h"
+#include "gromacs/utility/range.h"
+
+namespace nblib
+{
+
+//! \brief client-side provided system description data
+struct SystemDescription
+{
+    SystemDescription() = default;
+
+    SystemDescription(gmx::ArrayRef<int>     particleTypeIdOfAllParticles,
+                      gmx::ArrayRef<real>    nonBondedParams,
+                      gmx::ArrayRef<real>    charges,
+                      gmx::ArrayRef<int64_t> particleInteractionFlags)
+    {
+        std::array inputSizes{ particleTypeIdOfAllParticles.size(),
+                               charges.size(),
+                               particleInteractionFlags.size() };
+        if (static_cast<unsigned long>(std::count(begin(inputSizes), end(inputSizes), inputSizes[0]))
+            != inputSizes.size())
+        {
+            throw InputException("input array size inconsistent");
+        }
+
+        int numParticleTypes = int(std::round(std::sqrt(nonBondedParams.size() / 2)));
+        if (2 * numParticleTypes * numParticleTypes != int(nonBondedParams.size()))
+        {
+            throw InputException("Wrong size of nonBondedParams");
+        }
+
+        numParticles_     = particleTypeIdOfAllParticles.size();
+        numParticleTypes_ = numParticleTypes;
+
+        particleTypeIdOfAllParticles_ = std::vector<int>(particleTypeIdOfAllParticles.begin(),
+                                                         particleTypeIdOfAllParticles.end());
+
+        nonBondedParams_ = std::vector<real>(nonBondedParams.begin(), nonBondedParams.end());
+        charges_         = std::vector<real>(charges.begin(), charges.end());
+        particleInfo_ =
+                std::vector<int64_t>(particleInteractionFlags.begin(), particleInteractionFlags.end());
+    }
+
+    //! number of particles
+    size_t numParticles_{ 0 };
+
+    //! number of particle types
+    size_t numParticleTypes_{ 0 };
+
+    //! particle type id of all particles
+    std::vector<int> particleTypeIdOfAllParticles_;
+
+    //! Storage for parameters for short range interactions.
+    std::vector<real> nonBondedParams_;
+
+    //! electrostatic charges
+    std::vector<real> charges_;
+
+    //! flag for each particle to set LJ and Q interactions
+    std::vector<int64_t> particleInfo_;
+
+    //! Legacy matrix for box
+    Box box_{ 0 };
+};
+
+
+} // namespace nblib
+#endif // NBLIB_SYSTEMDESCRIPTION_H
index 141ce590234b6c155391b0f1a3c64f432c86d34c..e133bc96fcff87754c9430c5a1c8bdee85b889c8 100644 (file)
  */
 #include <gtest/gtest.h>
 
-#include "nblib/gmxsetup.h"
+#include "nblib/gmxcalculatorcpu.h"
 #include "nblib/kerneloptions.h"
 #include "nblib/simulationstate.h"
+#include "nblib/tests/testhelpers.h"
 #include "nblib/tests/testsystems.h"
 #include "gromacs/utility/arrayref.h"
 
@@ -61,23 +62,72 @@ TEST(NBlibTest, GmxForceCalculatorCanCompute)
     SimulationState             simState = argonSystemBuilder.setupSimulationState();
     NBKernelOptions             options  = NBKernelOptions();
     options.nbnxmSimd                    = SimdKernels::SimdNo;
-    std::unique_ptr<GmxForceCalculator> gmxForceCalculator = nblib::setupGmxForceCalculator(
-            simState.topology(), simState.coordinates(), simState.box(), options);
-    ;
-    EXPECT_NO_THROW(gmxForceCalculator->compute(simState.coordinates(), simState.forces()));
+
+    std::unique_ptr<GmxNBForceCalculatorCpu> gmxForceCalculator =
+            setupGmxForceCalculatorCpu(simState.topology(), options);
+    gmxForceCalculator->updatePairlist(simState.coordinates(), simState.box());
+
+    EXPECT_NO_THROW(gmxForceCalculator->compute(simState.coordinates(), simState.box(), simState.forces()));
 }
 
-TEST(NBlibTest, CanSetupStepWorkload)
+TEST(NBlibTest, ArgonVirialsAreCorrect)
 {
-    NBKernelOptions options;
-    EXPECT_NO_THROW(NbvSetupUtil{}.setupStepWorkload(options));
+    ArgonSimulationStateBuilder argonSystemBuilder(fftypes::OPLSA);
+    SimulationState             simState = argonSystemBuilder.setupSimulationState();
+    NBKernelOptions             options  = NBKernelOptions();
+    options.nbnxmSimd                    = SimdKernels::SimdNo;
+    std::unique_ptr<GmxNBForceCalculatorCpu> gmxForceCalculator =
+            setupGmxForceCalculatorCpu(simState.topology(), options);
+    gmxForceCalculator->updatePairlist(simState.coordinates(), simState.box());
+
+    std::vector<real> virialArray(9, 0.0);
+
+    gmxForceCalculator->compute(simState.coordinates(), simState.box(), simState.forces(), virialArray);
+
+    RefDataChecker virialsOutputTest(1e-7);
+    virialsOutputTest.testArrays<real>(virialArray, "Virials");
 }
 
-TEST(NBlibTest, GmxForceCalculatorCanSetupInteractionConst)
+TEST(NBlibTest, ArgonEnergiesAreCorrect)
 {
-    NBKernelOptions options;
-    EXPECT_NO_THROW(NbvSetupUtil{}.setupInteractionConst(options));
+    ArgonSimulationStateBuilder argonSystemBuilder(fftypes::OPLSA);
+    SimulationState             simState = argonSystemBuilder.setupSimulationState();
+    NBKernelOptions             options  = NBKernelOptions();
+    options.nbnxmSimd                    = SimdKernels::SimdNo;
+    std::unique_ptr<GmxNBForceCalculatorCpu> gmxForceCalculator =
+            setupGmxForceCalculatorCpu(simState.topology(), options);
+    gmxForceCalculator->updatePairlist(simState.coordinates(), simState.box());
+
+    // number of energy kinds is 5: COULSR, LJSR, BHAMSR, COUL14, LJ14,
+    std::vector<real> energies(5, 0.0);
+
+    gmxForceCalculator->compute(
+            simState.coordinates(), simState.box(), simState.forces(), gmx::ArrayRef<real>{}, energies);
+
+    RefDataChecker energiesOutputTest(5e-5);
+    energiesOutputTest.testArrays<real>(energies, "Argon energies");
+}
+
+TEST(NBlibTest, SpcMethanolEnergiesAreCorrect)
+{
+    SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder;
+    SimulationState                   simState = spcMethanolSystemBuilder.setupSimulationState();
+    NBKernelOptions                   options  = NBKernelOptions();
+    options.nbnxmSimd                          = SimdKernels::SimdNo;
+    std::unique_ptr<GmxNBForceCalculatorCpu> gmxForceCalculator =
+            setupGmxForceCalculatorCpu(simState.topology(), options);
+    gmxForceCalculator->updatePairlist(simState.coordinates(), simState.box());
+
+    // number of energy kinds is 5: COULSR, LJSR, BHAMSR, COUL14, LJ14,
+    std::vector<real> energies(5, 0.0);
+
+    gmxForceCalculator->compute(
+            simState.coordinates(), simState.box(), simState.forces(), gmx::ArrayRef<real>{}, energies);
+
+    RefDataChecker energiesOutputTest(5e-5);
+    energiesOutputTest.testArrays<real>(energies, "SPC-methanol energies");
 }
+
 } // namespace
 } // namespace test
 } // namespace nblib
index 15163e1384d60efdbd8c39f67cf64176440a7214..96c874a39361ebd49dbc23e75a830e4b940c2b2c 100644 (file)
@@ -45,8 +45,7 @@
 #include <gtest/gtest.h>
 
 #include "gromacs/topology/exclusionblocks.h"
-#include "nblib/gmxcalculator.h"
-#include "nblib/gmxsetup.h"
+#include "nblib/gmxcalculatorcpu.h"
 #include "nblib/integrator.h"
 #include "nblib/tests/testhelpers.h"
 #include "nblib/tests/testsystems.h"
@@ -68,12 +67,14 @@ TEST(NBlibTest, SpcMethanolForcesAreCorrect)
 
     SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder;
 
-    auto simState        = spcMethanolSystemBuilder.setupSimulationState();
-    auto forceCalculator = nblib::setupGmxForceCalculator(
-            simState.topology(), simState.coordinates(), simState.box(), options);
+    auto simState = spcMethanolSystemBuilder.setupSimulationState();
+
+    auto forceCalculator = setupGmxForceCalculatorCpu(simState.topology(), options);
+
+    forceCalculator->updatePairlist(simState.coordinates(), simState.box());
 
     gmx::ArrayRef<Vec3> forces(simState.forces());
-    ASSERT_NO_THROW(forceCalculator->compute(simState.coordinates(), forces));
+    ASSERT_NO_THROW(forceCalculator->compute(simState.coordinates(), simState.box(), forces));
 
     RefDataChecker forcesOutputTest(5e-5);
     forcesOutputTest.testArrays<Vec3>(forces, "SPC-methanol forces");
@@ -86,12 +87,14 @@ TEST(NBlibTest, ExpectedNumberOfForces)
 
     SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder;
 
-    auto simState        = spcMethanolSystemBuilder.setupSimulationState();
-    auto forceCalculator = nblib::setupGmxForceCalculator(
-            simState.topology(), simState.coordinates(), simState.box(), options);
+    auto simState = spcMethanolSystemBuilder.setupSimulationState();
+
+    auto forceCalculator = setupGmxForceCalculatorCpu(simState.topology(), options);
+
+    forceCalculator->updatePairlist(simState.coordinates(), simState.box());
 
     gmx::ArrayRef<Vec3> forces(simState.forces());
-    forceCalculator->compute(simState.coordinates(), forces);
+    forceCalculator->compute(simState.coordinates(), simState.box(), forces);
     EXPECT_EQ(simState.topology().numParticles(), forces.size());
 }
 
@@ -103,18 +106,19 @@ TEST(NBlibTest, CanIntegrateSystem)
 
     SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder;
 
-    auto simState        = spcMethanolSystemBuilder.setupSimulationState();
-    auto forceCalculator = nblib::setupGmxForceCalculator(
-            simState.topology(), simState.coordinates(), simState.box(), options);
+    auto simState = spcMethanolSystemBuilder.setupSimulationState();
+
+    auto forceCalculator = setupGmxForceCalculatorCpu(simState.topology(), options);
+
+    forceCalculator->updatePairlist(simState.coordinates(), simState.box());
 
     LeapFrog integrator(simState.topology(), simState.box());
 
     for (int iter = 0; iter < options.numIterations; iter++)
     {
         gmx::ArrayRef<Vec3> forces(simState.forces());
-        forceCalculator->compute(simState.coordinates(), simState.forces());
-        EXPECT_NO_THROW(integrator.integrate(
-                1.0, simState.coordinates(), simState.velocities(), simState.forces()));
+        forceCalculator->compute(simState.coordinates(), simState.box(), forces);
+        EXPECT_NO_THROW(integrator.integrate(1.0, simState.coordinates(), simState.velocities(), forces));
     }
 }
 
@@ -139,15 +143,18 @@ TEST(NBlibTest, UpdateChangesForces)
 
     SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder;
 
-    auto simState        = spcMethanolSystemBuilder.setupSimulationState();
-    auto forceCalculator = nblib::setupGmxForceCalculator(
-            simState.topology(), simState.coordinates(), simState.box(), options);
+    auto simState = spcMethanolSystemBuilder.setupSimulationState();
+
+    auto forceCalculator = setupGmxForceCalculatorCpu(simState.topology(), options);
+
+    forceCalculator->updatePairlist(simState.coordinates(), simState.box());
 
     LeapFrog integrator(simState.topology(), simState.box());
 
     // step 1
     gmx::ArrayRef<Vec3> forces(simState.forces());
-    forceCalculator->compute(simState.coordinates(), simState.forces());
+
+    forceCalculator->compute(simState.coordinates(), simState.box(), simState.forces());
 
     // copy computed forces to another array
     std::vector<Vec3> forces_1(forces.size());
@@ -157,7 +164,7 @@ TEST(NBlibTest, UpdateChangesForces)
     zeroCartesianArray(forces);
 
     // check if forces change without update step
-    forceCalculator->compute(simState.coordinates(), forces);
+    forceCalculator->compute(simState.coordinates(), simState.box(), forces);
 
     // check if forces change without update
     for (size_t i = 0; i < forces_1.size(); i++)
@@ -175,7 +182,8 @@ TEST(NBlibTest, UpdateChangesForces)
     zeroCartesianArray(forces);
 
     // step 2
-    forceCalculator->compute(simState.coordinates(), forces);
+    forceCalculator->compute(simState.coordinates(), simState.box(), forces);
+
     std::vector<Vec3> forces_2(forces.size());
     std::copy(forces.begin(), forces.end(), begin(forces_2));
 
@@ -197,12 +205,14 @@ TEST(NBlibTest, ArgonOplsaForcesAreCorrect)
 
     ArgonSimulationStateBuilder argonSystemBuilder(fftypes::OPLSA);
 
-    auto simState        = argonSystemBuilder.setupSimulationState();
-    auto forceCalculator = nblib::setupGmxForceCalculator(
-            simState.topology(), simState.coordinates(), simState.box(), options);
+    auto simState = argonSystemBuilder.setupSimulationState();
+
+    auto forceCalculator = setupGmxForceCalculatorCpu(simState.topology(), options);
+
+    forceCalculator->updatePairlist(simState.coordinates(), simState.box());
 
     gmx::ArrayRef<Vec3> testForces(simState.forces());
-    forceCalculator->compute(simState.coordinates(), simState.forces());
+    forceCalculator->compute(simState.coordinates(), simState.box(), simState.forces());
 
     RefDataChecker forcesOutputTest(1e-7);
     forcesOutputTest.testArrays<Vec3>(testForces, "Argon forces");
@@ -216,14 +226,16 @@ TEST(NBlibTest, ArgonGromos43A1ForcesAreCorrect)
 
     ArgonSimulationStateBuilder argonSystemBuilder(fftypes::GROMOS43A1);
 
-    auto simState        = argonSystemBuilder.setupSimulationState();
-    auto forceCalculator = nblib::setupGmxForceCalculator(
-            simState.topology(), simState.coordinates(), simState.box(), options);
+    auto simState = argonSystemBuilder.setupSimulationState();
+
+    auto forceCalculator = setupGmxForceCalculatorCpu(simState.topology(), options);
+
+    forceCalculator->updatePairlist(simState.coordinates(), simState.box());
 
     gmx::ArrayRef<Vec3> testForces(simState.forces());
-    forceCalculator->compute(simState.coordinates(), simState.forces());
+    forceCalculator->compute(simState.coordinates(), simState.box(), simState.forces());
 
-    RefDataChecker forcesOutputTest(1e-7);
+    RefDataChecker forcesOutputTest;
     forcesOutputTest.testArrays<Vec3>(testForces, "Argon forces");
 }
 
diff --git a/api/nblib/tests/refdata/NBlibTest_ArgonEnergiesAreCorrect.xml b/api/nblib/tests/refdata/NBlibTest_ArgonEnergiesAreCorrect.xml
new file mode 100644 (file)
index 0000000..08bfa0d
--- /dev/null
@@ -0,0 +1,12 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <Sequence Name="Argon energies">
+    <Int Name="Length">5</Int>
+    <Real>0</Real>
+    <Real>-1.0175924</Real>
+    <Real>0</Real>
+    <Real>0</Real>
+    <Real>0</Real>
+  </Sequence>
+</ReferenceData>
diff --git a/api/nblib/tests/refdata/NBlibTest_ArgonVirialsAreCorrect.xml b/api/nblib/tests/refdata/NBlibTest_ArgonVirialsAreCorrect.xml
new file mode 100644 (file)
index 0000000..d648ef1
--- /dev/null
@@ -0,0 +1,16 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+    <Sequence Name="Virials">
+        <Int Name="Length">9</Int>
+        <Real>0.094908442751139188</Real>
+        <Real>0.25238615516784413</Real>
+        <Real>0.02601194356883077</Real>
+        <Real>0.25238615516784413</Real>
+        <Real>0.67116022003893372</Real>
+        <Real>0.069172501786742557</Real>
+        <Real>0.026011943568830798</Real>
+        <Real>0.069172501786742391</Real>
+        <Real>0.0071291993484943611</Real>
+    </Sequence>
+</ReferenceData>
\ No newline at end of file