Fix nblib pairlist update function

Previously the function put_atoms_in_box was called only by the
nblib integrator, or when constructing a system via a SimulationState
object. In the case of the integrator, this is a needless performance
degradation. When using an nb force calculator without
first putting atoms in the box, this could lead to a cryptic
error from nbnxm grid search failure. Both of these issues are
rectified with this change, which also adds a member to the
non-bonded force calculator to hold the requested number of OMP
threads to use in a call to put_atoms_in_box_omp, which is faster
than the non OMP version.

diff --git a/api/nblib/gmxbackenddata.h b/api/nblib/gmxbackenddata.h
index 6794bb47debd747241e8fc99336e816c9e238330..d9528f4e35da177c79a73f2ec9d9158b00e92153 100644 (file)
--- a/api/nblib/gmxbackenddata.h
+++ b/api/nblib/gmxbackenddata.h
@@ -77,7 +77,8 @@ public:
     GmxBackendData(const NBKernelOptions& options,
                    int                    numEnergyGroups,
                    gmx::ArrayRef<int>     exclusionRanges,
-                   gmx::ArrayRef<int>     exclusionElements)
+                   gmx::ArrayRef<int>     exclusionElements) :
+        numThreads_(options.numOpenMPThreads)
     {
         // Set hardware params from the execution context
         setGmxNonBondedNThreads(options.numOpenMPThreads);
@@ -120,6 +121,9 @@ public:
     //! Non-bonded flop counter; currently only needed as an argument for dispatchNonbondedKernel
     t_nrnb nrnb_;
 
+    //! Number of OpenMP threads to use
+    int numThreads_;
+
     //! Keep track of whether updatePairlist has been called at least once
     bool updatePairlistCalled{ false };
 };

diff --git a/api/nblib/gmxcalculatorcpu.cpp b/api/nblib/gmxcalculatorcpu.cpp
index a976e5d370c9faee432d634376642f5660c05e44..5d8fa76fcb25d2eec9831c182d626cafbfc7e366 100644 (file)
--- a/api/nblib/gmxcalculatorcpu.cpp
+++ b/api/nblib/gmxcalculatorcpu.cpp
@@ -75,7 +75,7 @@ public:
             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);
+    void updatePairlist(gmx::ArrayRef<gmx::RVec> coordinates, const Box& box);
 
     //! Compute forces and return
     void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,
@@ -120,8 +120,7 @@ GmxNBForceCalculatorCpu::CpuImpl::CpuImpl(gmx::ArrayRef<int>     particleTypeIdO
                                    system_.nonBondedParams_);
 }
 
-void GmxNBForceCalculatorCpu::CpuImpl::updatePairlist(gmx::ArrayRef<const gmx::RVec> coordinates,
-                                                      const Box&                     box)
+void GmxNBForceCalculatorCpu::CpuImpl::updatePairlist(gmx::ArrayRef<gmx::RVec> coordinates, const Box& box)
 {
     if (coordinates.size() != system_.numParticles_)
     {
@@ -143,6 +142,9 @@ void GmxNBForceCalculatorCpu::CpuImpl::updatePairlist(gmx::ArrayRef<const gmx::R
 
     const real particleDensity = static_cast<real>(coordinates.size()) / det(legacyBox);
 
+    // If particles are too far outside the box, the grid setup can fail
+    put_atoms_in_box_omp(PbcType::Xyz, box.legacyMatrix(), coordinates, backend_.numThreads_);
+
     // Put particles on a grid based on bounds specified by the box
     nbnxn_put_on_grid(backend_.nbv_.get(),
                       legacyBox,
@@ -288,7 +290,7 @@ GmxNBForceCalculatorCpu::GmxNBForceCalculatorCpu(gmx::ArrayRef<int>  particleTyp
 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)
+void GmxNBForceCalculatorCpu::updatePairlist(gmx::ArrayRef<gmx::RVec> coordinates, const Box& box)
 {
     impl_->updatePairlist(coordinates, box);
 }

diff --git a/api/nblib/gmxcalculatorcpu.h b/api/nblib/gmxcalculatorcpu.h
index efc079896ef7d7baccbcfa5f9b39328b24a146c0..81783db7521e5dc5c94d2c6069389b63d002aa3d 100644 (file)
--- a/api/nblib/gmxcalculatorcpu.h
+++ b/api/nblib/gmxcalculatorcpu.h
@@ -79,7 +79,7 @@ public:
     ~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);
+    void updatePairlist(gmx::ArrayRef<gmx::RVec> coordinates, const Box& box);
 
     //! Compute forces and return
     void compute(gmx::ArrayRef<const gmx::RVec> coordinateInput,

diff --git a/api/nblib/integrator.cpp b/api/nblib/integrator.cpp
index a80338d5ffe939b093b6c983d786dec750a5b7e6..5168850eea5ed644ad4130be5f8768341355f701 100644 (file)
--- a/api/nblib/integrator.cpp
+++ b/api/nblib/integrator.cpp
@@ -60,6 +60,11 @@ LeapFrog::LeapFrog(const Topology& topology, const Box& box) : box_(box)
     }
 }
 
+LeapFrog::LeapFrog(gmx::ArrayRef<const real> inverseMasses, const Box& box) :
+    inverseMasses_(inverseMasses.begin(), inverseMasses.end()), box_(box)
+{
+}
+
 void LeapFrog::integrate(const real dt, gmx::ArrayRef<Vec3> x, gmx::ArrayRef<Vec3> v, gmx::ArrayRef<const Vec3> f)
 {
     for (size_t i = 0; i < x.size(); i++)
@@ -70,7 +75,6 @@ void LeapFrog::integrate(const real dt, gmx::ArrayRef<Vec3> x, gmx::ArrayRef<Vec
             x[i][dim] += v[i][dim] * dt;
         }
     }
-    put_atoms_in_box(PbcType::Xyz, box_.legacyMatrix(), x);
 }
 
 } // namespace nblib

diff --git a/api/nblib/integrator.h b/api/nblib/integrator.h
index f72d6b46d1fecc2059e7cba2479afcace1ca0a5d..02f3a907bf68e848f391bbb026d2f85d2106e8c3 100644 (file)
--- a/api/nblib/integrator.h
+++ b/api/nblib/integrator.h
@@ -74,6 +74,13 @@ public:
      */
     LeapFrog(const Topology& topology, const Box& box);
 
+    /*! \brief Constructor.
+     *
+     * \param[in] masses  List of inverse masses.
+     * \param[in] box     Box object for ensuring that coordinates remain within bounds
+     */
+    LeapFrog(gmx::ArrayRef<const real> inverseMasses, const Box& box);
+
     /*! \brief Integrate
      *
      * Integrates the equation of motion using Leap-Frog algorithm.