Merge pull request #1955 from glotzerlab/fix-hpmc-potential-energy

Fix hpmc potential energy calculation
glotzerlab · Nov 26, 2024 · 4eedb97 · 4eedb97
2 parents 324ff00 + 1102418
commit 4eedb97
Show file tree

Hide file tree

Showing 3 changed files with 23 additions and 10 deletions.
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -26,6 +26,9 @@ Change Log
   (`#1944 <https://github.com/glotzerlab/hoomd-blue/pull/1944>`__).
 * Read after write hazard in the GPU implementation of ``hoomd.md.mesh.conservation.Volume``
   (`#1953 <https://github.com/glotzerlab/hoomd-blue/pull/1953>`__).
+* ``hoomd.hpmc.pair.Pair.energy`` now computes the correct energy when there are multiple pair
+  potentials with different ``r_cut`` values
+  (`#1955 <https://github.com/glotzerlab/hoomd-blue/pull/1955>`__).
 
 *Added*
 

diff --git a/hoomd/hpmc/IntegratorHPMCMono.h b/hoomd/hpmc/IntegratorHPMCMono.h
@@ -1145,17 +1145,20 @@ double IntegratorHPMCMono<Shape>::computePairEnergy(uint64_t timestep,
                             if (h_tag.data[i] <= h_tag.data[j])
                                 {
                                 LongReal r_squared = dot(r_ij, r_ij);
-                                if (selected_pair
-                                    && r_squared < selected_pair->getRCutSquaredTotal(typ_i, typ_j))
+                                if (selected_pair)
                                     {
-                                    energy += selected_pair->energy(r_squared,
-                                                                    r_ij,
-                                                                    typ_i,
-                                                                    orientation_i,
-                                                                    h_charge.data[i],
-                                                                    typ_j,
-                                                                    orientation_j,
-                                                                    h_charge.data[j]);
+                                    if (r_squared
+                                        < selected_pair->getRCutSquaredTotal(typ_i, typ_j))
+                                        {
+                                        energy += selected_pair->energy(r_squared,
+                                                                        r_ij,
+                                                                        typ_i,
+                                                                        orientation_i,
+                                                                        h_charge.data[i],
+                                                                        typ_j,
+                                                                        orientation_j,
+                                                                        h_charge.data[j]);
+                                        }
                                     }
                                 else
                                     {

diff --git a/hoomd/hpmc/pytest/test_pair_lennard_jones.py b/hoomd/hpmc/pytest/test_pair_lennard_jones.py
@@ -205,6 +205,13 @@ def test_multiple_pair_potentials(mc_simulation_factory):
     assert simulation.operations.integrator.pair_energy == pytest.approx(
         expected=-3.0, rel=1e-5)
 
+    # check that individual energies are computed correctly with varying r_cut
+    lennard_jones_2.params[('A', 'A')] = dict(epsilon=2.0, sigma=1.0, r_cut=0)
+    assert simulation.operations.integrator.pair_energy == pytest.approx(
+        expected=-1.0, rel=1e-5)
+    assert lennard_jones_1.energy == pytest.approx(expected=-1.0, rel=1e-5)
+    assert lennard_jones_2.energy == pytest.approx(expected=0.0, abs=1e-5)
+
 
 def test_logging():
     hoomd.conftest.logging_check(