Skip to content

Implement smooth interpolation of thermodynamic force #49

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
J-Hizzle wants to merge 10 commits into
base: main
Choose a base branch
from
Open

Implement smooth interpolation of thermodynamic force #49

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
660b81b
reverted interpolation and implemented overlapping bins, Ref. #047
J-Hizzle Aug 29, 2025
e8037ff
updated restoring thermodynamic force, Ref. #047
J-Hizzle Aug 29, 2025
bb5e16c
adapted restoring thermodynamic force, Ref. #047
J-Hizzle Aug 29, 2025
a5d974c
clang-tidy fix, Ref. #047
J-Hizzle Sep 1, 2025
de7a976
added constraint to application region, Ref. #047
J-Hizzle Sep 3, 2025
8a8dc95
bugfixes, Ref. #047
J-Hizzle Sep 3, 2025
5897b0c
clang-format fixes, Ref. #047
J-Hizzle Sep 3, 2025
aae3a90
adapted test to new application region boundaries, Ref. #047
J-Hizzle Sep 3, 2025
cbb1067
Update mrmd/action/ThermodynamicForce.cpp
XzzX Dec 5, 2025
570e030
Update mrmd/analysis/AxialDensityProfile.cpp
XzzX Dec 5, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions examples/BinaryLennardJones/SPARTIAN.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -209,6 +209,7 @@ void spartian(YAML::Node& config,
thermodynamicForce.getDensityProfile().min,
thermodynamicForce.getDensityProfile().max,
thermodynamicForce.getDensityProfile().numBins,
thermodynamicForce.getDensityProfile().binSize,
COORD_X);
densityProfile.scale(
1_r / (densityProfile.binSize * subdomain.diameter[1] * subdomain.diameter[2]));
Expand Down
1 change: 1 addition & 0 deletions examples/DensityFluctuations/DensityFluctuations.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -247,6 +247,7 @@ void LJ(Config& config)
subdomain.minCorner[0],
subdomain.maxCorner[0],
10,
config.densityBinWidth,
COORD_X);
densityProfile.scale(1_r / (densityProfile.binSize * config.Ly * config.Lz * 4));
auto fluctuation = analysis::getFluctuation(densityProfile, rho, 0);
Expand Down
134 changes: 115 additions & 19 deletions mrmd/action/ThermodynamicForce.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,55 +21,121 @@ namespace mrmd
{
namespace action
{
ThermodynamicForce::ThermodynamicForce(const std::vector<real_t>& targetDensity,

ThermodynamicForce::ThermodynamicForce(const std::vector<real_t>& targetDensities,
const data::Subdomain& subdomain,
const real_t& requestedDensityGridSpacing,
const real_t& requestedDensityBinWidth,
const std::vector<real_t>& thermodynamicForceModulation,
const std::vector<real_t>& thermodynamicForceModulations,
const bool enforceSymmetry,
const bool usePeriodicity)
: force_("thermodynamic-force",
subdomain.minCorner[0],
subdomain.maxCorner[0],
idx_c(std::ceil(subdomain.diameter[0] / requestedDensityBinWidth)),
idx_c(targetDensity.size())),
densityProfile_("density-profile", force_),
idx_c(std::ceil(subdomain.diameter[0] / requestedDensityGridSpacing)),
idx_c(targetDensities.size())),
densityProfile_("density-profile",
subdomain.minCorner[0],
subdomain.maxCorner[0],
idx_c(std::ceil(subdomain.diameter[0] / requestedDensityGridSpacing)),
idx_c(targetDensities.size())),
binVolume_(subdomain.diameter[1] * subdomain.diameter[2] * densityProfile_.binSize),
targetDensity_(targetDensity),
thermodynamicForceModulation_(thermodynamicForceModulation),
forceFactor_("force-factor", targetDensity.size()),
targetDensities_(targetDensities),
densityBinWidth_(requestedDensityBinWidth),
thermodynamicForceModulations_(thermodynamicForceModulations),
forceFactor_("force-factor", targetDensities.size()),
enforceSymmetry_(enforceSymmetry),
usePeriodicity_(usePeriodicity)
{
MRMD_HOST_CHECK_LESSEQUAL(
force_.binSize, requestedDensityBinWidth, "requested bin size is not achieved");
MRMD_HOST_CHECK_LESSEQUAL(force_.binSize,
requestedDensityGridSpacing,
"requested thermodynamic force bin size is not achieved");
MRMD_HOST_CHECK_LESSEQUAL(densityProfile_.binSize,
requestedDensityGridSpacing,
"requested density profile bin size is not achieved");

MRMD_HOST_CHECK_EQUAL(targetDensity.size(), thermodynamicForceModulation.size());
numTypes_ = idx_c(targetDensity.size());
MRMD_HOST_CHECK_EQUAL(targetDensities.size(), thermodynamicForceModulations.size());
numTypes_ = idx_c(targetDensities.size());
MRMD_HOST_CHECK_GREATER(numTypes_, 0);

auto hForceFactor = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), forceFactor_);
for (auto i = 0; i < numTypes_; ++i)
{
hForceFactor(i) = thermodynamicForceModulation_[i] / targetDensity_[i];
hForceFactor(i) = thermodynamicForceModulations_[i] / targetDensities_[i];
}
Kokkos::deep_copy(forceFactor_, hForceFactor);

}

