diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ffdf67c1..aa2a2fde 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -22,7 +22,10 @@ set(python_files
set(_bond_evaluators )
# TODO: Add names of all pair evaluators`
-set(_pair_evaluators Hertz)
+set(_pair_evaluators
+ Hertz
+ PerturbedLennardJones
+ )
# TODO: Add names of all anisotropic pair evaluators
set(_aniso_pair_evaluators )
diff --git a/src/PairEvaluatorAshbaugh.h b/src/PairEvaluatorAshbaugh.h
deleted file mode 100644
index e4745875..00000000
--- a/src/PairEvaluatorAshbaugh.h
+++ /dev/null
@@ -1,164 +0,0 @@
-// Copyright (c) 2018-2020, Michael P. Howard
-// Copyright (c) 2021-2024, Auburn University
-// Part of azplugins, released under the BSD 3-Clause License.
-
-/*!
- * \file PairEvaluatorAshbaugh.h
- * \brief Defines the pair force evaluator class for Ashbaugh-Hatch potential
- */
-
-#ifndef AZPLUGINS_PAIR_EVALUATOR_ASHBAUGH_H_
-#define AZPLUGINS_PAIR_EVALUATOR_ASHBAUGH_H_
-
-#include "PairEvaluator.h"
-
-#ifdef NVCC
-#define DEVICE __device__
-#define HOSTDEVICE __host__ __device__
-#else
-#define DEVICE
-#define HOSTDEVICE
-#endif
-
-namespace azplugins
- {
-
-namespace detail
- {
-//! Ashbaugh-Hatch parameters
-/*!
- * \sa PairEvaluatorAshbaugh
- */
-struct ashbaugh_params
- {
- Scalar lj1; //!< The coefficient for 1/r^12
- Scalar lj2; //!< The coefficient for 1/r^6
- Scalar lambda; //!< Controls the attractive tail, between 0 and 1
- Scalar rwcasq; //!< The square of the location of the LJ potential minimum
- Scalar wca_shift; //!< The amount to shift the repulsive part by
- };
-
-//! Convenience function for making ashbaugh_params in python
-HOSTDEVICE inline ashbaugh_params
-make_ashbaugh_params(Scalar lj1, Scalar lj2, Scalar lambda, Scalar rwcasq, Scalar wca_shift)
- {
- ashbaugh_params p;
- p.lj1 = lj1;
- p.lj2 = lj2;
- p.lambda = lambda;
- p.rwcasq = rwcasq;
- p.wca_shift = wca_shift;
- return p;
- }
-
-//! Class for evaluating the Ashbaugh-Hatch pair potential
-/*!
- * PairEvaluatorAshbaugh evaluates the function:
- * \f{eqnarray*}{
- * V(r) = & V_{\mathrm{LJ}}(r, \varepsilon, \sigma, \alpha) + (1-\lambda)\varepsilon & r <
- * (2/\alpha)^{1/6}\sigma \\
- * = & \lambda V_{\mathrm{LJ}}(r, \varepsilon, \sigma, \alpha) & (2/\alpha)^{1/6}\sigma
- * \ge r < r_{\mathrm{cut}} \\ = & 0 & r \ge r_{\mathrm{cut}} \f}
- *
- * where \f$V_{\mathrm{LJ}}(r,\varepsilon,\sigma,\alpha)\f$ is the standard Lennard-Jones potential
- * (see EvaluatorPairLJ) with parameters \f$\varepsilon\f$, \f$\sigma\f$, and \f$\alpha\f$
- * (default: 1.0). This potential is implemented as given in H.S. Ashbaugh and H.W. Hatch, J. Am. Chem. Soc., 130,
- * 9536 (2008).
- *
- * The Ashbaugh-Hatch potential does not need diameter or charge. Five parameters are specified and
- * stored in a ashbaugh_params. These are related to the standard Lennard-Jones and Ashbaugh-Hatch
- * parameters by:
- * - \a lj1 = 4.0 * epsilon * pow(sigma,12.0)
- * - \a lj2 = alpha * 4.0 * epsilon * pow(sigma,6.0);
- * - \a lambda is the scale factor for the attraction (0 = WCA, 1 = LJ)
- * - \a rwcasq is the square of the location of the potential minimum (WCA cutoff),
- * pow(2.0/alpha,1./3.) * sigma * sigma
- * - \a wca_shift is the amount needed to shift the energy of the repulsive part to match the
- * attractive energy.
- */
-class PairEvaluatorAshbaugh : public PairEvaluator
- {
- public:
- //! Define the parameter type used by this pair potential evaluator
- typedef ashbaugh_params param_type;
-
- //! Constructor
- /*!
- * \param _rsq Squared distance between particles
- * \param _rcutsq Cutoff radius squared
- * \param _params Pair potential parameters, given by typedef above
- *
- * The functor initializes its members from \a _params.
- */
- DEVICE PairEvaluatorAshbaugh(Scalar _rsq, Scalar _rcutsq, const param_type& _params)
- : PairEvaluator(_rsq, _rcutsq), lj1(_params.lj1), lj2(_params.lj2), lambda(_params.lambda),
- rwcasq(_params.rwcasq), wca_shift(_params.wca_shift)
- {
- }
-
- //! Evaluate the force and energy
- /*!
- * \param force_divr Holds the computed force divided by r
- * \param pair_eng Holds the computed pair energy
- * \param energy_shift If true, the potential is shifted to zero at the cutoff
- *
- * \returns True if the energy calculation occurs
- *
- * The calculation does not occur if the pair distance is greater than the cutoff
- * or if the potential is scaled to zero.
- */
- DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& pair_eng, bool energy_shift)
- {
- if (rsq < rcutsq && lj1 != 0)
- {
- Scalar r2inv = Scalar(1.0) / rsq;
- Scalar r6inv = r2inv * r2inv * r2inv;
- force_divr = r2inv * r6inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2);
-
- pair_eng = r6inv * (lj1 * r6inv - lj2);
- if (rsq < rwcasq)
- {
- pair_eng += wca_shift;
- }
- else
- {
- force_divr *= lambda;
- pair_eng *= lambda;
- }
-
- if (energy_shift)
- {
- Scalar rcut2inv = Scalar(1.0) / rcutsq;
- Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv;
- pair_eng -= lambda * rcut6inv * (lj1 * rcut6inv - lj2);
- }
- return true;
- }
- else
- return false;
- }
-
-#ifndef NVCC
- //! Return the name of this potential
- static std::string getName()
- {
- return std::string("ashbaugh");
- }
-#endif
-
- protected:
- Scalar lj1; //!< lj1 parameter extracted from the params passed to the constructor
- Scalar lj2; //!< lj2 parameter extracted from the params passed to the constructor
- Scalar lambda; //!< lambda parameter
- Scalar rwcasq; //!< WCA cutoff radius squared
- Scalar wca_shift; //!< Energy shift for WCA part of the potential
- };
-
- } // end namespace detail
- } // end namespace azplugins
-
-#undef DEVICE
-#undef HOSTDEVICE
-
-#endif // AZPLUGINS_PAIR_EVALUATOR_ASHBAUGH_H_
diff --git a/src/PairEvaluatorPerturbedLennardJones.h b/src/PairEvaluatorPerturbedLennardJones.h
new file mode 100644
index 00000000..b9072896
--- /dev/null
+++ b/src/PairEvaluatorPerturbedLennardJones.h
@@ -0,0 +1,180 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2024, Auburn University
+// Part of azplugins, released under the BSD 3-Clause License.
+
+#ifndef AZPLUGINS_PAIR_EVALUATOR_PERTURBED_LENNARD_JONES_H_
+#define AZPLUGINS_PAIR_EVALUATOR_PERTURBED_LENNARD_JONES_H_
+
+#include "PairEvaluator.h"
+
+#ifdef __HIPCC__
+#define DEVICE __device__
+#define HOSTDEVICE __host__ __device__
+#else
+#define DEVICE
+#define HOSTDEVICE
+#endif
+
+namespace hoomd
+ {
+namespace azplugins
+ {
+namespace detail
+ {
+//! Define the parameter type used by this pair potential evaluator
+struct PairParametersPerturbedLennardJones : public PairParameters
+ {
+#ifndef __HIPCC__
+ PairParametersPerturbedLennardJones()
+ : sigma_6(0), epsilon_x_4(0), attraction_scale_factor(0), rwcasq(0)
+ {
+ }
+
+ PairParametersPerturbedLennardJones(pybind11::dict v, bool managed = false)
+ {
+ auto sigma(v["sigma"].cast());
+ auto epsilon(v["epsilon"].cast());
+
+ const Scalar sigma_2 = sigma * sigma;
+ const Scalar sigma_4 = sigma_2 * sigma_2;
+ sigma_6 = sigma_2 * sigma_4;
+
+ epsilon_x_4 = Scalar(4.0) * epsilon;
+ attraction_scale_factor = v["attraction_scale_factor"].cast();
+ rwcasq = pow(Scalar(2.), 1. / 3.) * sigma_2;
+ }
+
+ pybind11::dict asDict()
+ {
+ pybind11::dict v;
+ v["sigma"] = pow(sigma_6, 1. / 6.);
+ v["epsilon"] = epsilon_x_4 / Scalar(4.);
+ v["attraction_scale_factor"] = attraction_scale_factor;
+ return v;
+ }
+#endif // __HIPCC__
+
+ Scalar sigma_6; //!< The coefficient for 1/r^12
+ Scalar epsilon_x_4; //!< The coefficient for 1/r^6
+ Scalar attraction_scale_factor; //!< Controls the attractive tail, between 0 and 1
+ Scalar rwcasq; //!< The square of the location of the LJ potential minimum
+ }
+#if HOOMD_LONGREAL_SIZE == 32
+ __attribute__((aligned(16)));
+#else
+ __attribute__((aligned(32)));
+#endif
+
+//! Class for evaluating the perturbed Lennard-Jones pair potential
+/*!
+ * This class evaluates the function:
+ * \f{eqnarray*}{
+ * V(r) = & V_{\mathrm{LJ}}(r, \varepsilon, \sigma)
+ * + (1-\lambda)\varepsilon & r < 2^{1/6}\sigma \\
+ * = & \lambda V_{\mathrm{LJ}}(r, \varepsilon, \sigma) &
+ * 2^{1/6}\sigma \ge r < r_{\mathrm{cut}} \\
+ * = & 0 & r \ge r_{\mathrm{cut}}
+ * \f}
+ *
+ * where \f$V_{\mathrm{LJ}}(r,\varepsilon,\sigma)\f$ is the standard
+ * Lennard-Jones potential (see EvaluatorPairLJ) with parameters
+ * \f$\varepsilon\f$ and \f$\sigma\f$.
+ */
+class PairEvaluatorPerturbedLennardJones : public PairEvaluator
+ {
+ public:
+ typedef PairParametersPerturbedLennardJones param_type;
+
+ //! Constructor
+ /*!
+ * \param _rsq Squared distance between particles
+ * \param _rcutsq Cutoff radius squared
+ * \param _params Pair potential parameters, given by typedef above
+ *
+ * The functor initializes its members from \a _params.
+ */
+ DEVICE
+ PairEvaluatorPerturbedLennardJones(Scalar _rsq, Scalar _rcutsq, const param_type& _params)
+ : PairEvaluator(_rsq, _rcutsq),
+ lj1(_params.epsilon_x_4 * _params.sigma_6 * _params.sigma_6),
+ lj2(_params.epsilon_x_4 * _params.sigma_6),
+ attraction_scale_factor(_params.attraction_scale_factor), rwcasq(_params.rwcasq)
+ {
+ wca_shift = _params.epsilon_x_4 * (Scalar(1.0) - attraction_scale_factor) / Scalar(4.0);
+ }
+
+ //! Evaluate the force and energy
+ /*!
+ * \param force_divr Holds the computed force divided by r
+ * \param pair_eng Holds the computed pair energy
+ * \param energy_shift If true, the potential is shifted to zero at the cutoff
+ *
+ * \returns True if the energy calculation occurs
+ *
+ * The calculation does not occur if the pair distance is greater than the cutoff
+ * or if the potential is scaled to zero.
+ */
+ DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& pair_eng, bool energy_shift)
+ {
+ if (rsq < rcutsq && lj1 != 0)
+ {
+ const Scalar r2inv = Scalar(1.0) / rsq;
+ const Scalar r6inv = r2inv * r2inv * r2inv;
+ force_divr = r2inv * r6inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2);
+
+ pair_eng = r6inv * (lj1 * r6inv - lj2);
+ if (rsq < rwcasq)
+ {
+ pair_eng += wca_shift;
+ }
+ else
+ {
+ force_divr *= attraction_scale_factor;
+ pair_eng *= attraction_scale_factor;
+ }
+
+ if (energy_shift)
+ {
+ const Scalar rcut2inv = Scalar(1.0) / rcutsq;
+ const Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv;
+ Scalar pair_eng_shift = rcut6inv * (lj1 * rcut6inv - lj2);
+ if (rcutsq < rwcasq)
+ {
+ pair_eng_shift += wca_shift;
+ }
+ else
+ {
+ pair_eng_shift *= attraction_scale_factor;
+ }
+ pair_eng -= pair_eng_shift;
+ }
+ return true;
+ }
+ else
+ return false;
+ }
+
+#ifndef __HIPCC__
+ //! Return the name of this potential
+ static std::string getName()
+ {
+ return std::string("PerturbedLennardJones");
+ }
+#endif
+
+ private:
+ Scalar lj1; //!< lj1 parameter - 4 epsilon sigma^12
+ Scalar lj2; //!< lj2 parameter - 4 epsilon sigma^6
+ Scalar attraction_scale_factor; //!< attraction_scale_factor parameter
+ Scalar rwcasq; //!< WCA cutoff radius squared
+ Scalar wca_shift; //!< Energy shift for WCA part of the potential
+ };
+
+ } // end namespace detail
+ } // end namespace azplugins
+ } // end namespace hoomd
+
+#undef DEVICE
+#undef HOSTDEVICE
+
+#endif // AZPLUGINS_PAIR_EVALUATOR_PERTURBED_LENNARD_JONES_H_
diff --git a/src/module.cc b/src/module.cc
index 0f400bf8..09f17650 100644
--- a/src/module.cc
+++ b/src/module.cc
@@ -51,9 +51,11 @@ namespace detail
{
void export_PotentialPairHertz(pybind11::module&);
+void export_PotentialPairPerturbedLennardJones(pybind11::module&);
#ifdef ENABLE_HIP
void export_PotentialPairHertzGPU(pybind11::module&);
+void export_PotentialPairPerturbedLennardJonesGPU(pybind11::module&);
#endif // ENABLE_HIP
} // namespace detail
@@ -66,8 +68,10 @@ PYBIND11_MODULE(_azplugins, m)
using namespace hoomd::azplugins::detail;
export_PotentialPairHertz(m);
+ export_PotentialPairPerturbedLennardJones(m);
#ifdef ENABLE_HIP
export_PotentialPairHertzGPU(m);
+ export_PotentialPairPerturbedLennardJonesGPU(m);
#endif // ENABLE_HIP
}
diff --git a/src/pair.py b/src/pair.py
index eb5f74b8..71f90d8a 100644
--- a/src/pair.py
+++ b/src/pair.py
@@ -66,3 +66,80 @@ def __init__(self, nlist, default_r_cut=None, default_r_on=0, mode='none'):
TypeParameterDict(epsilon=float, len_keys=2),
)
self._add_typeparam(params)
+
+
+class PerturbedLennardJones(pair.Pair):
+ r"""Perturbed Lennard-Jones potential.
+
+ Args:
+ nlist (hoomd.md.nlist.NeighborList): Neighbor list.
+ default_r_cut (float): Default cutoff radius :math:`[\mathrm{length}]`.
+ default_r_on (float): Default turn-on radius :math:`[\mathrm{length}]`.
+ mode (str): Energy shifting/smoothing mode.
+
+ :py:class:`PerturbedLennardJones` is a Lennard-Jones perturbation potential.
+ The potential has a purely repulsive (Weeks-Chandler-Andersen) core, with a
+ parameter :math:`attraction_scale_factor` (\lambda) setting the strength of
+ the attractive tail. When
+ :math:`attraction_scale_factor` is 0, the potential is purely repulsive. When
+ :math:`attraction_scale_factor` is 1, the potential is the standard
+ Lennard-Jones potential (see :py:class:`hoomd.md.pair.LJ` for details).
+
+ .. math::
+ :nowrap:
+
+ \begin{eqnarray*}
+ V(r) = & V_{\mathrm{LJ}}(r, \varepsilon, \sigma) +
+ (1-\lambda)\varepsilon & r < 2^{1/6}\sigma \\
+ = & \lambda V_{\mathrm{LJ}}(
+ r, \varepsilon, \sigma) & 2^{1/6}\sigma \ge
+ r < r_{\mathrm{cut}} \\
+ = & 0 & r \ge r_{\mathrm{cut}}
+ \end{eqnarray*}
+
+
+ Example::
+
+ nl = nlist.Cell()
+ perturbed_lj = pair.PerturbedLennardJones(default_r_cut=3.0, nlist=nl)
+ perturbed_lj.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0,
+ attraction_scale_factor=0.5)
+ perturbed_lj.r_cut[('A', 'B')] = 3.0
+
+ .. py:attribute:: params
+
+ The Perturbed LJ potential parameters. The dictionary has the following
+ keys:
+
+ * ``epsilon`` (`float`, **required**) - energy parameter
+ :math:`\varepsilon` :math:`[\mathrm{energy}]`
+ * ``sigma`` (`float`, **required**) - particle size :math:`\sigma`
+ :math:`[\mathrm{length}]`
+ * ``attraction_scale_factor`` (`float`, **required**) - scale factor
+ for attraction, between 0 and 1 :math:`attraction_scale_factor`
+ :math:`[\mathrm{dimensionless}]`
+
+ Type: `TypeParameter` [`tuple` [``particle_type``, ``particle_type``],
+ `dict`]
+
+ .. py:attribute:: mode
+
+ Energy shifting/smoothing mode: ``"none"``, ``"shift"``, or ``"xplor"``.
+
+ Type: `str`
+ """
+
+ _ext_module = _azplugins
+ _cpp_class_name = 'PotentialPairPerturbedLennardJones'
+ _accepted_modes = ('none', 'shift', 'xplor')
+
+ def __init__(self, nlist, default_r_cut=None, default_r_on=0, mode='none'):
+ super().__init__(nlist, default_r_cut, default_r_on, mode)
+ params = TypeParameter(
+ 'params',
+ 'particle_types',
+ TypeParameterDict(
+ epsilon=float, sigma=float, attraction_scale_factor=float, len_keys=2
+ ),
+ )
+ self._add_typeparam(params)
diff --git a/src/pytest/test_pair.py b/src/pytest/test_pair.py
index 0a58e4b9..32c11a29 100644
--- a/src/pytest/test_pair.py
+++ b/src/pytest/test_pair.py
@@ -17,7 +17,6 @@
)
potential_tests = []
-
# Hertz
potential_tests += [
# test the calculation of force and potential
@@ -61,6 +60,95 @@
0,
),
]
+# PerturbedLennardJones
+potential_tests += [
+ # test the calculation of force and potential
+ # test when it's in the wca part, no potential shifting
+ PotentialTestCase(
+ hoomd.azplugins.pair.PerturbedLennardJones,
+ {'epsilon': 2.0, 'sigma': 1.05, 'attraction_scale_factor': 0.0},
+ 3.0,
+ False,
+ 1.05,
+ 2.0,
+ 45.7143,
+ ),
+ # change attraction_scale_factor to check for shifting
+ # of energy (force stays the same)
+ PotentialTestCase(
+ hoomd.azplugins.pair.PerturbedLennardJones,
+ {'epsilon': 2.0, 'sigma': 1.05, 'attraction_scale_factor': 0.5},
+ 3.0,
+ False,
+ 1.05,
+ 1.0,
+ 45.7143,
+ ),
+ # change sigma so that now the particle is in the LJ region
+ # when attraction_scale_factor = 0, then the potential and force are zero
+ PotentialTestCase(
+ hoomd.azplugins.pair.PerturbedLennardJones,
+ {'epsilon': 2.0, 'sigma': 0.5, 'attraction_scale_factor': 0.0},
+ 3.0,
+ False,
+ 1.05,
+ 0,
+ 0,
+ ),
+ # partially switch on the LJ with attraction_scale_factor = 0.5
+ PotentialTestCase(
+ hoomd.azplugins.pair.PerturbedLennardJones,
+ {'epsilon': 2.0, 'sigma': 0.5, 'attraction_scale_factor': 0.5},
+ 3.0,
+ False,
+ 1.05,
+ -0.0460947,
+ -0.260291,
+ ),
+ # test that energy shifting works (bump up sigma so that at
+ # rcut = 3 the shift is reasonable)
+ # check wca is shifted first
+ PotentialTestCase(
+ hoomd.azplugins.pair.PerturbedLennardJones,
+ {'epsilon': 2.0, 'sigma': 1.05, 'attraction_scale_factor': 0.5},
+ 3.0,
+ True,
+ 1.05,
+ 1.00734,
+ 45.7143,
+ ),
+ # and check lj
+ PotentialTestCase(
+ hoomd.azplugins.pair.PerturbedLennardJones,
+ {'epsilon': 2.0, 'sigma': 0.85, 'attraction_scale_factor': 0.5},
+ 3.0,
+ True,
+ 1.05,
+ -0.806849,
+ -2.81197,
+ ),
+ # test the cases where the potential should be zero
+ # outside cutoff
+ PotentialTestCase(
+ hoomd.azplugins.pair.PerturbedLennardJones,
+ {'epsilon': 1.0, 'sigma': 1.0, 'attraction_scale_factor': 0.5},
+ 1.0,
+ False,
+ 1.05,
+ 0,
+ 0,
+ ),
+ # inside cutoff but epsilon = 0
+ PotentialTestCase(
+ hoomd.azplugins.pair.PerturbedLennardJones,
+ {'epsilon': 0.0, 'sigma': 1.0, 'attraction_scale_factor': 0.5},
+ 3.0,
+ False,
+ 1.05,
+ 0,
+ 0,
+ ),
+]
@pytest.mark.parametrize(