From f7297c6b9de56884819379cfb0b3057762a9a1fc Mon Sep 17 00:00:00 2001 From: julianhille Date: Mon, 23 Feb 2026 12:28:43 +0100 Subject: [PATCH 1/4] added local force capping from #076, Ref. #072 --- mrmd/action/LennardJones.cpp | 61 +++--------------------- mrmd/action/LennardJones.hpp | 92 ++++++++++++++++++++++++++++++++++-- 2 files changed, 96 insertions(+), 57 deletions(-) diff --git a/mrmd/action/LennardJones.cpp b/mrmd/action/LennardJones.cpp index 8f49e5ff..db9f11d7 100644 --- a/mrmd/action/LennardJones.cpp +++ b/mrmd/action/LennardJones.cpp @@ -22,60 +22,13 @@ namespace mrmd::action { void LennardJones::apply(data::Atoms& atoms, HalfVerletList& verletList) { - energyAndVirial_ = data::EnergyAndVirialReducer(); - - pos_ = atoms.getPos(); - force_ = atoms.getForce(); - type_ = atoms.getType(); - verletList_ = verletList; - - auto policy = Kokkos::RangePolicy<>(0, atoms.numLocalAtoms); - Kokkos::parallel_reduce("LennardJones::applyForces", policy, *this, energyAndVirial_); - Kokkos::fence(); -} - -KOKKOS_FUNCTION -void LennardJones::operator()(const idx_t& idx, data::EnergyAndVirialReducer& energyAndVirial) const -{ - real_t posTmp[3]; - posTmp[0] = pos_(idx, 0); - posTmp[1] = pos_(idx, 1); - posTmp[2] = pos_(idx, 2); - - real_t forceTmp[3] = {0_r, 0_r, 0_r}; - - const auto numNeighbors = idx_c(HalfNeighborList::numNeighbor(verletList_, idx)); - for (idx_t n = 0; n < numNeighbors; ++n) - { - idx_t jdx = idx_c(HalfNeighborList::getNeighbor(verletList_, idx, n)); - assert(0 <= jdx); - - auto dx = posTmp[0] - pos_(jdx, 0); - auto dy = posTmp[1] - pos_(jdx, 1); - auto dz = posTmp[2] - pos_(jdx, 2); - - auto distSqr = dx * dx + dy * dy + dz * dz; - - if (distSqr > rcSqr_) continue; - - auto typeIdx = type_(idx) * numTypes_ + type_(jdx); - auto forceAndEnergy = LJ_.computeForceAndEnergy(distSqr, typeIdx); - assert(!std::isnan(forceAndEnergy.forceFactor)); - energyAndVirial.energy += forceAndEnergy.energy; - energyAndVirial.virial -= 0.5_r * forceAndEnergy.forceFactor * distSqr; - - forceTmp[0] += dx * forceAndEnergy.forceFactor; - forceTmp[1] += dy * forceAndEnergy.forceFactor; - forceTmp[2] += dz * forceAndEnergy.forceFactor; - - force_(jdx, 0) -= dx * forceAndEnergy.forceFactor; - force_(jdx, 1) -= dy * forceAndEnergy.forceFactor; - force_(jdx, 2) -= dz * forceAndEnergy.forceFactor; - } - - force_(idx, 0) += forceTmp[0]; - force_(idx, 1) += forceTmp[1]; - force_(idx, 2) += forceTmp[2]; + apply_if( + atoms, + verletList, + KOKKOS_LAMBDA( + const real_t, const real_t, const real_t, const real_t, const real_t, const real_t) { + return true; + }); } real_t LennardJones::getEnergy() const { return energyAndVirial_.energy; } diff --git a/mrmd/action/LennardJones.hpp b/mrmd/action/LennardJones.hpp index 4ab5cd94..5b6c0e43 100644 --- a/mrmd/action/LennardJones.hpp +++ b/mrmd/action/LennardJones.hpp @@ -85,14 +85,21 @@ class LennardJones data::EnergyAndVirialReducer energyAndVirial_; public: - KOKKOS_FUNCTION - void operator()(const idx_t& idx, data::EnergyAndVirialReducer& energyAndVirial) const; - real_t getEnergy() const; real_t getVirial() const; void apply(data::Atoms& atoms, HalfVerletList& verletList); + template BinaryPred> + void apply_if(const data::Atoms& atoms, + const HalfVerletList& verletList, + const BinaryPred& pred); + LennardJones(const real_t rc, const real_t& sigma, const real_t& epsilon, @@ -105,4 +112,83 @@ class LennardJones const idx_t& numTypes, const bool isShifted); }; + +template BinaryPred> +void LennardJones::apply_if(const data::Atoms& atoms, + const HalfVerletList& verletList, + const BinaryPred& pred) +{ + energyAndVirial_ = data::EnergyAndVirialReducer(); + pos_ = atoms.getPos(); + force_ = atoms.getForce(); + type_ = atoms.getType(); + verletList_ = verletList; + + auto policy = Kokkos::RangePolicy<>(0, atoms.numLocalAtoms); + + // avoid capturing this pointer + auto pos = pos_; + auto force = force_; + auto type = type_; + auto verletListLocal = verletList_; + auto rcSqr = rcSqr_; + auto LJ = LJ_; + auto numTypes = numTypes_; + auto predLocal = pred; + + auto kernel = KOKKOS_LAMBDA(const idx_t idx, data::EnergyAndVirialReducer& energyAndVirial) + { + real_t posTmp[3]; + posTmp[0] = pos(idx, 0); + posTmp[1] = pos(idx, 1); + posTmp[2] = pos(idx, 2); + + real_t forceTmp[3] = {0_r, 0_r, 0_r}; + + const auto numNeighbors = idx_c(HalfNeighborList::numNeighbor(verletListLocal, idx)); + for (idx_t n = 0; n < numNeighbors; ++n) + { + idx_t jdx = idx_c(HalfNeighborList::getNeighbor(verletListLocal, idx, n)); + assert(0 <= jdx); + + if (!predLocal(posTmp[0], posTmp[1], posTmp[2], pos(jdx, 0), pos(jdx, 1), pos(jdx, 2))) + continue; + + auto dx = posTmp[0] - pos(jdx, 0); + auto dy = posTmp[1] - pos(jdx, 1); + auto dz = posTmp[2] - pos(jdx, 2); + + auto distSqr = dx * dx + dy * dy + dz * dz; + + if (distSqr > rcSqr) continue; + + auto typeIdx = type(idx) * numTypes + type(jdx); + auto forceAndEnergy = LJ.computeForceAndEnergy(distSqr, typeIdx); + assert(!std::isnan(forceAndEnergy.forceFactor)); + energyAndVirial.energy += forceAndEnergy.energy; + energyAndVirial.virial -= 0.5_r * forceAndEnergy.forceFactor * distSqr; + + forceTmp[0] += dx * forceAndEnergy.forceFactor; + forceTmp[1] += dy * forceAndEnergy.forceFactor; + forceTmp[2] += dz * forceAndEnergy.forceFactor; + + force(jdx, 0) -= dx * forceAndEnergy.forceFactor; + force(jdx, 1) -= dy * forceAndEnergy.forceFactor; + force(jdx, 2) -= dz * forceAndEnergy.forceFactor; + } + + force(idx, 0) += forceTmp[0]; + force(idx, 1) += forceTmp[1]; + force(idx, 2) += forceTmp[2]; + }; + + Kokkos::parallel_reduce("LennardJones::apply_if", policy, kernel, energyAndVirial_); + Kokkos::fence(); +} + } // namespace mrmd::action \ No newline at end of file From 2fbbabb2b7ca409fd708afcd12e48a9d2ddb71c1 Mon Sep 17 00:00:00 2001 From: julianhille Date: Tue, 24 Feb 2026 15:59:19 +0100 Subject: [PATCH 2/4] generalized predicates to concepts in datatypes, Ref. #072 --- mrmd/action/LangevinThermostat.hpp | 6 ++---- mrmd/action/LennardJones.hpp | 16 ++-------------- mrmd/action/ThermodynamicForce.hpp | 7 +++---- mrmd/data/MultiHistogram.hpp | 8 +++----- mrmd/datatypes.hpp | 28 ++++++++++++++++++++++++++++ 5 files changed, 38 insertions(+), 27 deletions(-) diff --git a/mrmd/action/LangevinThermostat.hpp b/mrmd/action/LangevinThermostat.hpp index 9cb81cba..fcd86c84 100644 --- a/mrmd/action/LangevinThermostat.hpp +++ b/mrmd/action/LangevinThermostat.hpp @@ -38,8 +38,7 @@ class LangevinThermostat void apply(data::Atoms& atoms); - template UnaryPred> - void apply_if(data::Atoms& atoms, const UnaryPred& pred); + void apply_if(data::Atoms& atoms, const OnePositionPredicate auto& pred); void set(const real_t gamma, const real_t temperature, const real_t timestep) { @@ -53,8 +52,7 @@ class LangevinThermostat } }; -template UnaryPred> -void LangevinThermostat::apply_if(data::Atoms& atoms, const UnaryPred& pred) +void LangevinThermostat::apply_if(data::Atoms& atoms, const OnePositionPredicate auto& pred) { auto RNG = randPool_; auto pos = atoms.getPos(); diff --git a/mrmd/action/LennardJones.hpp b/mrmd/action/LennardJones.hpp index 5b6c0e43..8d487a97 100644 --- a/mrmd/action/LennardJones.hpp +++ b/mrmd/action/LennardJones.hpp @@ -90,15 +90,9 @@ class LennardJones void apply(data::Atoms& atoms, HalfVerletList& verletList); - template BinaryPred> void apply_if(const data::Atoms& atoms, const HalfVerletList& verletList, - const BinaryPred& pred); + const TwoPositionsPredicate auto& pred); LennardJones(const real_t rc, const real_t& sigma, @@ -113,15 +107,9 @@ class LennardJones const bool isShifted); }; -template BinaryPred> void LennardJones::apply_if(const data::Atoms& atoms, const HalfVerletList& verletList, - const BinaryPred& pred) + const TwoPositionsPredicate auto& pred) { energyAndVirial_ = data::EnergyAndVirialReducer(); pos_ = atoms.getPos(); diff --git a/mrmd/action/ThermodynamicForce.hpp b/mrmd/action/ThermodynamicForce.hpp index fc6e7d55..3613b71a 100644 --- a/mrmd/action/ThermodynamicForce.hpp +++ b/mrmd/action/ThermodynamicForce.hpp @@ -64,8 +64,7 @@ class ThermodynamicForce void update(const real_t& smoothingSigma, const real_t& smoothingIntensity); void apply(const data::Atoms& atoms) const; - template UnaryPred> - void apply_if(const data::Atoms& atoms, const UnaryPred& pred) const; + void apply_if(const data::Atoms& atoms, const OnePositionPredicate auto& pred) const; std::vector getMuLeft() const; std::vector getMuRight() const; @@ -85,8 +84,8 @@ class ThermodynamicForce const bool usePeriodicity = false); }; -template UnaryPred> -void ThermodynamicForce::apply_if(const data::Atoms& atoms, const UnaryPred& pred) const +void ThermodynamicForce::apply_if(const data::Atoms& atoms, + const OnePositionPredicate auto& pred) const { auto atomsPos = atoms.getPos(); auto atomsForce = atoms.getForce(); diff --git a/mrmd/data/MultiHistogram.hpp b/mrmd/data/MultiHistogram.hpp index 0c2ced26..5648f362 100644 --- a/mrmd/data/MultiHistogram.hpp +++ b/mrmd/data/MultiHistogram.hpp @@ -169,9 +169,6 @@ void transform(const MultiHistogram& input1, /** * Replaces histogram values with a new value if the bin position satisfies a predicate. * - * @tparam UnaryPred A unary predicate type that takes a real_t value and returns a boolean. - * Must satisfy std::predicate concept. - * * @param hist The MultiHistogram object whose values will be conditionally replaced. * @param pred A unary predicate function that is evaluated for each bin position. * If it returns true for a bin position, all histogram values at that bin @@ -179,8 +176,9 @@ void transform(const MultiHistogram& input1, * @param newValue The value to assign to histogram entries whose bin position satisfies * the predicate. */ -template UnaryPred> -void replace_if_bin_position(MultiHistogram& hist, const UnaryPred& pred, real_t newValue) +void replace_if_bin_position(MultiHistogram& hist, + const OneCoordinatePredicate auto& pred, + real_t newValue) { auto policy = Kokkos::MDRangePolicy>({0, 0}, {hist.numBins, hist.numHistograms}); diff --git a/mrmd/datatypes.hpp b/mrmd/datatypes.hpp index efbe38de..ebe04f5c 100644 --- a/mrmd/datatypes.hpp +++ b/mrmd/datatypes.hpp @@ -172,4 +172,32 @@ using NeighborList [[deprecated]] = Cabana::NeighborList; using HalfNeighborList = Cabana::NeighborList; using FullNeighborList = Cabana::NeighborList; +// Concept for predicates that take one coordinate as input, e.g. for spatially selective updating +// of histogram bins +template +concept OneCoordinatePredicate = std::predicate; + +// Concept for predicates that take one position as input, e.g. for spatially selective application +// of forces or thermostats +template +concept OnePositionPredicate = std::predicate; + +// Concept for predicates that take two positions as input, e.g. for conditions based on relative +// positions of two atoms +template +concept TwoPositionsPredicate = std::predicate; + } // namespace mrmd \ No newline at end of file From c657c7af42529a1c0ba4f6c1de0a2097267cf277 Mon Sep 17 00:00:00 2001 From: julianhille Date: Tue, 24 Feb 2026 16:53:59 +0100 Subject: [PATCH 3/4] made concepts device safe, Ref. #072 --- mrmd/action/LangevinThermostat.hpp | 8 ++++++-- mrmd/action/LennardJones.hpp | 10 ++++++---- mrmd/action/ThermodynamicForce.hpp | 9 ++++++--- mrmd/data/MultiHistogram.hpp | 6 +++--- 4 files changed, 21 insertions(+), 12 deletions(-) diff --git a/mrmd/action/LangevinThermostat.hpp b/mrmd/action/LangevinThermostat.hpp index fcd86c84..46590e20 100644 --- a/mrmd/action/LangevinThermostat.hpp +++ b/mrmd/action/LangevinThermostat.hpp @@ -38,7 +38,9 @@ class LangevinThermostat void apply(data::Atoms& atoms); - void apply_if(data::Atoms& atoms, const OnePositionPredicate auto& pred); + template + void apply_if(data::Atoms& atoms, const Pred& pred) + requires OnePositionPredicate; void set(const real_t gamma, const real_t temperature, const real_t timestep) { @@ -52,7 +54,9 @@ class LangevinThermostat } }; -void LangevinThermostat::apply_if(data::Atoms& atoms, const OnePositionPredicate auto& pred) +template +void LangevinThermostat::apply_if(data::Atoms& atoms, const Pred& pred) + requires OnePositionPredicate { auto RNG = randPool_; auto pos = atoms.getPos(); diff --git a/mrmd/action/LennardJones.hpp b/mrmd/action/LennardJones.hpp index 8d487a97..b54349ab 100644 --- a/mrmd/action/LennardJones.hpp +++ b/mrmd/action/LennardJones.hpp @@ -90,9 +90,9 @@ class LennardJones void apply(data::Atoms& atoms, HalfVerletList& verletList); - void apply_if(const data::Atoms& atoms, - const HalfVerletList& verletList, - const TwoPositionsPredicate auto& pred); + template + void apply_if(const data::Atoms& atoms, const HalfVerletList& verletList, const Pred& pred) + requires TwoPositionsPredicate; LennardJones(const real_t rc, const real_t& sigma, @@ -107,9 +107,11 @@ class LennardJones const bool isShifted); }; +template void LennardJones::apply_if(const data::Atoms& atoms, const HalfVerletList& verletList, - const TwoPositionsPredicate auto& pred) + const Pred& pred) + requires TwoPositionsPredicate { energyAndVirial_ = data::EnergyAndVirialReducer(); pos_ = atoms.getPos(); diff --git a/mrmd/action/ThermodynamicForce.hpp b/mrmd/action/ThermodynamicForce.hpp index 3613b71a..61b998f3 100644 --- a/mrmd/action/ThermodynamicForce.hpp +++ b/mrmd/action/ThermodynamicForce.hpp @@ -64,7 +64,9 @@ class ThermodynamicForce void update(const real_t& smoothingSigma, const real_t& smoothingIntensity); void apply(const data::Atoms& atoms) const; - void apply_if(const data::Atoms& atoms, const OnePositionPredicate auto& pred) const; + template + void apply_if(const data::Atoms& atoms, const Pred& pred) const + requires OnePositionPredicate; std::vector getMuLeft() const; std::vector getMuRight() const; @@ -84,8 +86,9 @@ class ThermodynamicForce const bool usePeriodicity = false); }; -void ThermodynamicForce::apply_if(const data::Atoms& atoms, - const OnePositionPredicate auto& pred) const +template +void ThermodynamicForce::apply_if(const data::Atoms& atoms, const Pred& pred) const + requires OnePositionPredicate { auto atomsPos = atoms.getPos(); auto atomsForce = atoms.getForce(); diff --git a/mrmd/data/MultiHistogram.hpp b/mrmd/data/MultiHistogram.hpp index 5648f362..aad68cf0 100644 --- a/mrmd/data/MultiHistogram.hpp +++ b/mrmd/data/MultiHistogram.hpp @@ -176,9 +176,9 @@ void transform(const MultiHistogram& input1, * @param newValue The value to assign to histogram entries whose bin position satisfies * the predicate. */ -void replace_if_bin_position(MultiHistogram& hist, - const OneCoordinatePredicate auto& pred, - real_t newValue) +template +void replace_if_bin_position(MultiHistogram& hist, const Pred& pred, real_t newValue) + requires OneCoordinatePredicate { auto policy = Kokkos::MDRangePolicy>({0, 0}, {hist.numBins, hist.numHistograms}); From 094ef7c816d1222a7aa282b2319b223866c8f194 Mon Sep 17 00:00:00 2001 From: Sebastian Eibl Date: Wed, 25 Feb 2026 17:45:25 +0100 Subject: [PATCH 4/4] Apply suggestions from code review --- mrmd/action/LangevinThermostat.hpp | 8 +++----- mrmd/action/LennardJones.hpp | 8 +++----- mrmd/action/ThermodynamicForce.hpp | 8 +++----- mrmd/data/MultiHistogram.hpp | 3 +-- 4 files changed, 10 insertions(+), 17 deletions(-) diff --git a/mrmd/action/LangevinThermostat.hpp b/mrmd/action/LangevinThermostat.hpp index 46590e20..1666cfaf 100644 --- a/mrmd/action/LangevinThermostat.hpp +++ b/mrmd/action/LangevinThermostat.hpp @@ -38,9 +38,8 @@ class LangevinThermostat void apply(data::Atoms& atoms); - template - void apply_if(data::Atoms& atoms, const Pred& pred) - requires OnePositionPredicate; + template + void apply_if(data::Atoms& atoms, const Pred& pred); void set(const real_t gamma, const real_t temperature, const real_t timestep) { @@ -54,9 +53,8 @@ class LangevinThermostat } }; -template +template void LangevinThermostat::apply_if(data::Atoms& atoms, const Pred& pred) - requires OnePositionPredicate { auto RNG = randPool_; auto pos = atoms.getPos(); diff --git a/mrmd/action/LennardJones.hpp b/mrmd/action/LennardJones.hpp index b54349ab..e3919c20 100644 --- a/mrmd/action/LennardJones.hpp +++ b/mrmd/action/LennardJones.hpp @@ -90,9 +90,8 @@ class LennardJones void apply(data::Atoms& atoms, HalfVerletList& verletList); - template - void apply_if(const data::Atoms& atoms, const HalfVerletList& verletList, const Pred& pred) - requires TwoPositionsPredicate; + template + void apply_if(const data::Atoms& atoms, const HalfVerletList& verletList, const Pred& pred); LennardJones(const real_t rc, const real_t& sigma, @@ -107,11 +106,10 @@ class LennardJones const bool isShifted); }; -template +template void LennardJones::apply_if(const data::Atoms& atoms, const HalfVerletList& verletList, const Pred& pred) - requires TwoPositionsPredicate { energyAndVirial_ = data::EnergyAndVirialReducer(); pos_ = atoms.getPos(); diff --git a/mrmd/action/ThermodynamicForce.hpp b/mrmd/action/ThermodynamicForce.hpp index 61b998f3..7ec46084 100644 --- a/mrmd/action/ThermodynamicForce.hpp +++ b/mrmd/action/ThermodynamicForce.hpp @@ -64,9 +64,8 @@ class ThermodynamicForce void update(const real_t& smoothingSigma, const real_t& smoothingIntensity); void apply(const data::Atoms& atoms) const; - template - void apply_if(const data::Atoms& atoms, const Pred& pred) const - requires OnePositionPredicate; + template + void apply_if(const data::Atoms& atoms, const Pred& pred) const; std::vector getMuLeft() const; std::vector getMuRight() const; @@ -86,9 +85,8 @@ class ThermodynamicForce const bool usePeriodicity = false); }; -template +template void ThermodynamicForce::apply_if(const data::Atoms& atoms, const Pred& pred) const - requires OnePositionPredicate { auto atomsPos = atoms.getPos(); auto atomsForce = atoms.getForce(); diff --git a/mrmd/data/MultiHistogram.hpp b/mrmd/data/MultiHistogram.hpp index aad68cf0..317d54cd 100644 --- a/mrmd/data/MultiHistogram.hpp +++ b/mrmd/data/MultiHistogram.hpp @@ -176,9 +176,8 @@ void transform(const MultiHistogram& input1, * @param newValue The value to assign to histogram entries whose bin position satisfies * the predicate. */ -template +template void replace_if_bin_position(MultiHistogram& hist, const Pred& pred, real_t newValue) - requires OneCoordinatePredicate { auto policy = Kokkos::MDRangePolicy>({0, 0}, {hist.numBins, hist.numHistograms});