From 52a2063fcaaea787e2258c5e8b48123ab1cb44c2 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Wed, 12 Jun 2024 14:19:36 -0400
Subject: [PATCH 01/32] Update CMake lists

---
 hoomd/hpmc/CMakeLists.txt        | 2 ++
 hoomd/hpmc/pair/CMakeLists.txt   | 1 +
 hoomd/hpmc/pytest/CMakeLists.txt | 1 +
 3 files changed, 4 insertions(+)

diff --git a/hoomd/hpmc/CMakeLists.txt b/hoomd/hpmc/CMakeLists.txt
index 32d556890a..9ec7caea7f 100644
--- a/hoomd/hpmc/CMakeLists.txt
+++ b/hoomd/hpmc/CMakeLists.txt
@@ -39,6 +39,7 @@ set(_hpmc_sources   module.cc
                     ExternalFieldWall.cc
                     PairPotential.cc
                     PairPotentialLennardJones.cc
+                    PairPotentialExpandedGaussian.cc
                     PairPotentialStep.cc
                     PairPotentialUnion.cc
                     PairPotentialAngularStep.cc
@@ -81,6 +82,7 @@ set(_hpmc_headers
     OBBTree.h
     PairPotential.h
     PairPotentialLennardJones.h
+    PairPotentialExpandedGaussian.h
     PairPotentialStep.h
     PairPotentialUnion.h
     PairPotentialAngularStep.h
diff --git a/hoomd/hpmc/pair/CMakeLists.txt b/hoomd/hpmc/pair/CMakeLists.txt
index 5d1601d8f6..4170d377f3 100644
--- a/hoomd/hpmc/pair/CMakeLists.txt
+++ b/hoomd/hpmc/pair/CMakeLists.txt
@@ -1,5 +1,6 @@
 set(files __init__.py
         lennard_jones.py
+        expanded_gaussian.py
         pair.py
         step.py
         union.py
diff --git a/hoomd/hpmc/pytest/CMakeLists.txt b/hoomd/hpmc/pytest/CMakeLists.txt
index e085c9603d..2bfff5d32c 100644
--- a/hoomd/hpmc/pytest/CMakeLists.txt
+++ b/hoomd/hpmc/pytest/CMakeLists.txt
@@ -14,6 +14,7 @@ set(files __init__.py
           test_shape_utils.py
           test_move_size_tuner.py
           test_pair_lennard_jones.py
+          test_pair_expanded_gaussian.py
           test_pair_step.py
           test_pair_user.py
           test_pair_union.py

From 866db90cb882a3f7447852f93f9b760200b0b11d Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Wed, 12 Jun 2024 14:22:34 -0400
Subject: [PATCH 02/32] Create PairPotentialExpandedGaussian.cc

---
 hoomd/hpmc/PairPotentialExpandedGaussian.cc | 96 +++++++++++++++++++++
 1 file changed, 96 insertions(+)
 create mode 100644 hoomd/hpmc/PairPotentialExpandedGaussian.cc

diff --git a/hoomd/hpmc/PairPotentialExpandedGaussian.cc b/hoomd/hpmc/PairPotentialExpandedGaussian.cc
new file mode 100644
index 0000000000..3cfaa585ad
--- /dev/null
+++ b/hoomd/hpmc/PairPotentialExpandedGaussian.cc
@@ -0,0 +1,96 @@
+// Copyright (c) 2009-2024 The Regents of the University of Michigan.
+// Part of HOOMD-blue, released under the BSD 3-Clause License.
+
+#include "PairPotentialLennardJones.h"
+
+namespace hoomd
+    {
+namespace hpmc
+    {
+
+PairPotentialLennardJones::PairPotentialLennardJones(std::shared_ptr<SystemDefinition> sysdef)
+    : PairPotential(sysdef), m_params(m_type_param_index.getNumElements())
+    {
+    }
+
+LongReal PairPotentialLennardJones::energy(const LongReal r_squared,
+                                           const vec3<LongReal>& r_ij,
+                                           const unsigned int type_i,
+                                           const quat<LongReal>& q_i,
+                                           const LongReal charge_i,
+                                           const unsigned int type_j,
+                                           const quat<LongReal>& q_j,
+                                           const LongReal charge_j) const
+    {
+    unsigned int param_index = m_type_param_index(type_i, type_j);
+    const auto& param = m_params[param_index];
+
+    LongReal lj2 = param.epsilon_x_4 * param.sigma_6;
+    LongReal lj1 = lj2 * param.sigma_6;
+
+    LongReal r_2_inverse = LongReal(1.0) / r_squared;
+    LongReal r_6_inverse = r_2_inverse * r_2_inverse * r_2_inverse;
+
+    LongReal energy = r_6_inverse * (lj1 * r_6_inverse - lj2);
+
+    if (m_mode == shift || (m_mode == xplor && param.r_on_squared >= param.r_cut_squared))
+        {
+        LongReal r_cut_2_inverse = LongReal(1.0) / param.r_cut_squared;
+        LongReal r_cut_6_inverse = r_cut_2_inverse * r_cut_2_inverse * r_cut_2_inverse;
+        energy -= r_cut_6_inverse * (lj1 * r_cut_6_inverse - lj2);
+        }
+
+    if (m_mode == xplor && r_squared > param.r_on_squared)
+        {
+        LongReal a = param.r_cut_squared - param.r_on_squared;
+        LongReal denominator = a * a * a;
+
+        LongReal b = param.r_cut_squared - r_squared;
+        LongReal numerator = b * b
+                             * (param.r_cut_squared + LongReal(2.0) * r_squared
+                                - LongReal(3.0) * param.r_on_squared);
+        energy *= numerator / denominator;
+        }
+
+    return energy;
+    }
+
+void PairPotentialLennardJones::setParamsPython(pybind11::tuple typ, pybind11::dict params)
+    {
+    auto pdata = m_sysdef->getParticleData();
+    auto type_i = pdata->getTypeByName(typ[0].cast<std::string>());
+    auto type_j = pdata->getTypeByName(typ[1].cast<std::string>());
+    unsigned int param_index_1 = m_type_param_index(type_i, type_j);
+    m_params[param_index_1] = ParamType(params);
+    unsigned int param_index_2 = m_type_param_index(type_j, type_i);
+    m_params[param_index_2] = ParamType(params);
+
+    notifyRCutChanged();
+    }
+
+pybind11::dict PairPotentialLennardJones::getParamsPython(pybind11::tuple typ)
+    {
+    auto pdata = m_sysdef->getParticleData();
+    auto type_i = pdata->getTypeByName(typ[0].cast<std::string>());
+    auto type_j = pdata->getTypeByName(typ[1].cast<std::string>());
+    unsigned int param_index = m_type_param_index(type_i, type_j);
+    return m_params[param_index].asDict();
+    }
+
+namespace detail
+    {
+void exportPairPotentialLennardJones(pybind11::module& m)
+    {
+    pybind11::class_<PairPotentialLennardJones,
+                     PairPotential,
+                     std::shared_ptr<PairPotentialLennardJones>>(m, "PairPotentialLennardJones")
+        .def(pybind11::init<std::shared_ptr<SystemDefinition>>())
+        .def("setParams", &PairPotentialLennardJones::setParamsPython)
+        .def("getParams", &PairPotentialLennardJones::getParamsPython)
+        .def_property("mode",
+                      &PairPotentialLennardJones::getMode,
+                      &PairPotentialLennardJones::setMode);
+    }
+    } // end namespace detail
+    } // end namespace hpmc
+    } // end namespace hoomd

From 13bcab18529ded3540498652ba72e185dcf3eeaf Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Wed, 12 Jun 2024 14:22:44 -0400
Subject: [PATCH 03/32] Create PairPotentialExpandedGaussian.h

---
 hoomd/hpmc/PairPotentialExpandedGaussian.h | 138 +++++++++++++++++++++
 1 file changed, 138 insertions(+)
 create mode 100644 hoomd/hpmc/PairPotentialExpandedGaussian.h

diff --git a/hoomd/hpmc/PairPotentialExpandedGaussian.h b/hoomd/hpmc/PairPotentialExpandedGaussian.h
new file mode 100644
index 0000000000..50c43ee556
--- /dev/null
+++ b/hoomd/hpmc/PairPotentialExpandedGaussian.h
@@ -0,0 +1,138 @@
+// Copyright (c) 2009-2024 The Regents of the University of Michigan.
+// Part of HOOMD-blue, released under the BSD 3-Clause License.
+
+#pragma once
+
+#include "PairPotential.h"
+
+namespace hoomd
+    {
+namespace hpmc
+    {
+
+/*** Compute Lennard-Jones energy between two particles.
+
+For use with HPMC simulations.
+*/
+class PairPotentialLennardJones : public hpmc::PairPotential
+    {
+    public:
+    PairPotentialLennardJones(std::shared_ptr<SystemDefinition> sysdef);
+    virtual ~PairPotentialLennardJones() { }
+
+    virtual LongReal energy(const LongReal r_squared,
+                            const vec3<LongReal>& r_ij,
+                            const unsigned int type_i,
+                            const quat<LongReal>& q_i,
+                            const LongReal charge_i,
+                            const unsigned int type_j,
+                            const quat<LongReal>& q_j,
+                            const LongReal charge_j) const;
+
+    /// Compute the non-additive cuttoff radius
+    virtual LongReal computeRCutNonAdditive(unsigned int type_i, unsigned int type_j) const
+        {
+        unsigned int param_index = m_type_param_index(type_i, type_j);
+        return slow::sqrt(m_params[param_index].r_cut_squared);
+        }
+
+    /// Set type pair dependent parameters to the potential.
+    void setParamsPython(pybind11::tuple typ, pybind11::dict params);
+
+    /// Get type pair dependent parameters.
+    pybind11::dict getParamsPython(pybind11::tuple typ);
+
+    void setMode(const std::string& mode_str)
+        {
+        if (mode_str == "none")
+            {
+            m_mode = no_shift;
+            }
+        else if (mode_str == "shift")
+            {
+            m_mode = shift;
+            }
+        else if (mode_str == "xplor")
+            {
+            m_mode = xplor;
+            }
+        else
+            {
+            throw std::domain_error("Invalid mode " + mode_str);
+            }
+        }
+
+    std::string getMode()
+        {
+        std::string result = "none";
+
+        if (m_mode == shift)
+            {
+            result = "shift";
+            }
+        if (m_mode == xplor)
+            {
+            result = "xplor";
+            }
+
+        return result;
+        }
+
+    protected:
+    /// Shifting modes that can be applied to the energy
+    enum EnergyShiftMode
+        {
+        no_shift = 0,
+        shift,
+        xplor
+        };
+
+    /// Type pair parameters of LJ potential
+    struct ParamType
+        {
+        ParamType()
+            {
+            sigma_6 = 0;
+            epsilon_x_4 = 0;
+            r_cut_squared = 0;
+            r_on_squared = 0;
+            }
+
+        ParamType(pybind11::dict v)
+            {
+            auto sigma(v["sigma"].cast<LongReal>());
+            auto epsilon(v["epsilon"].cast<LongReal>());
+            auto r_cut(v["r_cut"].cast<LongReal>());
+            auto r_on(v["r_on"].cast<LongReal>());
+
+            sigma_6 = sigma * sigma * sigma * sigma * sigma * sigma;
+            epsilon_x_4 = LongReal(4.0) * epsilon;
+            r_cut_squared = r_cut * r_cut;
+            r_on_squared = r_on * r_on;
+            }
+
+        pybind11::dict asDict()
+            {
+            pybind11::dict result;
+
+            result["sigma"] = pow(sigma_6, 1. / 6.);
+            result["epsilon"] = epsilon_x_4 / 4.0;
+            result["r_cut"] = slow::sqrt(r_cut_squared);
+            result["r_on"] = slow::sqrt(r_on_squared);
+
+            return result;
+            }
+
+        LongReal sigma_6;
+        LongReal epsilon_x_4;
+        LongReal r_cut_squared;
+        LongReal r_on_squared;
+        };
+
+    std::vector<ParamType> m_params;
+
+    EnergyShiftMode m_mode = no_shift;
+    };
+
+    } // end namespace hpmc
+    } // end namespace hoomd

From 8383b84fa257e529b8f69ccbca0e3402c2845b80 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Wed, 12 Jun 2024 14:22:57 -0400
Subject: [PATCH 04/32] Create expanded_gaussian.py

---
 hoomd/hpmc/pair/expanded_gaussian.py | 99 ++++++++++++++++++++++++++++
 1 file changed, 99 insertions(+)
 create mode 100644 hoomd/hpmc/pair/expanded_gaussian.py

diff --git a/hoomd/hpmc/pair/expanded_gaussian.py b/hoomd/hpmc/pair/expanded_gaussian.py
new file mode 100644
index 0000000000..ac680ab1e1
--- /dev/null
+++ b/hoomd/hpmc/pair/expanded_gaussian.py
@@ -0,0 +1,99 @@
+# Copyright (c) 2009-2024 The Regents of the University of Michigan.
+# Part of HOOMD-blue, released under the BSD 3-Clause License.
+
+"""Lennard-Jones pair potential.
+
+.. invisible-code-block: python
+
+    simulation = hoomd.util.make_example_simulation()
+    sphere = hoomd.hpmc.integrate.Sphere()
+    sphere.shape['A'] = dict(diameter=0.0)
+    simulation.operations.integrator = sphere
+"""
+
+import hoomd
+
+from .pair import Pair
+
+
+@hoomd.logging.modify_namespace(('hpmc', 'pair', 'LennardJones'))
+class LennardJones(Pair):
+    """Lennard-Jones pair potential (HPMC).
+
+    Args:
+        default_r_cut (float): Default cutoff radius :math:`[\\mathrm{length}]`.
+        default_r_on (float): Default XPLOR on radius
+          :math:`[\\mathrm{length}]`.
+        mode (str): Energy shifting/smoothing mode.
+
+    `LennardJones` computes the Lennard-Jones pair potential between every pair
+    of particles in the simulation state. The functional form of the potential,
+    including its behavior under shifting modes, is identical to that in
+    the MD pair potential `hoomd.md.pair.LJ`.
+
+    See Also:
+        `hoomd.md.pair.LJ`
+
+        `hoomd.md.pair`
+
+    .. rubric:: Example
+
+    .. code-block:: python
+
+        lennard_jones =  hoomd.hpmc.pair.LennardJones()
+        lennard_jones.params[('A', 'A')] = dict(epsilon=1, sigma=1, r_cut=2.5)
+        simulation.operations.integrator.pair_potentials = [lennard_jones]
+
+    .. py:attribute:: params
+
+        The potential parameters. The dictionary has the following keys:
+
+        * ``epsilon`` (`float`, **required**) -
+          Energy well depth :math:`\\varepsilon` :math:`[\\mathrm{energy}]`.
+        * ``sigma`` (`float`, **required**) -
+          Characteristic length scale :math:`\\sigma`
+          :math:`[\\mathrm{length}]`.
+        * ``r_cut`` (`float`): Cutoff radius :math:`[\\mathrm{length}]`.
+          Defaults to the value given in ``default_r_cut`` on construction.
+        * ``r_on`` (`float`): XPLOR on radius :math:`[\\mathrm{length}]`.
+          Defaults to the value given in ``default_r_on`` on construction.
+
+        Type: `TypeParameter` [`tuple` [``particle_type``, ``particle_type``],
+        `dict`]
+
+    .. py:attribute:: mode
+
+        The energy shifting/smoothing mode: Possible values are:
+        ``"none"``, ``"shift"``, and ``"xplor"``.
+
+        .. rubric:: Example
+
+        .. code-block:: python
+
+            lennard_jones.mode = 'shift'
+
+        Type: `str`
+    """
+    _cpp_class_name = "PairPotentialLennardJones"
+
+    def __init__(self, default_r_cut=None, default_r_on=0.0, mode='none'):
+        if default_r_cut is None:
+            default_r_cut = float
+        else:
+            default_r_cut = float(default_r_cut)
+
+        params = hoomd.data.typeparam.TypeParameter(
+            'params', 'particle_types',
+            hoomd.data.parameterdicts.TypeParameterDict(
+                epsilon=float,
+                sigma=float,
+                r_cut=default_r_cut,
+                r_on=float(default_r_on),
+                len_keys=2))
+        self._add_typeparam(params)
+
+        self._param_dict.update(
+            hoomd.data.parameterdicts.ParameterDict(
+                mode=hoomd.data.typeconverter.OnlyFrom(("none", "shift",
+                                                        "xplor"))))
+        self.mode = mode

From 9659523df01955b8a00530e1bfb0157e4d6f273b Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Wed, 12 Jun 2024 16:02:32 -0400
Subject: [PATCH 05/32] Create test_pair_expanded_gaussian.py

---
 .../pytest/test_pair_expanded_gaussian.py     | 216 ++++++++++++++++++
 1 file changed, 216 insertions(+)
 create mode 100644 hoomd/hpmc/pytest/test_pair_expanded_gaussian.py

diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
new file mode 100644
index 0000000000..65473d0ca0
--- /dev/null
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -0,0 +1,216 @@
+# Copyright (c) 2009-2024 The Regents of the University of Michigan.
+# Part of HOOMD-blue, released under the BSD 3-Clause License.
+
+"""Test hoomd.hpmc.pair.LennardJones and HPMC pair infrastructure."""
+
+import hoomd
+import pytest
+
+valid_constructor_args = [
+    {},
+    dict(default_r_cut=2.5),
+    dict(default_r_on=2.0),
+    dict(mode='shift'),
+]
+
+
+@pytest.mark.parametrize("constructor_args", valid_constructor_args)
+def test_valid_construction(device, constructor_args):
+    """Test that LennardJones can be constructed with valid arguments."""
+    hoomd.hpmc.pair.LennardJones(**constructor_args)
+
+
+@pytest.fixture(scope='session')
+def mc_simulation_factory(simulation_factory, two_particle_snapshot_factory):
+    """Make a MC simulation with two particles separate dy by a distance d."""
+
+    def make_simulation(d=1):
+        snapshot = two_particle_snapshot_factory(d=d)
+        simulation = simulation_factory(snapshot)
+
+        sphere = hoomd.hpmc.integrate.Sphere()
+        sphere.shape['A'] = dict(diameter=0)
+        simulation.operations.integrator = sphere
+
+        return simulation
+
+    return make_simulation
+
+
+@pytest.mark.cpu
+def test_attaching(mc_simulation_factory):
+    """Test that LennardJones attaches."""
+    lennard_jones = hoomd.hpmc.pair.LennardJones()
+    lennard_jones.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0, r_cut=2.5)
+
+    simulation = mc_simulation_factory()
+    simulation.operations.integrator.pair_potentials = [lennard_jones]
+    simulation.run(0)
+
+    assert simulation.operations.integrator._attached
+    assert lennard_jones._attached
+
+    simulation.operations.integrator.pair_potentials.remove(lennard_jones)
+    assert not lennard_jones._attached
+
+
+invalid_parameters = [
+    {},
+    dict(epsilon=1.0),
+    dict(epsilon=1.0, sigma=1.0),
+    dict(epsilon=1.0, sigma=1.0, r_cut='invalid'),
+    dict(epsilon=1.0, sigma=1.0, r_cut=2.5, r_on='invalid'),
+    dict(epsilon=1.0, sigma=1.0, r_cut=2.5, r_on=2.0, invalid=10),
+]
+
+
+@pytest.mark.parametrize("parameters", invalid_parameters)
+@pytest.mark.cpu
+def test_invalid_params_on_attach(mc_simulation_factory, parameters):
+    """Test that LennardJones validates parameters."""
+    lennard_jones = hoomd.hpmc.pair.LennardJones()
+    lennard_jones.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0, r_cut=2.5)
+
+    # Some parameters are validated only after attaching.
+    simulation = mc_simulation_factory()
+    simulation.operations.integrator.pair_potentials = [lennard_jones]
+    simulation.run(0)
+
+    with pytest.raises((
+            RuntimeError,
+            hoomd.error.TypeConversionError,
+            KeyError,
+    )):
+        lennard_jones.params[('A', 'A')] = parameters
+
+
+def xplor_factor(r, r_on, r_cut):
+    """Compute the XPLOR smoothing factor."""
+    if r < r_on:
+        return 1
+    if r < r_cut:
+        denominator = (r_cut**2 - r_on**2)**3
+        numerator = (r_cut**2 - r**2)**2 * (r_cut**2 + 2 * r**2 - 3 * r_on**2)
+        return numerator / denominator
+
+    return 0
+
+
+def lj(r, r_cut, epsilon, sigma):
+    """Compute the lj energy."""
+    return 4 * epsilon * ((sigma / r)**12 - (sigma / r)**6)
+
+
+# (pair params,
+#  distance between particles,
+#  expected energy)
+lennard_jones_test_parameters = [
+    (
+        dict(epsilon=2.0, sigma=1.5, r_cut=2.5),
+        'none',
+        3.0,
+        0.0,
+    ),
+    (
+        dict(epsilon=2.0, sigma=1.5, r_cut=2.5),
+        'none',
+        1.5 * 2**(1 / 6),
+        -2.0,
+    ),
+    (
+        dict(epsilon=3.0, sigma=0.5, r_cut=2.5),
+        'none',
+        0.5,
+        0,
+    ),
+    (
+        dict(epsilon=5.0, sigma=1.1, r_cut=2),
+        'none',
+        1.5,
+        lj(1.5, 2, 5, 1.1),
+    ),
+    (
+        dict(epsilon=5.0, sigma=1.1, r_cut=2),
+        'shift',
+        1.5,
+        lj(1.5, 2, 5, 1.1) - lj(2, 2, 5, 1.1),
+    ),
+    (
+        dict(epsilon=1.0, sigma=1, r_cut=3),
+        'shift',
+        3.2,
+        0,
+    ),
+    (
+        dict(epsilon=5.0, sigma=1.1, r_cut=2.5, r_on=2.0),
+        'xplor',
+        1.5,
+        lj(1.5, 2.5, 5.0, 1.1) * xplor_factor(1.5, 2.0, 2.5),
+    ),
+    (
+        dict(epsilon=5.0, sigma=1.1, r_cut=2.5, r_on=2.0),
+        'xplor',
+        2.3,
+        lj(2.3, 2.5, 5, 1.1) * xplor_factor(2.3, 2.0, 2.5),
+    ),
+    (
+        dict(epsilon=5.0, sigma=1.1, r_cut=2, r_on=3),
+        'xplor',
+        1.5,
+        lj(1.5, 2, 5, 1.1) - lj(2, 2, 5, 1.1),
+    ),
+    (
+        dict(epsilon=1.0, sigma=1, r_cut=3, r_on=4),
+        'xplor',
+        3.2,
+        0,
+    ),
+]
+
+
+@pytest.mark.parametrize('params, mode, d, expected_energy',
+                         lennard_jones_test_parameters)
+@pytest.mark.cpu
+def test_energy(mc_simulation_factory, params, mode, d, expected_energy):
+    """Test that LennardJones computes the correct energies for 1 pair."""
+    lennard_jones = hoomd.hpmc.pair.LennardJones(mode=mode)
+    lennard_jones.params[('A', 'A')] = params
+
+    simulation = mc_simulation_factory(d=d)
+    simulation.operations.integrator.pair_potentials = [lennard_jones]
+    simulation.run(0)
+
+    assert lennard_jones.energy == pytest.approx(expected=expected_energy,
+                                                 rel=1e-5)
+
+
+@pytest.mark.cpu
+def test_multiple_pair_potentials(mc_simulation_factory):
+    """Test that energy operates correctly with multiple pair potentials."""
+    lennard_jones_1 = hoomd.hpmc.pair.LennardJones()
+    lennard_jones_1.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0, r_cut=2.5)
+
+    lennard_jones_2 = hoomd.hpmc.pair.LennardJones()
+    lennard_jones_2.params[('A', 'A')] = dict(epsilon=2.0, sigma=1.0, r_cut=2.5)
+
+    # Some parameters are validated only after attaching.
+    simulation = mc_simulation_factory(2**(1 / 6))
+    simulation.operations.integrator.pair_potentials = [
+        lennard_jones_1, lennard_jones_2
+    ]
+    simulation.run(0)
+
+    assert lennard_jones_1.energy == pytest.approx(expected=-1.0, rel=1e-5)
+    assert lennard_jones_2.energy == pytest.approx(expected=-2.0, rel=1e-5)
+    assert simulation.operations.integrator.pair_energy == pytest.approx(
+        expected=-3.0, rel=1e-5)
+
+
+def test_logging():
+    hoomd.conftest.logging_check(
+        hoomd.hpmc.pair.LennardJones, ('hpmc', 'pair'), {
+            'energy': {
+                'category': hoomd.logging.LoggerCategories.scalar,
+                'default': True
+            }
+        })

From 891b7305e61f2325d5e4bff0ad23d7516e9d3620 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 10:35:21 -0400
Subject: [PATCH 06/32] Replace Lennard-Jones variable names and text with
 Expanded Gaussian

---
 hoomd/hpmc/PairPotentialExpandedGaussian.cc   | 24 +++----
 hoomd/hpmc/PairPotentialExpandedGaussian.h    | 10 +--
 hoomd/hpmc/pair/expanded_gaussian.py          | 24 +++----
 .../pytest/test_pair_expanded_gaussian.py     | 66 ++++++++++---------
 4 files changed, 65 insertions(+), 59 deletions(-)

diff --git a/hoomd/hpmc/PairPotentialExpandedGaussian.cc b/hoomd/hpmc/PairPotentialExpandedGaussian.cc
index 3cfaa585ad..00435227f6 100644
--- a/hoomd/hpmc/PairPotentialExpandedGaussian.cc
+++ b/hoomd/hpmc/PairPotentialExpandedGaussian.cc
@@ -1,19 +1,19 @@
 // Copyright (c) 2009-2024 The Regents of the University of Michigan.
 // Part of HOOMD-blue, released under the BSD 3-Clause License.
 
-#include "PairPotentialLennardJones.h"
+#include "PairPotentialExpandedGaussian.h"
 
 namespace hoomd
     {
 namespace hpmc
     {
 
-PairPotentialLennardJones::PairPotentialLennardJones(std::shared_ptr<SystemDefinition> sysdef)
+PairPotentialExpandedGaussian::PairPotentialExpandedGaussian(std::shared_ptr<SystemDefinition> sysdef)
     : PairPotential(sysdef), m_params(m_type_param_index.getNumElements())
     {
     }
 
-LongReal PairPotentialLennardJones::energy(const LongReal r_squared,
+LongReal PairPotentialExpandedGaussian::energy(const LongReal r_squared,
                                            const vec3<LongReal>& r_ij,
                                            const unsigned int type_i,
                                            const quat<LongReal>& q_i,
@@ -55,7 +55,7 @@ LongReal PairPotentialLennardJones::energy(const LongReal r_squared,
     return energy;
     }
 
-void PairPotentialLennardJones::setParamsPython(pybind11::tuple typ, pybind11::dict params)
+void PairPotentialExpandedGaussian::setParamsPython(pybind11::tuple typ, pybind11::dict params)
     {
     auto pdata = m_sysdef->getParticleData();
     auto type_i = pdata->getTypeByName(typ[0].cast<std::string>());
@@ -68,7 +68,7 @@ void PairPotentialLennardJones::setParamsPython(pybind11::tuple typ, pybind11::d
     notifyRCutChanged();
     }
 
-pybind11::dict PairPotentialLennardJones::getParamsPython(pybind11::tuple typ)
+pybind11::dict PairPotentialExpandedGaussian::getParamsPython(pybind11::tuple typ)
     {
     auto pdata = m_sysdef->getParticleData();
     auto type_i = pdata->getTypeByName(typ[0].cast<std::string>());
@@ -79,17 +79,17 @@ pybind11::dict PairPotentialLennardJones::getParamsPython(pybind11::tuple typ)
 
 namespace detail
     {
-void exportPairPotentialLennardJones(pybind11::module& m)
+void exportPairPotentialExpandedGaussian(pybind11::module& m)
     {
-    pybind11::class_<PairPotentialLennardJones,
+    pybind11::class_<PairPotentialExpandedGaussian,
                      PairPotential,
-                     std::shared_ptr<PairPotentialLennardJones>>(m, "PairPotentialLennardJones")
+                     std::shared_ptr<PairPotentialExpandedGaussian>>(m, "PairPotentialExpandedGaussian")
         .def(pybind11::init<std::shared_ptr<SystemDefinition>>())
-        .def("setParams", &PairPotentialLennardJones::setParamsPython)
-        .def("getParams", &PairPotentialLennardJones::getParamsPython)
+        .def("setParams", &PairPotentialExpandedGaussian::setParamsPython)
+        .def("getParams", &PairPotentialExpandedGaussian::getParamsPython)
         .def_property("mode",
-                      &PairPotentialLennardJones::getMode,
-                      &PairPotentialLennardJones::setMode);
+                      &PairPotentialExpandedGaussian::getMode,
+                      &PairPotentialExpandedGaussian::setMode);
     }
     } // end namespace detail
     } // end namespace hpmc
diff --git a/hoomd/hpmc/PairPotentialExpandedGaussian.h b/hoomd/hpmc/PairPotentialExpandedGaussian.h
index 50c43ee556..84cb5d73a8 100644
--- a/hoomd/hpmc/PairPotentialExpandedGaussian.h
+++ b/hoomd/hpmc/PairPotentialExpandedGaussian.h
@@ -10,15 +10,15 @@ namespace hoomd
 namespace hpmc
     {
 
-/*** Compute Lennard-Jones energy between two particles.
+/*** Compute Expanded Gaussian energy between two particles.
 
 For use with HPMC simulations.
 */
-class PairPotentialLennardJones : public hpmc::PairPotential
+class PairPotentialExpandedGaussian : public hpmc::PairPotential
     {
     public:
-    PairPotentialLennardJones(std::shared_ptr<SystemDefinition> sysdef);
-    virtual ~PairPotentialLennardJones() { }
+    PairPotentialExpandedGaussian(std::shared_ptr<SystemDefinition> sysdef);
+    virtual ~PairPotentialExpandedGaussian() { }
 
     virtual LongReal energy(const LongReal r_squared,
                             const vec3<LongReal>& r_ij,
@@ -87,7 +87,7 @@ class PairPotentialLennardJones : public hpmc::PairPotential
         xplor
         };
 
-    /// Type pair parameters of LJ potential
+    /// Type pair parameters of EG potential
     struct ParamType
         {
         ParamType()
diff --git a/hoomd/hpmc/pair/expanded_gaussian.py b/hoomd/hpmc/pair/expanded_gaussian.py
index ac680ab1e1..5ce65ff92e 100644
--- a/hoomd/hpmc/pair/expanded_gaussian.py
+++ b/hoomd/hpmc/pair/expanded_gaussian.py
@@ -1,7 +1,7 @@
 # Copyright (c) 2009-2024 The Regents of the University of Michigan.
 # Part of HOOMD-blue, released under the BSD 3-Clause License.
 
-"""Lennard-Jones pair potential.
+"""Expanded Gaussian pair potential.
 
 .. invisible-code-block: python
 
@@ -16,9 +16,9 @@
 from .pair import Pair
 
 
-@hoomd.logging.modify_namespace(('hpmc', 'pair', 'LennardJones'))
-class LennardJones(Pair):
-    """Lennard-Jones pair potential (HPMC).
+@hoomd.logging.modify_namespace(('hpmc', 'pair', 'ExpandedGaussian'))
+class ExpandedGaussian(Pair):
+    """Expanded Gaussian pair potential (HPMC).
 
     Args:
         default_r_cut (float): Default cutoff radius :math:`[\\mathrm{length}]`.
@@ -26,13 +26,13 @@ class LennardJones(Pair):
           :math:`[\\mathrm{length}]`.
         mode (str): Energy shifting/smoothing mode.
 
-    `LennardJones` computes the Lennard-Jones pair potential between every pair
+    `ExpandedGaussian` computes the Expanded Gaussian pair potential between every pair
     of particles in the simulation state. The functional form of the potential,
     including its behavior under shifting modes, is identical to that in
-    the MD pair potential `hoomd.md.pair.LJ`.
+    the MD pair potential `hoomd.md.pair.ExpandedGaussian`.
 
     See Also:
-        `hoomd.md.pair.LJ`
+        `hoomd.md.pair.ExpandedGaussian`
 
         `hoomd.md.pair`
 
@@ -40,9 +40,9 @@ class LennardJones(Pair):
 
     .. code-block:: python
 
-        lennard_jones =  hoomd.hpmc.pair.LennardJones()
-        lennard_jones.params[('A', 'A')] = dict(epsilon=1, sigma=1, r_cut=2.5)
-        simulation.operations.integrator.pair_potentials = [lennard_jones]
+        expanded_gaussian =  hoomd.hpmc.pair.ExpandedGaussian()
+        expanded_gaussian.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0, delta=1.0, r_cut=2.5)
+        simulation.operations.integrator.pair_potentials = [expanded_gaussian]
 
     .. py:attribute:: params
 
@@ -70,11 +70,11 @@ class LennardJones(Pair):
 
         .. code-block:: python
 
-            lennard_jones.mode = 'shift'
+            expanded_gaussian.mode = 'shift'
 
         Type: `str`
     """
-    _cpp_class_name = "PairPotentialLennardJones"
+    _cpp_class_name = "PairPotentialExpandedGaussian"
 
     def __init__(self, default_r_cut=None, default_r_on=0.0, mode='none'):
         if default_r_cut is None:
diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index 65473d0ca0..ff2242002e 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -1,7 +1,7 @@
 # Copyright (c) 2009-2024 The Regents of the University of Michigan.
 # Part of HOOMD-blue, released under the BSD 3-Clause License.
 
-"""Test hoomd.hpmc.pair.LennardJones and HPMC pair infrastructure."""
+"""Test hoomd.hpmc.pair.ExpandedGaussian and HPMC pair infrastructure."""
 
 import hoomd
 import pytest
@@ -39,19 +39,19 @@ def make_simulation(d=1):
 
 @pytest.mark.cpu
 def test_attaching(mc_simulation_factory):
-    """Test that LennardJones attaches."""
-    lennard_jones = hoomd.hpmc.pair.LennardJones()
-    lennard_jones.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0, r_cut=2.5)
+    """Test that ExpandedGaussian attaches."""
+    expanded_gauss = hoomd.hpmc.pair.ExpandedGaussian()
+    expanded_gauss.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0, r_cut=2.5)
 
     simulation = mc_simulation_factory()
-    simulation.operations.integrator.pair_potentials = [lennard_jones]
+    simulation.operations.integrator.pair_potentials = [expanded_gauss]
     simulation.run(0)
 
     assert simulation.operations.integrator._attached
-    assert lennard_jones._attached
+    assert expanded_gauss._attached
 
-    simulation.operations.integrator.pair_potentials.remove(lennard_jones)
-    assert not lennard_jones._attached
+    simulation.operations.integrator.pair_potentials.remove(expanded_gauss)
+    assert not expanded_gauss._attached
 
 
 invalid_parameters = [
@@ -67,13 +67,13 @@ def test_attaching(mc_simulation_factory):
 @pytest.mark.parametrize("parameters", invalid_parameters)
 @pytest.mark.cpu
 def test_invalid_params_on_attach(mc_simulation_factory, parameters):
-    """Test that LennardJones validates parameters."""
-    lennard_jones = hoomd.hpmc.pair.LennardJones()
-    lennard_jones.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0, r_cut=2.5)
+    """Test that ExpandedGaussian validates parameters."""
+    expanded_gauss = hoomd.hpmc.pair.ExpandedGaussian()
+    expanded_gauss.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0, r_cut=2.5)
 
     # Some parameters are validated only after attaching.
     simulation = mc_simulation_factory()
-    simulation.operations.integrator.pair_potentials = [lennard_jones]
+    simulation.operations.integrator.pair_potentials = [expanded_gauss]
     simulation.run(0)
 
     with pytest.raises((
@@ -81,7 +81,7 @@ def test_invalid_params_on_attach(mc_simulation_factory, parameters):
             hoomd.error.TypeConversionError,
             KeyError,
     )):
-        lennard_jones.params[('A', 'A')] = parameters
+        expanded_gauss.params[('A', 'A')] = parameters
 
 
 def xplor_factor(r, r_on, r_cut):
@@ -104,7 +104,7 @@ def lj(r, r_cut, epsilon, sigma):
 # (pair params,
 #  distance between particles,
 #  expected energy)
-lennard_jones_test_parameters = [
+expanded_gauss_test_parameters = [
     (
         dict(epsilon=2.0, sigma=1.5, r_cut=2.5),
         'none',
@@ -169,46 +169,52 @@ def lj(r, r_cut, epsilon, sigma):
 
 
 @pytest.mark.parametrize('params, mode, d, expected_energy',
-                         lennard_jones_test_parameters)
+                         expanded_gauss_test_parameters)
 @pytest.mark.cpu
 def test_energy(mc_simulation_factory, params, mode, d, expected_energy):
-    """Test that LennardJones computes the correct energies for 1 pair."""
-    lennard_jones = hoomd.hpmc.pair.LennardJones(mode=mode)
-    lennard_jones.params[('A', 'A')] = params
+    """Test that ExpandedGaussian computes the correct energies for 1 pair."""
+    expanded_gauss = hoomd.hpmc.pair.ExpandedGaussian(mode=mode)
+    expanded_gauss.params[('A', 'A')] = params
 
     simulation = mc_simulation_factory(d=d)
-    simulation.operations.integrator.pair_potentials = [lennard_jones]
+    simulation.operations.integrator.pair_potentials = [expanded_gauss]
     simulation.run(0)
 
-    assert lennard_jones.energy == pytest.approx(expected=expected_energy,
-                                                 rel=1e-5)
+    assert expanded_gauss.energy == pytest.approx(expected=expected_energy,
+                                                  rel=1e-5)
 
 
 @pytest.mark.cpu
 def test_multiple_pair_potentials(mc_simulation_factory):
     """Test that energy operates correctly with multiple pair potentials."""
-    lennard_jones_1 = hoomd.hpmc.pair.LennardJones()
-    lennard_jones_1.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0, r_cut=2.5)
-
-    lennard_jones_2 = hoomd.hpmc.pair.LennardJones()
-    lennard_jones_2.params[('A', 'A')] = dict(epsilon=2.0, sigma=1.0, r_cut=2.5)
+    expanded_gauss_1 = hoomd.hpmc.pair.ExpandedGaussian()
+    expanded_gauss_1.params[('A', 'A')] = dict(epsilon=1.0,
+                                               sigma=1.0,
+                                               delta=1.0,
+                                               r_cut=2.5)
+
+    expanded_gauss_2 = hoomd.hpmc.pair.ExpandedGaussian()
+    expanded_gauss_2.params[('A', 'A')] = dict(epsilon=2.0,
+                                               sigma=1.0,
+                                               delta=1.0,
+                                               r_cut=2.5)
 
     # Some parameters are validated only after attaching.
     simulation = mc_simulation_factory(2**(1 / 6))
     simulation.operations.integrator.pair_potentials = [
-        lennard_jones_1, lennard_jones_2
+        expanded_gauss_1, expanded_gauss_2
     ]
     simulation.run(0)
 
-    assert lennard_jones_1.energy == pytest.approx(expected=-1.0, rel=1e-5)
-    assert lennard_jones_2.energy == pytest.approx(expected=-2.0, rel=1e-5)
+    assert expanded_gauss_1.energy == pytest.approx(expected=-1.0, rel=1e-5)
+    assert expanded_gauss_2.energy == pytest.approx(expected=-2.0, rel=1e-5)
     assert simulation.operations.integrator.pair_energy == pytest.approx(
         expected=-3.0, rel=1e-5)
 
 
 def test_logging():
     hoomd.conftest.logging_check(
-        hoomd.hpmc.pair.LennardJones, ('hpmc', 'pair'), {
+        hoomd.hpmc.pair.ExpandedGaussian, ('hpmc', 'pair'), {
             'energy': {
                 'category': hoomd.logging.LoggerCategories.scalar,
                 'default': True

From 5cac99a153278e6aa958005aa12f8d08d9df2675 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 10:37:01 -0400
Subject: [PATCH 07/32] Update C++ class header params to match EG definition

---
 hoomd/hpmc/PairPotentialExpandedGaussian.h | 33 ++++++++++++++--------
 1 file changed, 22 insertions(+), 11 deletions(-)

diff --git a/hoomd/hpmc/PairPotentialExpandedGaussian.h b/hoomd/hpmc/PairPotentialExpandedGaussian.h
index 84cb5d73a8..5926fd3a46 100644
--- a/hoomd/hpmc/PairPotentialExpandedGaussian.h
+++ b/hoomd/hpmc/PairPotentialExpandedGaussian.h
@@ -33,7 +33,7 @@ class PairPotentialExpandedGaussian : public hpmc::PairPotential
     virtual LongReal computeRCutNonAdditive(unsigned int type_i, unsigned int type_j) const
         {
         unsigned int param_index = m_type_param_index(type_i, type_j);
-        return slow::sqrt(m_params[param_index].r_cut_squared);
+        return m_params[param_index].r_cut;
         }
 
     /// Set type pair dependent parameters to the potential.
@@ -92,8 +92,11 @@ class PairPotentialExpandedGaussian : public hpmc::PairPotential
         {
         ParamType()
             {
-            sigma_6 = 0;
-            epsilon_x_4 = 0;
+            sigma_2 = 0;
+            epsilon = 0;
+            delta = 0;
+            r_cut = 0;
+            r_on = 0;
             r_cut_squared = 0;
             r_on_squared = 0;
             }
@@ -102,11 +105,15 @@ class PairPotentialExpandedGaussian : public hpmc::PairPotential
             {
             auto sigma(v["sigma"].cast<LongReal>());
             auto epsilon(v["epsilon"].cast<LongReal>());
+            auto delta(v["delta"].cast<LongReal>());
             auto r_cut(v["r_cut"].cast<LongReal>());
             auto r_on(v["r_on"].cast<LongReal>());
 
-            sigma_6 = sigma * sigma * sigma * sigma * sigma * sigma;
-            epsilon_x_4 = LongReal(4.0) * epsilon;
+            sigma_2 = sigma * sigma;
+            epsilon = epsilon;
+            delta = delta;
+            r_cut = r_cut;
+            r_on = r_on;
             r_cut_squared = r_cut * r_cut;
             r_on_squared = r_on * r_on;
             }
@@ -115,16 +122,20 @@ class PairPotentialExpandedGaussian : public hpmc::PairPotential
             {
             pybind11::dict result;
 
-            result["sigma"] = pow(sigma_6, 1. / 6.);
-            result["epsilon"] = epsilon_x_4 / 4.0;
-            result["r_cut"] = slow::sqrt(r_cut_squared);
-            result["r_on"] = slow::sqrt(r_on_squared);
+            result["sigma"] = pow(sigma_2, 1. / 2.);
+            result["epsilon"] = epsilon;
+            result["delta"] = delta;
+            result["r_cut"] = r_cut;
+            result["r_on"] = r_on;
 
             return result;
             }
 
-        LongReal sigma_6;
-        LongReal epsilon_x_4;
+        LongReal sigma_2;
+        LongReal epsilon;
+        LongReal delta;
+        LongReal r_cut;
+        LongReal r_on;
         LongReal r_cut_squared;
         LongReal r_on_squared;
         };

From b684f3351cc867bf16619598c114e7a3997d0bca Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 10:37:37 -0400
Subject: [PATCH 08/32] Update C++ class energy method to match EG definition

---
 hoomd/hpmc/PairPotentialExpandedGaussian.cc | 17 +++++++----------
 1 file changed, 7 insertions(+), 10 deletions(-)

diff --git a/hoomd/hpmc/PairPotentialExpandedGaussian.cc b/hoomd/hpmc/PairPotentialExpandedGaussian.cc
index 00435227f6..9172678c1d 100644
--- a/hoomd/hpmc/PairPotentialExpandedGaussian.cc
+++ b/hoomd/hpmc/PairPotentialExpandedGaussian.cc
@@ -25,19 +25,16 @@ LongReal PairPotentialExpandedGaussian::energy(const LongReal r_squared,
     unsigned int param_index = m_type_param_index(type_i, type_j);
     const auto& param = m_params[param_index];
 
-    LongReal lj2 = param.epsilon_x_4 * param.sigma_6;
-    LongReal lj1 = lj2 * param.sigma_6;
-
-    LongReal r_2_inverse = LongReal(1.0) / r_squared;
-    LongReal r_6_inverse = r_2_inverse * r_2_inverse * r_2_inverse;
-
-    LongReal energy = r_6_inverse * (lj1 * r_6_inverse - lj2);
+    LongReal r = fast::sqrt(r_squared); // Scalar or LongReal?
+    LongReal num = r - param.delta;
+    LongReal exponent = (LongReal(-0.5) * (num * num)) / param.sigma_2;
+    LongReal energy = param.epsilon * fast::exp(exponent);
 
     if (m_mode == shift || (m_mode == xplor && param.r_on_squared >= param.r_cut_squared))
         {
-        LongReal r_cut_2_inverse = LongReal(1.0) / param.r_cut_squared;
-        LongReal r_cut_6_inverse = r_cut_2_inverse * r_cut_2_inverse * r_cut_2_inverse;
-        energy -= r_cut_6_inverse * (lj1 * r_cut_6_inverse - lj2);
+        LongReal num_r_cut = param.r_cut - param.delta;
+        LongReal exponent_r_cut = (LongReal(-0.5) * (num_r_cut * num_r_cut) / param.sigma_2);
+        energy -= param.epsilon * fast::exp(exponent_r_cut);
         }
 
     if (m_mode == xplor && r_squared > param.r_on_squared)

From a24f6e4c7cd44fd0aea2902ac610fa4ca37ebb15 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 10:38:07 -0400
Subject: [PATCH 09/32] Update python class docstring

---
 hoomd/hpmc/pair/expanded_gaussian.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/hoomd/hpmc/pair/expanded_gaussian.py b/hoomd/hpmc/pair/expanded_gaussian.py
index 5ce65ff92e..2638b9b7a6 100644
--- a/hoomd/hpmc/pair/expanded_gaussian.py
+++ b/hoomd/hpmc/pair/expanded_gaussian.py
@@ -53,6 +53,9 @@ class ExpandedGaussian(Pair):
         * ``sigma`` (`float`, **required**) -
           Characteristic length scale :math:`\\sigma`
           :math:`[\\mathrm{length}]`.
+        * ``delta`` (`float`, **required**) -
+          Characteristic length scale :math:`\\delta`
+          :math:`[\\mathrm{length}]`.
         * ``r_cut`` (`float`): Cutoff radius :math:`[\\mathrm{length}]`.
           Defaults to the value given in ``default_r_cut`` on construction.
         * ``r_on`` (`float`): XPLOR on radius :math:`[\\mathrm{length}]`.

From b8ac72ed5ef8730cb6ee11efab475674b91c7e93 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 10:38:35 -0400
Subject: [PATCH 10/32] Update python class params to match EG definition

---
 hoomd/hpmc/pair/expanded_gaussian.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/hoomd/hpmc/pair/expanded_gaussian.py b/hoomd/hpmc/pair/expanded_gaussian.py
index 2638b9b7a6..8edfa84672 100644
--- a/hoomd/hpmc/pair/expanded_gaussian.py
+++ b/hoomd/hpmc/pair/expanded_gaussian.py
@@ -90,6 +90,7 @@ def __init__(self, default_r_cut=None, default_r_on=0.0, mode='none'):
             hoomd.data.parameterdicts.TypeParameterDict(
                 epsilon=float,
                 sigma=float,
+                delta=float,
                 r_cut=default_r_cut,
                 r_on=float(default_r_on),
                 len_keys=2))

From 32238db013b1dd58d9766be1146365a5063b6ef8 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 10:39:20 -0400
Subject: [PATCH 11/32] Remove duplicate test

---
 .../pytest/test_pair_expanded_gaussian.py     | 22 ++-----------------
 1 file changed, 2 insertions(+), 20 deletions(-)

diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index ff2242002e..d54bf43139 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -16,26 +16,8 @@
 
 @pytest.mark.parametrize("constructor_args", valid_constructor_args)
 def test_valid_construction(device, constructor_args):
-    """Test that LennardJones can be constructed with valid arguments."""
-    hoomd.hpmc.pair.LennardJones(**constructor_args)
-
-
-@pytest.fixture(scope='session')
-def mc_simulation_factory(simulation_factory, two_particle_snapshot_factory):
-    """Make a MC simulation with two particles separate dy by a distance d."""
-
-    def make_simulation(d=1):
-        snapshot = two_particle_snapshot_factory(d=d)
-        simulation = simulation_factory(snapshot)
-
-        sphere = hoomd.hpmc.integrate.Sphere()
-        sphere.shape['A'] = dict(diameter=0)
-        simulation.operations.integrator = sphere
-
-        return simulation
-
-    return make_simulation
-
+    """Test that ExpandedGaussian can be constructed with valid arguments."""
+    hoomd.hpmc.pair.ExpandedGaussian(**constructor_args)
 
 @pytest.mark.cpu
 def test_attaching(mc_simulation_factory):

From 4321731209267b9ec487b57b52d706cde5c25c9b Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 10:39:37 -0400
Subject: [PATCH 12/32] Update invalid params for tests

---
 hoomd/hpmc/pytest/test_pair_expanded_gaussian.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index d54bf43139..cd76f1bd35 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -40,9 +40,9 @@ def test_attaching(mc_simulation_factory):
     {},
     dict(epsilon=1.0),
     dict(epsilon=1.0, sigma=1.0),
-    dict(epsilon=1.0, sigma=1.0, r_cut='invalid'),
-    dict(epsilon=1.0, sigma=1.0, r_cut=2.5, r_on='invalid'),
-    dict(epsilon=1.0, sigma=1.0, r_cut=2.5, r_on=2.0, invalid=10),
+    dict(epsilon=1.0, sigma=1.0, delta=1.0, r_cut='invalid'),
+    dict(epsilon=1.0, sigma=1.0, delta=1.0, r_cut=2.5, r_on='invalid'),
+    dict(epsilon=1.0, sigma=1.0, delta=1.0, r_cut=2.5, r_on=2.0, invalid=10),
 ]
 
 

From 35e1ae81020d60490c0557db938abbc93334f7ac Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 10:40:32 -0400
Subject: [PATCH 13/32] Update tests' energy calculation to match EG definition

---
 hoomd/hpmc/pytest/test_pair_expanded_gaussian.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index cd76f1bd35..760ad75ca6 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -5,6 +5,7 @@
 
 import hoomd
 import pytest
+import math
 
 valid_constructor_args = [
     {},
@@ -78,9 +79,9 @@ def xplor_factor(r, r_on, r_cut):
     return 0
 
 
-def lj(r, r_cut, epsilon, sigma):
-    """Compute the lj energy."""
-    return 4 * epsilon * ((sigma / r)**12 - (sigma / r)**6)
+def eg(r, epsilon, sigma, delta):
+    """Compute the eg energy."""
+    return epsilon * math.exp(-0.5*(((r-delta)/sigma)**2))
 
 
 # (pair params,

From 494e0ed2f8220600ba99b997cbe25a61cade88d7 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 10:41:20 -0400
Subject: [PATCH 14/32] Update test params to suit EG definition

---
 .../pytest/test_pair_expanded_gaussian.py     | 43 +++++++------------
 1 file changed, 16 insertions(+), 27 deletions(-)

diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index 760ad75ca6..636426d469 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -85,65 +85,54 @@ def eg(r, epsilon, sigma, delta):
 
 
 # (pair params,
+#  mode,
 #  distance between particles,
 #  expected energy)
 expanded_gauss_test_parameters = [
     (
-        dict(epsilon=2.0, sigma=1.5, r_cut=2.5),
+        dict(epsilon=2.0, sigma=1.5, delta=1.0, r_cut=2.5),
         'none',
         3.0,
         0.0,
     ),
     (
-        dict(epsilon=2.0, sigma=1.5, r_cut=2.5),
+        dict(epsilon=2.0, sigma=1.5, delta=1.0, r_cut=2.5),
         'none',
-        1.5 * 2**(1 / 6),
-        -2.0,
+        1.0,
+        2.0,
     ),
     (
-        dict(epsilon=3.0, sigma=0.5, r_cut=2.5),
-        'none',
-        0.5,
-        0,
-    ),
-    (
-        dict(epsilon=5.0, sigma=1.1, r_cut=2),
+        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=6.0),
         'none',
         1.5,
-        lj(1.5, 2, 5, 1.1),
+        eg(1.5, 5, 1.1, 1.0),
     ),
     (
-        dict(epsilon=5.0, sigma=1.1, r_cut=2),
+        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=6.0),
         'shift',
         1.5,
-        lj(1.5, 2, 5, 1.1) - lj(2, 2, 5, 1.1),
-    ),
-    (
-        dict(epsilon=1.0, sigma=1, r_cut=3),
-        'shift',
-        3.2,
-        0,
+        eg(1.5, 5, 1.1, 1.0) - eg(2, 5, 1.1, 1.0),
     ),
     (
-        dict(epsilon=5.0, sigma=1.1, r_cut=2.5, r_on=2.0),
+        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=6.0, r_on=2.0),
         'xplor',
         1.5,
-        lj(1.5, 2.5, 5.0, 1.1) * xplor_factor(1.5, 2.0, 2.5),
+        eg(1.5, 5.0, 1.1, 1.0) * xplor_factor(1.5, 2.0, 6.0),
     ),
     (
-        dict(epsilon=5.0, sigma=1.1, r_cut=2.5, r_on=2.0),
+        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=6.0, r_on=2.0),
         'xplor',
         2.3,
-        lj(2.3, 2.5, 5, 1.1) * xplor_factor(2.3, 2.0, 2.5),
+        eg(2.3, 5, 1.1, 1.0) * xplor_factor(2.3, 2.0, 6.0),
     ),
     (
-        dict(epsilon=5.0, sigma=1.1, r_cut=2, r_on=3),
+        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=2.0, r_on=3),
         'xplor',
         1.5,
-        lj(1.5, 2, 5, 1.1) - lj(2, 2, 5, 1.1),
+        eg(1.5, 5, 1.1, 1.0) - eg(2, 5, 1.1, 1.0),
     ),
     (
-        dict(epsilon=1.0, sigma=1, r_cut=3, r_on=4),
+        dict(epsilon=1.0, sigma=1, delta=1.0, r_cut=3.0, r_on=4),
         'xplor',
         3.2,
         0,

From 40a4117237c775ab2867ad537d5881d921995c4d Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 11:38:54 -0400
Subject: [PATCH 15/32] Add EG import to __init__.py

---
 hoomd/hpmc/pair/__init__.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/hoomd/hpmc/pair/__init__.py b/hoomd/hpmc/pair/__init__.py
index db1fe0bbe2..ffb72ccbb1 100644
--- a/hoomd/hpmc/pair/__init__.py
+++ b/hoomd/hpmc/pair/__init__.py
@@ -30,6 +30,7 @@
 from . import user
 from .pair import Pair
 from .lennard_jones import LennardJones
+from .expanded_gaussian import ExpandedGaussian
 from .union import Union
 from .angular_step import AngularStep
 from .step import Step

From a0600c123653df08b6e939f7148ad7e2b23df10c Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 11:43:27 -0400
Subject: [PATCH 16/32] Add EG to documentation

---
 sphinx-doc/module-hpmc-pair.rst | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/sphinx-doc/module-hpmc-pair.rst b/sphinx-doc/module-hpmc-pair.rst
index af00e3ce67..05e643ffc8 100644
--- a/sphinx-doc/module-hpmc-pair.rst
+++ b/sphinx-doc/module-hpmc-pair.rst
@@ -12,6 +12,7 @@ hoomd.hpmc.pair
     :nosignatures:
 
     AngularStep
+    ExpandedGaussian
     LennardJones
     Pair
     Step
@@ -23,6 +24,7 @@ hoomd.hpmc.pair
     :synopsis: Pair potentials for HPMC.
     :show-inheritance:
     :members: AngularStep,
+              ExpandedGaussian,
               LennardJones,
               Pair,
               Step,

From 6452645f9f2ef9d013bebbba063853178a107360 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 14:34:13 -0400
Subject: [PATCH 17/32] Add missing auxiliary function

---
 hoomd/hpmc/pytest/test_pair_expanded_gaussian.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index 636426d469..81c97c6368 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -19,6 +19,22 @@
 def test_valid_construction(device, constructor_args):
     """Test that ExpandedGaussian can be constructed with valid arguments."""
     hoomd.hpmc.pair.ExpandedGaussian(**constructor_args)
+    
+@pytest.fixture(scope='session')
+def mc_simulation_factory(simulation_factory, two_particle_snapshot_factory):
+    """Make a MC simulation with two particles separate dy by a distance d."""
+
+    def make_simulation(d=1):
+        snapshot = two_particle_snapshot_factory(d=d)
+        simulation = simulation_factory(snapshot)
+
+        sphere = hoomd.hpmc.integrate.Sphere()
+        sphere.shape['A'] = dict(diameter=0)
+        simulation.operations.integrator = sphere
+
+        return simulation
+
+    return make_simulation
 
 @pytest.mark.cpu
 def test_attaching(mc_simulation_factory):

From edfca6e5abec1fc5c1a230f20d960b7fbe251ac2 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 14:34:30 -0400
Subject: [PATCH 18/32] Add missing export

---
 hoomd/hpmc/module.cc | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/hoomd/hpmc/module.cc b/hoomd/hpmc/module.cc
index 8869164f6c..7e0ee1dae4 100644
--- a/hoomd/hpmc/module.cc
+++ b/hoomd/hpmc/module.cc
@@ -46,6 +46,8 @@ void exportPairPotential(pybind11::module& m);
 
 void exportPairPotentialLennardJones(pybind11::module& m);
 
+void exportPairPotentialExpandedGaussian(pybind11::module& m);
+
 void exportPairPotentialAngularStep(pybind11::module& m);
 
 void exportPairPotentialStep(pybind11::module& m);

From 6a42b659b588a0e5d2d8e0c9416394da6a13c9d9 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 14:38:21 -0400
Subject: [PATCH 19/32] Fix missing export

---
 hoomd/hpmc/module.cc | 1 +
 1 file changed, 1 insertion(+)

diff --git a/hoomd/hpmc/module.cc b/hoomd/hpmc/module.cc
index 7e0ee1dae4..4b1f179e54 100644
--- a/hoomd/hpmc/module.cc
+++ b/hoomd/hpmc/module.cc
@@ -150,6 +150,7 @@ PYBIND11_MODULE(_hpmc, m)
 
     exportPairPotential(m);
     exportPairPotentialLennardJones(m);
+    exportPairPotentialExpandedGaussian(m);
     exportPairPotentialAngularStep(m);
     exportPairPotentialStep(m);
     exportPairPotentialUnion(m);

From 1d330323ee84334b90253b799331c5b3b2960619 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 14:42:07 -0400
Subject: [PATCH 20/32] Add missing parameter in test

---
 hoomd/hpmc/pytest/test_pair_expanded_gaussian.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index 81c97c6368..7a795f312b 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -40,7 +40,10 @@ def make_simulation(d=1):
 def test_attaching(mc_simulation_factory):
     """Test that ExpandedGaussian attaches."""
     expanded_gauss = hoomd.hpmc.pair.ExpandedGaussian()
-    expanded_gauss.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0, r_cut=2.5)
+    expanded_gauss.params[('A', 'A')] = dict(epsilon=1.0, 
+                                             sigma=1.0, 
+                                             delta=1.0, 
+                                             r_cut=2.5)
 
     simulation = mc_simulation_factory()
     simulation.operations.integrator.pair_potentials = [expanded_gauss]

From c1dac23386e92752ea790e5b115035ffe04b3a1d Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 14:43:55 -0400
Subject: [PATCH 21/32] Add missing parameter to test

---
 hoomd/hpmc/pytest/test_pair_expanded_gaussian.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index 7a795f312b..514731d1a5 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -71,7 +71,10 @@ def test_attaching(mc_simulation_factory):
 def test_invalid_params_on_attach(mc_simulation_factory, parameters):
     """Test that ExpandedGaussian validates parameters."""
     expanded_gauss = hoomd.hpmc.pair.ExpandedGaussian()
-    expanded_gauss.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0, r_cut=2.5)
+    expanded_gauss.params[('A', 'A')] = dict(epsilon=1.0,
+                                             sigma=1.0,
+                                             delta=1.0,
+                                             r_cut=2.5)
 
     # Some parameters are validated only after attaching.
     simulation = mc_simulation_factory()

From 61ebd04b1f20d4bb92e9b665591c11043891a419 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 16:35:24 -0400
Subject: [PATCH 22/32] Refactor pre-computed class attributes to match LJ

---
 hoomd/hpmc/PairPotentialExpandedGaussian.h | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/hoomd/hpmc/PairPotentialExpandedGaussian.h b/hoomd/hpmc/PairPotentialExpandedGaussian.h
index 5926fd3a46..63cd15f6af 100644
--- a/hoomd/hpmc/PairPotentialExpandedGaussian.h
+++ b/hoomd/hpmc/PairPotentialExpandedGaussian.h
@@ -33,7 +33,7 @@ class PairPotentialExpandedGaussian : public hpmc::PairPotential
     virtual LongReal computeRCutNonAdditive(unsigned int type_i, unsigned int type_j) const
         {
         unsigned int param_index = m_type_param_index(type_i, type_j);
-        return m_params[param_index].r_cut;
+        return slow::sqrt(m_params[param_index].r_cut_squared);
         }
 
     /// Set type pair dependent parameters to the potential.
@@ -125,8 +125,8 @@ class PairPotentialExpandedGaussian : public hpmc::PairPotential
             result["sigma"] = pow(sigma_2, 1. / 2.);
             result["epsilon"] = epsilon;
             result["delta"] = delta;
-            result["r_cut"] = r_cut;
-            result["r_on"] = r_on;
+            result["r_cut"] = slow::sqrt(r_cut_squared);
+            result["r_on"] = slow::sqrt(r_on_squared);
 
             return result;
             }
@@ -134,8 +134,6 @@ class PairPotentialExpandedGaussian : public hpmc::PairPotential
         LongReal sigma_2;
         LongReal epsilon;
         LongReal delta;
-        LongReal r_cut;
-        LongReal r_on;
         LongReal r_cut_squared;
         LongReal r_on_squared;
         };

From 6f6c04d31eba17b0064c8d5983e20836a91bbdae Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 16:37:30 -0400
Subject: [PATCH 23/32] Refactor energy calculation to match MD

---
 hoomd/hpmc/PairPotentialExpandedGaussian.cc | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/hoomd/hpmc/PairPotentialExpandedGaussian.cc b/hoomd/hpmc/PairPotentialExpandedGaussian.cc
index 9172678c1d..c0c8ead13f 100644
--- a/hoomd/hpmc/PairPotentialExpandedGaussian.cc
+++ b/hoomd/hpmc/PairPotentialExpandedGaussian.cc
@@ -25,16 +25,18 @@ LongReal PairPotentialExpandedGaussian::energy(const LongReal r_squared,
     unsigned int param_index = m_type_param_index(type_i, type_j);
     const auto& param = m_params[param_index];
 
-    LongReal r = fast::sqrt(r_squared); // Scalar or LongReal?
-    LongReal num = r - param.delta;
-    LongReal exponent = (LongReal(-0.5) * (num * num)) / param.sigma_2;
-    LongReal energy = param.epsilon * fast::exp(exponent);
+    LongReal r = fast::sqrt(r_squared);
+    LongReal rmd_2 = (r - param.delta) * (r - param.delta);
+    LongReal rmd_over_sigma_2 = rmd_2 / param.sigma_2;
+    LongReal exp_val = fast::exp(-LongReal(1.0) / LongReal(2.0) * rmd_over_sigma_2);
+    LongReal energy = param.epsilon * exp_val;
 
     if (m_mode == shift || (m_mode == xplor && param.r_on_squared >= param.r_cut_squared))
         {
-        LongReal num_r_cut = param.r_cut - param.delta;
-        LongReal exponent_r_cut = (LongReal(-0.5) * (num_r_cut * num_r_cut) / param.sigma_2);
-        energy -= param.epsilon * fast::exp(exponent_r_cut);
+        LongReal r_cut = fast::sqrt(param.r_cut_squared);
+        LongReal rcutmd_2 = (r_cut - param.delta) * (r_cut - param.delta);
+        LongReal rcutmd_over_sigma_2 = rcutmd_2 / param.sigma_2;
+        energy -= param.epsilon * fast::exp(-LongReal(1.0) / LongReal(2.0) * rcutmd_over_sigma_2);
         }
 
     if (m_mode == xplor && r_squared > param.r_on_squared)

From 120e78c5a02061b51ee26bb60b4e79a0f2162822 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 16:38:27 -0400
Subject: [PATCH 24/32] Fix typo in single pair potential test

---
 hoomd/hpmc/pytest/test_pair_expanded_gaussian.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index 514731d1a5..4052ce8b03 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -133,7 +133,7 @@ def eg(r, epsilon, sigma, delta):
         dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=6.0),
         'shift',
         1.5,
-        eg(1.5, 5, 1.1, 1.0) - eg(2, 5, 1.1, 1.0),
+        eg(1.5, 5, 1.1, 1.0) - eg(6.0, 5, 1.1, 1.0),
     ),
     (
         dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=6.0, r_on=2.0),

From b8df89dfc325de9c7ac15b4627ffa02b09234bcd Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 16:39:02 -0400
Subject: [PATCH 25/32] Revise multi-potential test to match EG definition

---
 hoomd/hpmc/pytest/test_pair_expanded_gaussian.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index 4052ce8b03..052cf98f7b 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -194,16 +194,18 @@ def test_multiple_pair_potentials(mc_simulation_factory):
                                                r_cut=2.5)
 
     # Some parameters are validated only after attaching.
-    simulation = mc_simulation_factory(2**(1 / 6))
+    simulation = mc_simulation_factory(1.5)
     simulation.operations.integrator.pair_potentials = [
         expanded_gauss_1, expanded_gauss_2
     ]
     simulation.run(0)
 
-    assert expanded_gauss_1.energy == pytest.approx(expected=-1.0, rel=1e-5)
-    assert expanded_gauss_2.energy == pytest.approx(expected=-2.0, rel=1e-5)
+    assert expanded_gauss_1.energy == pytest.approx(
+        expected=eg(1.5, 1.0, 1.0, 1.0), rel=1e-5)
+    assert expanded_gauss_2.energy == pytest.approx(
+        expected=eg(1.5, 2.0, 1.0, 1.0), rel=1e-5)
     assert simulation.operations.integrator.pair_energy == pytest.approx(
-        expected=-3.0, rel=1e-5)
+        expected=eg(1.5, 1.0, 1.0, 1.0) + eg(1.5, 2.0, 1.0, 1.0), rel=1e-5)
 
 
 def test_logging():

From 0d11fdb29c9c0aedf6a71c9937038ce69a05b5b0 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 16:39:24 -0400
Subject: [PATCH 26/32] Fix variable scoping error

---
 hoomd/hpmc/PairPotentialExpandedGaussian.h | 10 ++--------
 1 file changed, 2 insertions(+), 8 deletions(-)

diff --git a/hoomd/hpmc/PairPotentialExpandedGaussian.h b/hoomd/hpmc/PairPotentialExpandedGaussian.h
index 63cd15f6af..dd237d8826 100644
--- a/hoomd/hpmc/PairPotentialExpandedGaussian.h
+++ b/hoomd/hpmc/PairPotentialExpandedGaussian.h
@@ -95,8 +95,6 @@ class PairPotentialExpandedGaussian : public hpmc::PairPotential
             sigma_2 = 0;
             epsilon = 0;
             delta = 0;
-            r_cut = 0;
-            r_on = 0;
             r_cut_squared = 0;
             r_on_squared = 0;
             }
@@ -104,16 +102,12 @@ class PairPotentialExpandedGaussian : public hpmc::PairPotential
         ParamType(pybind11::dict v)
             {
             auto sigma(v["sigma"].cast<LongReal>());
-            auto epsilon(v["epsilon"].cast<LongReal>());
-            auto delta(v["delta"].cast<LongReal>());
             auto r_cut(v["r_cut"].cast<LongReal>());
             auto r_on(v["r_on"].cast<LongReal>());
 
             sigma_2 = sigma * sigma;
-            epsilon = epsilon;
-            delta = delta;
-            r_cut = r_cut;
-            r_on = r_on;
+            epsilon = v["epsilon"].cast<LongReal>();
+            delta = v["delta"].cast<LongReal>();
             r_cut_squared = r_cut * r_cut;
             r_on_squared = r_on * r_on;
             }

From 91495a175cfef3b06c959509787660488b06f8cc Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 13 Jun 2024 21:26:59 +0000
Subject: [PATCH 27/32] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 hoomd/hpmc/PairPotentialExpandedGaussian.cc   | 21 ++++++++++--------
 .../pytest/test_pair_expanded_gaussian.py     | 22 +++++++++++--------
 2 files changed, 25 insertions(+), 18 deletions(-)

diff --git a/hoomd/hpmc/PairPotentialExpandedGaussian.cc b/hoomd/hpmc/PairPotentialExpandedGaussian.cc
index c0c8ead13f..9d8ee2063c 100644
--- a/hoomd/hpmc/PairPotentialExpandedGaussian.cc
+++ b/hoomd/hpmc/PairPotentialExpandedGaussian.cc
@@ -8,19 +8,20 @@ namespace hoomd
 namespace hpmc
     {
 
-PairPotentialExpandedGaussian::PairPotentialExpandedGaussian(std::shared_ptr<SystemDefinition> sysdef)
+PairPotentialExpandedGaussian::PairPotentialExpandedGaussian(
+    std::shared_ptr<SystemDefinition> sysdef)
     : PairPotential(sysdef), m_params(m_type_param_index.getNumElements())
     {
     }
 
 LongReal PairPotentialExpandedGaussian::energy(const LongReal r_squared,
-                                           const vec3<LongReal>& r_ij,
-                                           const unsigned int type_i,
-                                           const quat<LongReal>& q_i,
-                                           const LongReal charge_i,
-                                           const unsigned int type_j,
-                                           const quat<LongReal>& q_j,
-                                           const LongReal charge_j) const
+                                               const vec3<LongReal>& r_ij,
+                                               const unsigned int type_i,
+                                               const quat<LongReal>& q_i,
+                                               const LongReal charge_i,
+                                               const unsigned int type_j,
+                                               const quat<LongReal>& q_j,
+                                               const LongReal charge_j) const
     {
     unsigned int param_index = m_type_param_index(type_i, type_j);
     const auto& param = m_params[param_index];
@@ -82,7 +83,9 @@ void exportPairPotentialExpandedGaussian(pybind11::module& m)
     {
     pybind11::class_<PairPotentialExpandedGaussian,
                      PairPotential,
-                     std::shared_ptr<PairPotentialExpandedGaussian>>(m, "PairPotentialExpandedGaussian")
+                     std::shared_ptr<PairPotentialExpandedGaussian>>(
+        m,
+        "PairPotentialExpandedGaussian")
         .def(pybind11::init<std::shared_ptr<SystemDefinition>>())
         .def("setParams", &PairPotentialExpandedGaussian::setParamsPython)
         .def("getParams", &PairPotentialExpandedGaussian::getParamsPython)
diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index 052cf98f7b..a816a5e150 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -19,7 +19,8 @@
 def test_valid_construction(device, constructor_args):
     """Test that ExpandedGaussian can be constructed with valid arguments."""
     hoomd.hpmc.pair.ExpandedGaussian(**constructor_args)
-    
+
+
 @pytest.fixture(scope='session')
 def mc_simulation_factory(simulation_factory, two_particle_snapshot_factory):
     """Make a MC simulation with two particles separate dy by a distance d."""
@@ -36,13 +37,14 @@ def make_simulation(d=1):
 
     return make_simulation
 
+
 @pytest.mark.cpu
 def test_attaching(mc_simulation_factory):
     """Test that ExpandedGaussian attaches."""
     expanded_gauss = hoomd.hpmc.pair.ExpandedGaussian()
-    expanded_gauss.params[('A', 'A')] = dict(epsilon=1.0, 
-                                             sigma=1.0, 
-                                             delta=1.0, 
+    expanded_gauss.params[('A', 'A')] = dict(epsilon=1.0,
+                                             sigma=1.0,
+                                             delta=1.0,
                                              r_cut=2.5)
 
     simulation = mc_simulation_factory()
@@ -103,7 +105,7 @@ def xplor_factor(r, r_on, r_cut):
 
 def eg(r, epsilon, sigma, delta):
     """Compute the eg energy."""
-    return epsilon * math.exp(-0.5*(((r-delta)/sigma)**2))
+    return epsilon * math.exp(-0.5 * (((r - delta) / sigma)**2))
 
 
 # (pair params,
@@ -200,10 +202,12 @@ def test_multiple_pair_potentials(mc_simulation_factory):
     ]
     simulation.run(0)
 
-    assert expanded_gauss_1.energy == pytest.approx(
-        expected=eg(1.5, 1.0, 1.0, 1.0), rel=1e-5)
-    assert expanded_gauss_2.energy == pytest.approx(
-        expected=eg(1.5, 2.0, 1.0, 1.0), rel=1e-5)
+    assert expanded_gauss_1.energy == pytest.approx(expected=eg(
+        1.5, 1.0, 1.0, 1.0),
+                                                    rel=1e-5)
+    assert expanded_gauss_2.energy == pytest.approx(expected=eg(
+        1.5, 2.0, 1.0, 1.0),
+                                                    rel=1e-5)
     assert simulation.operations.integrator.pair_energy == pytest.approx(
         expected=eg(1.5, 1.0, 1.0, 1.0) + eg(1.5, 2.0, 1.0, 1.0), rel=1e-5)
 

From 05a0f84aa48eee9665b0854401331e73bc3aaf63 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Thu, 13 Jun 2024 17:32:26 -0400
Subject: [PATCH 28/32] Reformat docstrings to match line length restrictions

---
 hoomd/hpmc/pair/expanded_gaussian.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/hoomd/hpmc/pair/expanded_gaussian.py b/hoomd/hpmc/pair/expanded_gaussian.py
index 8edfa84672..0571806db1 100644
--- a/hoomd/hpmc/pair/expanded_gaussian.py
+++ b/hoomd/hpmc/pair/expanded_gaussian.py
@@ -26,10 +26,10 @@ class ExpandedGaussian(Pair):
           :math:`[\\mathrm{length}]`.
         mode (str): Energy shifting/smoothing mode.
 
-    `ExpandedGaussian` computes the Expanded Gaussian pair potential between every pair
-    of particles in the simulation state. The functional form of the potential,
-    including its behavior under shifting modes, is identical to that in
-    the MD pair potential `hoomd.md.pair.ExpandedGaussian`.
+    `ExpandedGaussian` computes the Expanded Gaussian pair potential between 
+    every pair of particles in the simulation state. The functional form of the 
+    potential, including its behavior under shifting modes, is identical to that
+    in the MD pair potential `hoomd.md.pair.ExpandedGaussian`.
 
     See Also:
         `hoomd.md.pair.ExpandedGaussian`
@@ -41,7 +41,10 @@ class ExpandedGaussian(Pair):
     .. code-block:: python
 
         expanded_gaussian =  hoomd.hpmc.pair.ExpandedGaussian()
-        expanded_gaussian.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0, delta=1.0, r_cut=2.5)
+        expanded_gaussian.params[('A', 'A')] = dict(epsilon=1.0, 
+                                                    sigma=1.0, 
+                                                    delta=1.0, 
+                                                    r_cut=2.5)
         simulation.operations.integrator.pair_potentials = [expanded_gaussian]
 
     .. py:attribute:: params

From f20215f2b07f978c01d1a99a1d6d01b5b343a773 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 13 Jun 2024 21:36:52 +0000
Subject: [PATCH 29/32] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 hoomd/hpmc/pair/expanded_gaussian.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/hoomd/hpmc/pair/expanded_gaussian.py b/hoomd/hpmc/pair/expanded_gaussian.py
index 0571806db1..3b886f994f 100644
--- a/hoomd/hpmc/pair/expanded_gaussian.py
+++ b/hoomd/hpmc/pair/expanded_gaussian.py
@@ -26,8 +26,8 @@ class ExpandedGaussian(Pair):
           :math:`[\\mathrm{length}]`.
         mode (str): Energy shifting/smoothing mode.
 
-    `ExpandedGaussian` computes the Expanded Gaussian pair potential between 
-    every pair of particles in the simulation state. The functional form of the 
+    `ExpandedGaussian` computes the Expanded Gaussian pair potential between
+    every pair of particles in the simulation state. The functional form of the
     potential, including its behavior under shifting modes, is identical to that
     in the MD pair potential `hoomd.md.pair.ExpandedGaussian`.
 
@@ -41,9 +41,9 @@ class ExpandedGaussian(Pair):
     .. code-block:: python
 
         expanded_gaussian =  hoomd.hpmc.pair.ExpandedGaussian()
-        expanded_gaussian.params[('A', 'A')] = dict(epsilon=1.0, 
-                                                    sigma=1.0, 
-                                                    delta=1.0, 
+        expanded_gaussian.params[('A', 'A')] = dict(epsilon=1.0,
+                                                    sigma=1.0,
+                                                    delta=1.0,
                                                     r_cut=2.5)
         simulation.operations.integrator.pair_potentials = [expanded_gaussian]
 

From 65645aa54770ba0c5c062a6bd879dc2a96966dbb Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Fri, 14 Jun 2024 10:11:09 -0400
Subject: [PATCH 30/32] Shorten interaction parameters in test to satisfy MPI
 requirements

---
 hoomd/hpmc/pytest/test_pair_expanded_gaussian.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index a816a5e150..5da084aed2 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -126,28 +126,28 @@ def eg(r, epsilon, sigma, delta):
         2.0,
     ),
     (
-        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=6.0),
+        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=2.0),
         'none',
         1.5,
         eg(1.5, 5, 1.1, 1.0),
     ),
     (
-        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=6.0),
+        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=2.0),
         'shift',
         1.5,
-        eg(1.5, 5, 1.1, 1.0) - eg(6.0, 5, 1.1, 1.0),
+        eg(1.5, 5, 1.1, 1.0) - eg(2.0, 5, 1.1, 1.0),
     ),
     (
-        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=6.0, r_on=2.0),
+        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=2.5, r_on=2.0),
         'xplor',
         1.5,
-        eg(1.5, 5.0, 1.1, 1.0) * xplor_factor(1.5, 2.0, 6.0),
+        eg(1.5, 5.0, 1.1, 1.0) * xplor_factor(1.5, 2.0, 2.5),
     ),
     (
-        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=6.0, r_on=2.0),
+        dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=2.5, r_on=2.0),
         'xplor',
         2.3,
-        eg(2.3, 5, 1.1, 1.0) * xplor_factor(2.3, 2.0, 6.0),
+        eg(2.3, 5, 1.1, 1.0) * xplor_factor(2.3, 2.0, 2.5),
     ),
     (
         dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=2.0, r_on=3),

From 3fed0c74f49a36a47d66aa166bbc9065d7b547b7 Mon Sep 17 00:00:00 2001
From: josephburkhart <joaburkh@umich.edu>
Date: Wed, 10 Jul 2024 15:49:07 -0400
Subject: [PATCH 31/32] Save expected results to variables

---
 hoomd/hpmc/pytest/test_pair_expanded_gaussian.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
index 5da084aed2..e3040e8b5d 100644
--- a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
+++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py
@@ -194,6 +194,8 @@ def test_multiple_pair_potentials(mc_simulation_factory):
                                                sigma=1.0,
                                                delta=1.0,
                                                r_cut=2.5)
+    expected_1 = eg(1.5, 1.0, 1.0, 1.0)
+    expected_2 = eg(1.5, 2.0, 1.0, 1.0)
 
     # Some parameters are validated only after attaching.
     simulation = mc_simulation_factory(1.5)
@@ -202,14 +204,12 @@ def test_multiple_pair_potentials(mc_simulation_factory):
     ]
     simulation.run(0)
 
-    assert expanded_gauss_1.energy == pytest.approx(expected=eg(
-        1.5, 1.0, 1.0, 1.0),
+    assert expanded_gauss_1.energy == pytest.approx(expected=expected_1,
                                                     rel=1e-5)
-    assert expanded_gauss_2.energy == pytest.approx(expected=eg(
-        1.5, 2.0, 1.0, 1.0),
+    assert expanded_gauss_2.energy == pytest.approx(expected=expected_2,
                                                     rel=1e-5)
     assert simulation.operations.integrator.pair_energy == pytest.approx(
-        expected=eg(1.5, 1.0, 1.0, 1.0) + eg(1.5, 2.0, 1.0, 1.0), rel=1e-5)
+        expected=expected_1 + expected_2, rel=1e-5)
 
 
 def test_logging():

From bcb6d6830324d80e3f72a616f5a6d3f61d5750f2 Mon Sep 17 00:00:00 2001
From: "Joshua A. Anderson" <joaander@umich.edu>
Date: Thu, 11 Jul 2024 16:11:33 -0400
Subject: [PATCH 32/32] Ensure that sigma is positive.

---
 hoomd/hpmc/pair/expanded_gaussian.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/hoomd/hpmc/pair/expanded_gaussian.py b/hoomd/hpmc/pair/expanded_gaussian.py
index 3b886f994f..8f611259f7 100644
--- a/hoomd/hpmc/pair/expanded_gaussian.py
+++ b/hoomd/hpmc/pair/expanded_gaussian.py
@@ -14,6 +14,7 @@
 import hoomd
 
 from .pair import Pair
+from hoomd.data.typeconverter import positive_real
 
 
 @hoomd.logging.modify_namespace(('hpmc', 'pair', 'ExpandedGaussian'))
@@ -92,7 +93,7 @@ def __init__(self, default_r_cut=None, default_r_on=0.0, mode='none'):
             'params', 'particle_types',
             hoomd.data.parameterdicts.TypeParameterDict(
                 epsilon=float,
-                sigma=float,
+                sigma=positive_real,
                 delta=float,
                 r_cut=default_r_cut,
                 r_on=float(default_r_on),