Remove mdatoms from forceprovider call signature
authorPaul Bauer <paul.bauer.q@gmail.com>
Fri, 19 Mar 2021 11:11:49 +0000 (11:11 +0000)
committerMark Abraham <mark.j.abraham@gmail.com>
Fri, 19 Mar 2021 11:11:49 +0000 (11:11 +0000)
This makes it possible to test the forceprovider without having to build
a partial mdatoms datastructure first. Needed for future refactoring of
mdatoms.

src/gromacs/applied_forces/densityfitting/densityfittingamplitudelookup.cpp
src/gromacs/applied_forces/densityfitting/densityfittingamplitudelookup.h
src/gromacs/applied_forces/densityfitting/densityfittingforceprovider.cpp
src/gromacs/applied_forces/densityfitting/tests/densityfittingamplitudelookup.cpp
src/gromacs/applied_forces/electricfield.cpp
src/gromacs/applied_forces/tests/electricfield.cpp
src/gromacs/mdlib/sim_util.cpp
src/gromacs/mdtypes/iforceprovider.h
src/gromacs/restraint/restraintmdmodule.cpp

index 32dfd24aaefdbcd2f993f9521d2f66dbc572a8c1..a3bcbc1fb23b5b123c9792fb787cf03e6b66356a 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2019,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.
@@ -46,7 +46,6 @@
 #include <algorithm>
 #include <vector>
 
-#include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/utility/arrayref.h"
 
 namespace gmx
@@ -62,8 +61,10 @@ public:
     DensityFittingAmplitudeLookupImpl(const DensityFittingAmplitudeLookupImpl&) = default;
     virtual ~DensityFittingAmplitudeLookupImpl()                                = default;
 
-    virtual const std::vector<real>& operator()(const t_mdatoms& atoms, ArrayRef<const int> localIndex) = 0;
-    virtual std::unique_ptr<DensityFittingAmplitudeLookupImpl> clone() = 0;
+    virtual const std::vector<real>& operator()(ArrayRef<const real> /*chargeA*/,
+                                                ArrayRef<const real> /*massT*/,
+                                                ArrayRef<const int> localIndex) = 0;
+    virtual std::unique_ptr<DensityFittingAmplitudeLookupImpl> clone()          = 0;
 };
 
 namespace
@@ -75,7 +76,9 @@ public:
     UnitAmplitudes(const UnitAmplitudes&) = default;
     ~UnitAmplitudes() override            = default;
     std::unique_ptr<DensityFittingAmplitudeLookupImpl> clone() override;
-    const std::vector<real>& operator()(const t_mdatoms& atoms, ArrayRef<const int> localIndex) override;
+    const std::vector<real>&                           operator()(ArrayRef<const real> /*chargeA*/,
+                                        ArrayRef<const real> /*massT*/,
+                                        ArrayRef<const int> localIndex) override;
 
 private:
     std::vector<real> amplitude_;
@@ -86,7 +89,9 @@ std::unique_ptr<DensityFittingAmplitudeLookupImpl> UnitAmplitudes::clone()
     return std::make_unique<UnitAmplitudes>(*this);
 };
 