ThermodynamicForce::ThermodynamicForce(const real_t targetDensity,
ThermodynamicForce::ThermodynamicForce(const std::vector<real_t>& targetDensities,
const data::Subdomain& subdomain,
const real_t& requestedDensityBinWidth,
const real_t thermodynamicForceModulation,
const real_t& requestedDensityGridSpacing,
const std::vector<real_t>& thermodynamicForceModulations,
const bool enforceSymmetry,
const bool usePeriodicity)
: ThermodynamicForce(targetDensities,
subdomain,
requestedDensityGridSpacing,
requestedDensityGridSpacing,
thermodynamicForceModulations,
enforceSymmetry,
usePeriodicity)
{
}

ThermodynamicForce::ThermodynamicForce(const real_t& targetDensity,
const data::Subdomain& subdomain,
const real_t& requestedDensityGridSpacing,
const real_t& thermodynamicForceModulation,
const bool enforceSymmetry,
const bool usePeriodicity)
: ThermodynamicForce(std::vector<real_t>{targetDensity},
: ThermodynamicForce({targetDensity},
subdomain,
requestedDensityBinWidth,
requestedDensityGridSpacing,
requestedDensityGridSpacing,
{thermodynamicForceModulation},
enforceSymmetry,
usePeriodicity)
{
}

ThermodynamicForce::ThermodynamicForce(const std::vector<real_t>& targetDensities,
const data::Subdomain& subdomain,
const idx_t& requestedDensityBinNumber,
const real_t& requestedDensityBinWidth,
const std::vector<real_t>& thermodynamicForceModulations,
const bool enforceSymmetry,
const bool usePeriodicity)
: force_("thermodynamic-force",
subdomain.minCorner[0],
subdomain.maxCorner[0],
requestedDensityBinNumber,
idx_c(targetDensities.size())),
densityProfile_("density-profile",
subdomain.minCorner[0],
subdomain.maxCorner[0],
requestedDensityBinNumber,
idx_c(targetDensities.size())),
binVolume_(subdomain.diameter[1] * subdomain.diameter[2] * densityProfile_.binSize),
targetDensities_(targetDensities),
densityBinWidth_(requestedDensityBinWidth),
thermodynamicForceModulations_(thermodynamicForceModulations),
forceFactor_("force-factor", targetDensities.size()),
enforceSymmetry_(enforceSymmetry),
usePeriodicity_(usePeriodicity)
{
MRMD_HOST_CHECK_EQUAL(targetDensities.size(), thermodynamicForceModulations.size());
numTypes_ = idx_c(targetDensities.size());
MRMD_HOST_CHECK_GREATER(numTypes_, 0);

auto hForceFactor = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), forceFactor_);
for (auto i = 0; i < numTypes_; ++i)
{
hForceFactor(i) = thermodynamicForceModulations_[i] / targetDensities_[i];
}
Kokkos::deep_copy(forceFactor_, hForceFactor);
}

void ThermodynamicForce::sample(data::Atoms& atoms)
{
densityProfile_ += analysis::getAxialDensityProfile(atoms.numLocalAtoms,
Expand All @@ -79,6 +145,7 @@ void ThermodynamicForce::sample(data::Atoms& atoms)
densityProfile_.min,
densityProfile_.max,
densityProfile_.numBins,
densityBinWidth_,
COORD_X);

++densityProfileSamples_;
Expand All @@ -101,7 +168,36 @@ void ThermodynamicForce::update(const real_t& smoothingSigma, const real_t& smoo
auto smoothedDensityGradient = data::gradient(smoothedDensityProfile, usePeriodicity_);
smoothedDensityGradient.scale(forceFactor_);

force_ -= smoothedDensityGradient;
force_ -= util::interpolate(smoothedDensityGradient, force_.createGrid_d());

// reset sampling data
Kokkos::deep_copy(densityProfile_.data, 0_r);
densityProfileSamples_ = 0;
}

void ThermodynamicForce::update(const real_t& smoothingSigma,
const real_t& smoothingIntensity,
const util::ApplicationRegion& applicationRegion)
{
MRMD_HOST_CHECK_GREATER(densityProfileSamples_, 0);

if (enforceSymmetry_)
{
densityProfile_.makeSymmetric();
}

auto normalizationFactor = 1_r / (binVolume_ * real_c(densityProfileSamples_));
densityProfile_.scale(normalizationFactor);

auto smoothedDensityProfile =
data::smoothen(densityProfile_, smoothingSigma, smoothingIntensity, usePeriodicity_);
auto smoothedDensityGradient = data::gradient(smoothedDensityProfile, usePeriodicity_);
smoothedDensityGradient.scale(forceFactor_);

auto constrainedDensityGradient =
util::constrainToApplicationRegion(smoothedDensityGradient, applicationRegion);

force_ -= util::interpolate(constrainedDensityGradient, force_.createGrid_d());

// reset sampling data
Kokkos::deep_copy(densityProfile_.data, 0_r);
Expand Down
33 changes: 27 additions & 6 deletions mrmd/action/ThermodynamicForce.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include "data/Subdomain.hpp"
#include "datatypes.hpp"
#include "util/ApplicationRegion.hpp"
#include "util/interpolation.hpp"
#include "weighting_function/Slab.hpp"

namespace mrmd
Expand All @@ -33,8 +34,9 @@ class ThermodynamicForce
data::MultiHistogram densityProfile_;
idx_t densityProfileSamples_ = 0;
real_t binVolume_;
const std::vector<real_t> targetDensity_;
const std::vector<real_t> thermodynamicForceModulation_;
const std::vector<real_t> targetDensities_;
const real_t densityBinWidth_;
const std::vector<real_t> thermodynamicForceModulations_;
idx_t numTypes_;

ScalarView forceFactor_; ///< precalculated prefactor for force calculation
Expand All @@ -61,23 +63,42 @@ class ThermodynamicForce

void sample(data::Atoms& atoms);
void update(const real_t& smoothingSigma, const real_t& smoothingIntensity);
void update(const real_t& smoothingSigma,
const real_t& smoothingIntensity,
const util::ApplicationRegion& applicationRegion);
void apply(const data::Atoms& atoms) const;
void apply(const data::Atoms& atoms, const weighting_function::Slab& slab) const;
void apply(const data::Atoms& atoms, const util::ApplicationRegion& applicationRegion) const;
std::vector<real_t> getMuLeft() const;
std::vector<real_t> getMuRight() const;

ThermodynamicForce(const std::vector<real_t>& targetDensity,
ThermodynamicForce(const std::vector<real_t>& targetDensities,
const data::Subdomain& subdomain,
const real_t& requestedDensityGridSpacing,
const real_t& requestedDensityBinWidth,
const std::vector<real_t>& thermodynamicForceModulation,
const std::vector<real_t>& thermodynamicForceModulations,
const bool enforceSymmetry = false,
const bool usePeriodicity = false);
Comment on lines +75 to 81
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing documentation to distinguish requestedDensityGridSpacing from requestedDensityBinWidth. The ThermodynamicForce constructors use both parameters but their distinct purposes are not documented. Based on the implementation, gridSpacing appears to define the spacing between bin centers in the force grid, while binWidth defines the width of bins used during density profile sampling. This distinction should be documented in comments for the constructor parameters to avoid confusion.

Copilot uses AI. Check for mistakes.
Copy link
Contributor Author

@J-Hizzle J-Hizzle Dec 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point, could be a naming inconsistency that slipped my attention.


ThermodynamicForce(const real_t targetDensity,
ThermodynamicForce(const std::vector<real_t>& targetDensities,
const data::Subdomain& subdomain,
const real_t& requestedDensityGridSpacing,
const std::vector<real_t>& thermodynamicForceModulations,
const bool enforceSymmetry = false,
const bool usePeriodicity = false);

ThermodynamicForce(const real_t& targetDensity,
const data::Subdomain& subdomain,
const real_t& requestedDensityGridSpacing,
const real_t& thermodynamicForceModulation,
const bool enforceSymmetry = false,
const bool usePeriodicity = false);

ThermodynamicForce(const std::vector<real_t>& targetDensities,
const data::Subdomain& subdomain,
const idx_t& requestedDensityBinNumber,
const real_t& requestedDensityBinWidth,
const real_t thermodynamicForceModulation,
const std::vector<real_t>& thermodynamicForceModulations,
const bool enforceSymmetry = false,
const bool usePeriodicity = false);
};
Expand Down
19 changes: 12 additions & 7 deletions mrmd/analysis/AxialDensityProfile.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,13 +20,15 @@ namespace mrmd
{
namespace analysis
{

data::MultiHistogram getAxialDensityProfile(const idx_t numAtoms,
const data::Atoms::pos_t& positions,
const data::Atoms::type_t& types,
const int64_t numTypes,
const real_t min,
const real_t max,
const int64_t numBins,
const real_t binWidth,
const int64_t axis)
{
MRMD_HOST_CHECK_GREATEREQUAL(max, min);
Expand All @@ -37,15 +39,18 @@ data::MultiHistogram getAxialDensityProfile(const idx_t numAtoms,
data::MultiHistogram histogram("density-profile", min, max, numBins, numTypes);
MultiScatterView scatter(histogram.data);

auto policy = Kokkos::RangePolicy<>(0, numAtoms);
auto kernel = KOKKOS_LAMBDA(const idx_t idx)
auto policy = Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {numAtoms, numBins});
auto kernel = KOKKOS_LAMBDA(const idx_t atomIdx, const idx_t binIdx)
{
MRMD_DEVICE_ASSERT_GREATEREQUAL(types(idx), 0);
MRMD_DEVICE_ASSERT_LESS(types(idx), numTypes);
auto bin = histogram.getBin(positions(idx, axis));
if (bin == -1) return;
MRMD_DEVICE_ASSERT_GREATEREQUAL(types(atomIdx), 0);
MRMD_DEVICE_ASSERT_LESS(types(atomIdx), numTypes);
auto binPos = histogram.getBinPosition(binIdx);
auto atomPos = positions(atomIdx, axis);
auto access = scatter.access();
access(bin, types(idx)) += 1_r;
if (atomPos >= binPos - binWidth / 2 && atomPos <= binPos + binWidth / 2)
{
access(binIdx, types(atomIdx)) += 1_r;
}
};
Comment on lines +42 to 54
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The new implementation has significantly worse computational complexity. The algorithm changed from O(N) to O(N*M) where N is the number of atoms and M is the number of bins. The old implementation (commented out below) iterated only over atoms and directly computed which bin each atom belonged to. The new implementation uses a 2D iteration over all (atom, bin) pairs and checks if each atom falls within each bin's range. For large systems, this could cause significant performance degradation. Consider optimizing by iterating only over atoms and identifying which bins each atom should contribute to based on the binWidth parameter.

Copilot uses AI. Check for mistakes.
Copy link
Contributor Author

@J-Hizzle J-Hizzle Dec 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a great point. We do indeed observe a significant loss in performance for this implementation. Should be changed accordingly.

Copy link
Owner

@XzzX XzzX Dec 9, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What was the problem with the initial implementation? I do not get what your changes do differently? Has it to do with the out-of-bounds access you mentioned?

Kokkos::parallel_for(policy, kernel);
Kokkos::Experimental::contribute(histogram.data, scatter);
Expand Down
1 change: 1 addition & 0 deletions mrmd/analysis/AxialDensityProfile.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ data::MultiHistogram getAxialDensityProfile(const idx_t numAtoms,
const real_t min,
const real_t max,
const int64_t numBins,
const real_t binWidth,
const int64_t axis);

} // namespace analysis
Expand Down
31 changes: 30 additions & 1 deletion mrmd/analysis/AxialDensityProfile.test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ TEST(AxialDensityProfile, histogram)
auto atoms = initAtoms();

auto histogram = getAxialDensityProfile(
atoms.numLocalAtoms, atoms.getPos(), atoms.getType(), 2, 0_r, 10_r, 10, COORD_X);
atoms.numLocalAtoms, atoms.getPos(), atoms.getType(), 2, 0_r, 10_r, 10, 1_r, COORD_X);
auto h_data = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogram.data);

for (auto i = 0; i < 10; ++i)
Expand All @@ -69,5 +69,34 @@ TEST(AxialDensityProfile, histogram)
EXPECT_FLOAT_EQ(h_data(i, 1), real_c(11 - (i + 1)));
}
}

TEST(AxialDensityProfile, overlappingBins)
{
auto atoms = initAtoms();

auto histogram = getAxialDensityProfile(
atoms.numLocalAtoms, atoms.getPos(), atoms.getType(), 2, 0_r, 10_r, 10, 3_r, COORD_X);
auto h_data = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogram.data);

for (auto i = 0; i < 10; ++i)
{
if (i == 0)
{
EXPECT_FLOAT_EQ(h_data(i, 0), (real_c(i) + real_c(i + 1) + real_c(i + 2)));
EXPECT_FLOAT_EQ(h_data(i, 1), (real_c(11 - (i + 1)) + real_c(11 - (i + 2))));
continue;
}
if (i == 9)
{
EXPECT_FLOAT_EQ(h_data(i, 0), (real_c(i) + real_c(i + 1)));
EXPECT_FLOAT_EQ(h_data(i, 1), (real_c(11 - i) + real_c(11 - (i + 1))));
continue;
}
EXPECT_FLOAT_EQ(h_data(i, 0), (real_c(i) + real_c(i + 1) + real_c(i + 2)));
EXPECT_FLOAT_EQ(h_data(i, 1),
(real_c(11 - i) + real_c(11 - (i + 1)) + real_c(11 - (i + 2))));
}
}

} // namespace analysis
} // namespace mrmd
18 changes: 18 additions & 0 deletions mrmd/data/MultiHistogram.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,24 @@ ScalarView::HostMirror MultiHistogram::createGrid() const
return grid;
}

ScalarView MultiHistogram::createGrid_d() const
{
MRMD_HOST_CHECK_GREATEREQUAL(numBins, 0);

ScalarView grid("grid", numBins);
const real_t min_ = min;
const real_t binSize_ = binSize;
auto policy = Kokkos::RangePolicy<>(0, numBins);
auto kernel = KOKKOS_LAMBDA(const idx_t idx)
Copy link
Owner

@XzzX XzzX Dec 9, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not use getBinPosition?

Copy link
Contributor Author

@J-Hizzle J-Hizzle Dec 17, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did it to avoid capturing the this pointer therefore I opted for this redundant but inline implementation. How can I use getBinPosition here without capturing the this pointer?

{
grid[idx] = min_ + (real_c(idx) + 0.5_r) * binSize_;
};
Kokkos::parallel_for("MultiHistogram::scale", policy, kernel);
Copy link
Owner

@XzzX XzzX Dec 9, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wrong label

Kokkos::fence();

return grid;
}

MultiHistogram& MultiHistogram::operator+=(const MultiHistogram& rhs)
{
transform(*this, rhs, *this, bin_op::Add());
Expand Down
Loading