-const std::vector<real>& UnitAmplitudes::operator()(const t_mdatoms& /*atoms*/, ArrayRef<const int> localIndex)
+const std::vector<real>& UnitAmplitudes::operator()(ArrayRef<const real> /*chargeA*/,
+                                                    ArrayRef<const real> /*massT*/,
+                                                    ArrayRef<const int> localIndex)
 {
     if (amplitude_.size() != localIndex.size())
     {
@@ -103,7 +108,9 @@ public:
     ChargesAsAmplitudes(const ChargesAsAmplitudes&) = default;
     ~ChargesAsAmplitudes() override                 = default;
     std::unique_ptr<DensityFittingAmplitudeLookupImpl> clone() override;
-    const std::vector<real>& operator()(const t_mdatoms& atoms, ArrayRef<const int> localIndex) override;
+    const std::vector<real>&                           operator()(ArrayRef<const real> chargeA,
+                                        ArrayRef<const real> /*massT*/,
+                                        ArrayRef<const int> localIndex) override;
 
 private:
     std::vector<real> amplitude_;
@@ -114,7 +121,9 @@ std::unique_ptr<DensityFittingAmplitudeLookupImpl> ChargesAsAmplitudes::clone()
     return std::make_unique<ChargesAsAmplitudes>(*this);
 };
 
-const std::vector<real>& ChargesAsAmplitudes::operator()(const t_mdatoms& atoms, ArrayRef<const int> localIndex)
+const std::vector<real>& ChargesAsAmplitudes::operator()(ArrayRef<const real> chargeA,
+                                                         ArrayRef<const real> /*massT*/,
+                                                         ArrayRef<const int> localIndex)
 {
     if (amplitude_.size() != localIndex.size())
     {
@@ -124,7 +133,7 @@ const std::vector<real>& ChargesAsAmplitudes::operator()(const t_mdatoms& atoms,
     std::transform(std::begin(localIndex),
                    std::end(localIndex),
                    std::begin(amplitude_),
-                   [&atoms](gmx::index index) { return atoms.chargeA[index]; });
+                   [&chargeA](gmx::index index) { return chargeA[index]; });
     return amplitude_;
 }
 
@@ -135,7 +144,9 @@ public:
     MassesAsAmplitudes(const MassesAsAmplitudes&) = default;
     ~MassesAsAmplitudes() override                = default;
     std::unique_ptr<DensityFittingAmplitudeLookupImpl> clone() override;
-    const std::vector<real>& operator()(const t_mdatoms& atoms, ArrayRef<const int> localIndex) override;
+    const std::vector<real>&                           operator()(ArrayRef<const real> /*chargeA*/,
+                                        ArrayRef<const real> massT,
+                                        ArrayRef<const int>  localIndex) override;
 
 private:
     std::vector<real> amplitude_;
@@ -146,7 +157,9 @@ std::unique_ptr<DensityFittingAmplitudeLookupImpl> MassesAsAmplitudes::clone()
     return std::make_unique<MassesAsAmplitudes>(*this);
 };
 
-const std::vector<real>& MassesAsAmplitudes::operator()(const t_mdatoms& atoms, ArrayRef<const int> localIndex)
+const std::vector<real>& MassesAsAmplitudes::operator()(ArrayRef<const real> /*chargeA*/,
+                                                        ArrayRef<const real> massT,
+                                                        ArrayRef<const int>  localIndex)
 {
     if (amplitude_.size() != localIndex.size())
     {
@@ -156,7 +169,7 @@ const std::vector<real>& MassesAsAmplitudes::operator()(const t_mdatoms& atoms,
     std::transform(std::begin(localIndex),
                    std::end(localIndex),
                    std::begin(amplitude_),
-                   [&atoms](gmx::index index) { return atoms.massT[index]; });
+                   [&massT](gmx::index index) { return massT[index]; });
     return amplitude_;
 }
 
@@ -184,10 +197,11 @@ DensityFittingAmplitudeLookup::DensityFittingAmplitudeLookup(const DensityFittin
     }
 }
 
-const std::vector<real>& DensityFittingAmplitudeLookup::operator()(const t_mdatoms&    atoms,
-                                                                   ArrayRef<const int> localIndex)
+const std::vector<real>& DensityFittingAmplitudeLookup::operator()(ArrayRef<const real> chargeA,
+                                                                   ArrayRef<const real> massT,
+                                                                   ArrayRef<const int>  localIndex)
 {
-    return (*impl_)(atoms, localIndex);
+    return (*impl_)(chargeA, massT, localIndex);
 }
 
 DensityFittingAmplitudeLookup::~DensityFittingAmplitudeLookup() = default;
index 86675bcec649c882f9021e06867e2d569bcd3561..92b21f1e2b5505055631f33159c4a23c094c628d 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2019,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.
@@ -49,8 +49,6 @@
 #include "gromacs/utility/enumerationhelpers.h"
 #include "gromacs/utility/real.h"
 
-struct t_mdatoms;
-
 namespace gmx
 {
 
@@ -88,11 +86,14 @@ public:
     //! Move assignment
     DensityFittingAmplitudeLookup& operator=(DensityFittingAmplitudeLookup&& other) noexcept;
     /*! \brief Return the amplitudes for spreading atoms of a given local index.
-     * \param[in] atoms the atom information
+     * \param[in] chargeA Atom charges
+     * \param[in] massT   Atom masses.
      * \param[in] localIndex the local atom indices
      * \returns amplitudes
      */
-    const std::vector<real>& operator()(const t_mdatoms& atoms, ArrayRef<const int> localIndex);
+    const std::vector<real>& operator()(ArrayRef<const real> chargeA,
+                                        ArrayRef<const real> massT,
+                                        ArrayRef<const int>  localIndex);
 
 private:
     std::unique_ptr<DensityFittingAmplitudeLookupImpl> impl_;
index 4bd0a47df4c1fcd6bdff52e379fafc7c49fef974..eec0e3c8b45375d4c5387c28c1545ca5e3a98a58 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2019,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.
@@ -298,8 +298,8 @@ void DensityFittingForceProvider::Impl::calculateForces(const ForceProviderInput
     // spread atoms on grid
     gaussTransform_.setZero();
 
-    std::vector<real> amplitudes =
-            amplitudeLookup_(forceProviderInput.mdatoms_, localAtomSet_.localIndex());
+    std::vector<real> amplitudes = amplitudeLookup_(
+            forceProviderInput.chargeA_, forceProviderInput.massT_, localAtomSet_.localIndex());
 
     if (parameters_.normalizeDensities_)
     {
index d0eeb3a9e3df47523225265a2395f97fc85d05e5..7ffc527d22bbf8753d85be0f5524363ac3356ae6 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2019,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.
@@ -55,14 +55,6 @@ namespace gmx
 
 class DensityFittingAmplitudeLookupTest : public ::testing::Test
 {
-public:
-    DensityFittingAmplitudeLookupTest()
-    {
-        atoms_.nr      = numberOfAtoms_;
-        atoms_.massT   = masses_.data();
-        atoms_.chargeA = charges_.data();
-    }
-
 protected:
     int               numberOfAtoms_ = 3;
     std::vector<real> masses_        = { 2, 3, 4 };
@@ -74,7 +66,7 @@ protected:
 TEST_F(DensityFittingAmplitudeLookupTest, Unity)
 {
     DensityFittingAmplitudeLookup lookup(DensityFittingAmplitudeMethod::Unity);
-    const auto                    lookupResult = lookup(atoms_, lookupIndices_);
+    const auto                    lookupResult = lookup(charges_, masses_, lookupIndices_);
     EXPECT_EQ(lookupResult[0], 1);
     EXPECT_EQ(lookupResult[1], 1);
 }
@@ -82,7 +74,7 @@ TEST_F(DensityFittingAmplitudeLookupTest, Unity)
 TEST_F(DensityFittingAmplitudeLookupTest, Charge)
 {
     DensityFittingAmplitudeLookup lookup(DensityFittingAmplitudeMethod::Charge);
-    const auto                    lookupResult = lookup(atoms_, lookupIndices_);
+    const auto                    lookupResult = lookup(charges_, masses_, lookupIndices_);
     EXPECT_EQ(lookupResult[0], 30);
     EXPECT_EQ(lookupResult[1], 40);
 }
@@ -90,7 +82,7 @@ TEST_F(DensityFittingAmplitudeLookupTest, Charge)
 TEST_F(DensityFittingAmplitudeLookupTest, Masses)
 {
     DensityFittingAmplitudeLookup lookup(DensityFittingAmplitudeMethod::Mass);
-    const auto                    lookupResult = lookup(atoms_, lookupIndices_);
+    const auto                    lookupResult = lookup(charges_, masses_, lookupIndices_);
     EXPECT_EQ(lookupResult[0], 3);
     EXPECT_EQ(lookupResult[1], 4);
 }
@@ -99,7 +91,7 @@ TEST_F(DensityFittingAmplitudeLookupTest, CanCopyAssign)
 {
     DensityFittingAmplitudeLookup lookup(DensityFittingAmplitudeMethod::Unity);
     DensityFittingAmplitudeLookup lookupCopied = lookup;
-    const auto                    lookupResult = lookupCopied(atoms_, lookupIndices_);
+    const auto                    lookupResult = lookupCopied(charges_, masses_, lookupIndices_);
     EXPECT_EQ(lookupResult[0], 1);
     EXPECT_EQ(lookupResult[1], 1);
 }
@@ -108,7 +100,7 @@ TEST_F(DensityFittingAmplitudeLookupTest, CanCopyConstruct)
 {
     DensityFittingAmplitudeLookup lookup(DensityFittingAmplitudeMethod::Unity);
     DensityFittingAmplitudeLookup lookupCopied(lookup);
-    const auto                    lookupResult = lookupCopied(atoms_, lookupIndices_);
+    const auto                    lookupResult = lookupCopied(charges_, masses_, lookupIndices_);
     EXPECT_EQ(lookupResult[0], 1);
     EXPECT_EQ(lookupResult[1], 1);
 }
@@ -117,7 +109,7 @@ TEST_F(DensityFittingAmplitudeLookupTest, CanMoveAssign)
 {
     DensityFittingAmplitudeLookup lookup(DensityFittingAmplitudeMethod::Unity);
     DensityFittingAmplitudeLookup lookupCopied = std::move(lookup);
-    const auto                    lookupResult = lookupCopied(atoms_, lookupIndices_);
+    const auto                    lookupResult = lookupCopied(charges_, masses_, lookupIndices_);
     EXPECT_EQ(lookupResult[0], 1);
     EXPECT_EQ(lookupResult[1], 1);
 }
@@ -126,7 +118,7 @@ TEST_F(DensityFittingAmplitudeLookupTest, CanMoveConstruct)
 {
     DensityFittingAmplitudeLookup lookup(DensityFittingAmplitudeMethod::Unity);
     DensityFittingAmplitudeLookup lookupCopied(std::move(lookup));
-    const auto                    lookupResult = lookupCopied(atoms_, lookupIndices_);
+    const auto                    lookupResult = lookupCopied(charges_, masses_, lookupIndices_);
     EXPECT_EQ(lookupResult[0], 1);
     EXPECT_EQ(lookupResult[1], 1);
 }
index 6c4381747aa9b5c4dda5caa4883eb25c333db04f..2140b6226324a98f831890b8a8bbf7138e30555a 100644 (file)
@@ -59,7 +59,6 @@
 #include "gromacs/mdtypes/imdmodule.h"
 #include "gromacs/mdtypes/imdoutputprovider.h"
 #include "gromacs/mdtypes/imdpoptionprovider.h"
-#include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/options/basicoptions.h"
 #include "gromacs/options/ioptionscontainerwithsections.h"
 #include "gromacs/options/optionsection.h"
@@ -303,13 +302,13 @@ void ElectricField::calculateForces(const ForceProviderInput& forceProviderInput
 {
     if (isActive())
     {
-        const t_mdatoms& mdatoms = forceProviderInput.mdatoms_;
-        const double     t       = forceProviderInput.t_;
-        const t_commrec& cr      = forceProviderInput.cr_;
+        const double     t  = forceProviderInput.t_;
+        const t_commrec& cr = forceProviderInput.cr_;
 
         // NOTE: The non-conservative electric field does not have a virial
         ArrayRef<RVec> f = forceProviderOutput->forceWithVirial_.force_;
 
+        auto chargeA = forceProviderInput.chargeA_;
         for (int m = 0; (m < DIM); m++)
         {
             const real fieldStrength = gmx::c_fieldfac * field(m, t);
@@ -317,10 +316,10 @@ void ElectricField::calculateForces(const ForceProviderInput& forceProviderInput
             if (fieldStrength != 0)
             {
                 // TODO: Check parallellism
-                for (int i = 0; i < mdatoms.homenr; ++i)
+                for (int i = 0; i < forceProviderInput.homenr_; ++i)
                 {
                     // NOTE: Not correct with perturbed charges
-                    f[i][m] += mdatoms.chargeA[i] * fieldStrength;
+                    f[i][m] += chargeA[i] * fieldStrength;
                 }
             }
         }
index 8e19955ae37654f42471b1793933e171251ec962..fc514edae3e9d995fe3f6b9f82740b3f87013bb0 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2015,2016,2017,2018,2019,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.
@@ -106,13 +106,10 @@ public:
         }
 
         // Prepare a ForceProviderInput
-        t_mdatoms         md;
-        std::vector<real> chargeA{ 1 };
-        md.homenr  = ssize(chargeA);
-        md.chargeA = chargeA.data();
+        std::vector<real>  chargeA{ 1 };
         t_commrec          cr;
         matrix             boxDummy = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } };
-        ForceProviderInput forceProviderInput({}, md, 0.0, boxDummy, cr);
+        ForceProviderInput forceProviderInput({}, ssize(chargeA), chargeA, {}, 0.0, boxDummy, cr);
 
         // Prepare a ForceProviderOutput
         PaddedVector<RVec>  f = { { 0, 0, 0 } };
index 5208ffc5c1c5c40e877551e139f97d7d3367d622..99dc4758c4070546394ece82f252eecca21dc6c4 100644 (file)
@@ -641,7 +641,14 @@ static void computeSpecialForces(FILE*                          fplog,
      */
     if (stepWork.computeForces)
     {
-        gmx::ForceProviderInput  forceProviderInput(x, *mdatoms, t, box, *cr);
+        gmx::ForceProviderInput forceProviderInput(
+                x,
+                mdatoms->homenr,
+                gmx::arrayRefFromArray(mdatoms->chargeA, mdatoms->homenr),
+                gmx::arrayRefFromArray(mdatoms->massT, mdatoms->homenr),
+                t,
+                box,
+                *cr);
         gmx::ForceProviderOutput forceProviderOutput(forceWithVirialMtsLevel0, enerd);
 
         /* Collect forces from modules */
index 748241617bea4daa88b6bc8eb9b764393e9b89e8..3006bf88e971914fe5786a54ed72b97cde1b2a5f 100644 (file)
@@ -55,7 +55,6 @@
 struct gmx_enerdata_t;
 struct t_commrec;
 struct t_forcerec;
-struct t_mdatoms;
 
 namespace gmx
 {
@@ -82,28 +81,36 @@ class ForceProviderInput
 public:
     /*! \brief Constructor assembles all necessary force provider input data
      *
-     * \param[in]  x        Atomic positions
-     * \param[in]  cr       Communication record structure
-     * \param[in]  box      The simulation box
-     * \param[in]  time     The current time in the simulation
-     * \param[in]  mdatoms  The atomic data
+     * \param[in]  x        Atomic positions.
+     * \param[in]  homenr   Number of atoms on the domain.
+     * \param[in]  chargeA  Atomic charges for atoms on the domain.
+     * \param[in]  massT    Atomic masses for atoms on the domain.
+     * \param[in]  time     The current time in the simulation.
+     * \param[in]  box      The simulation box.
+     * \param[in]  cr       Communication record structure.
      */
     ForceProviderInput(ArrayRef<const RVec> x,
-                       const t_mdatoms&     mdatoms,
+                       int                  homenr,
+                       ArrayRef<const real> chargeA,
+                       ArrayRef<const real> massT,
                        double               time,
                        const matrix         box,
                        const t_commrec&     cr) :
         x_(x),
-        mdatoms_(mdatoms),
+        homenr_(homenr),
+        chargeA_(chargeA),
+        massT_(massT),
         t_(time),
         cr_(cr)
     {
         copy_mat(box, box_);
     }
 
-    ArrayRef<const RVec> x_;       //!< The atomic positions
-    const t_mdatoms&     mdatoms_; //!< Atomic data
-    double               t_;       //!< The current time in the simulation
+    ArrayRef<const RVec> x_; //!< The atomic positions
+    int                  homenr_;
+    ArrayRef<const real> chargeA_;
+    ArrayRef<const real> massT_;
+    double               t_; //!< The current time in the simulation
     matrix               box_ = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } }; //!< The simulation box
     const t_commrec&     cr_; //!< Communication record structure
 };
index b1dcb03d09c0f3eef63b75d5b6e2062fc581d262..5ea9cdd07c91b856f427acd5daecb42e9d2ac7ce 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,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.
@@ -68,8 +68,8 @@ void RestraintForceProvider::calculateForces(const ForceProviderInput& forceProv
 {
     GMX_ASSERT(restraint_, "Restraint must be initialized.");
 
-    const auto& mdatoms = forceProviderInput.mdatoms_;
-    GMX_ASSERT(mdatoms.homenr >= 0, "number of home atoms must be non-negative.");
+    const int homenr = forceProviderInput.homenr_;
+    GMX_ASSERT(homenr >= 0, "number of home atoms must be non-negative.");
 
     const auto& box = forceProviderInput.box_;
     GMX_ASSERT(check_box(PbcType::Unset, box) == nullptr, "Invalid box.");
@@ -80,7 +80,7 @@ void RestraintForceProvider::calculateForces(const ForceProviderInput& forceProv
     const auto& cr = forceProviderInput.cr_;
     const auto& t  = forceProviderInput.t_;
     // Cooperatively get Cartesian coordinates for center of mass of each site
-    RVec r1 = sites_[0].centerOfMass(cr, static_cast<size_t>(mdatoms.homenr), x, t);
+    RVec r1 = sites_[0].centerOfMass(cr, static_cast<size_t>(homenr), x, t);
     // r2 is to be constructed as
     // r2 = (site[N] - site[N-1]) + (site_{N-1} - site_{N-2}) + ... + (site_2 - site_1) + site_1
     // where the minimum image convention is applied to each path but not to the overall sum.
@@ -96,8 +96,8 @@ void RestraintForceProvider::calculateForces(const ForceProviderInput& forceProv
     // a big molecule in a small box.
     for (size_t i = 0; i < sites_.size() - 1; ++i)
     {
-        RVec a = sites_[i].centerOfMass(cr, static_cast<size_t>(mdatoms.homenr), x, t);
-        RVec b = sites_[i + 1].centerOfMass(cr, static_cast<size_t>(mdatoms.homenr), x, t);
+        RVec a = sites_[i].centerOfMass(cr, static_cast<size_t>(homenr), x, t);
+        RVec b = sites_[i + 1].centerOfMass(cr, static_cast<size_t>(homenr), x, t);
         // dr = minimum_image_vector(b - a)
         pbc_dx(&pbc, b, a, dr);
         r2[0] += dr[0];