Union Potentials in HPMC

hoomd/hpmc/CMakeLists.txt

@@ -39,6 +39,7 @@ set(_hpmc_sources   module.cc
                     ExternalFieldWall.cc
                     PairPotential.cc
                     PairPotentialLennardJones.cc
+                    PairPotentialUnion.cc
                     ShapeUtils.cc
                     UpdaterBoxMC.cc
                     UpdaterQuickCompress.cc
@@ -78,6 +79,7 @@ set(_hpmc_headers
     OBBTree.h
     PairPotential.h
     PairPotentialLennardJones.h
+    PairPotentialUnion.h
     ShapeConvexPolygon.h
     ShapeConvexPolyhedron.h
     ShapeEllipsoid.h

hoomd/hpmc/PairPotential.cc

@@ -11,7 +11,8 @@ namespace hoomd::hpmc::detail
 void exportPairPotential(pybind11::module& m)
     {
     pybind11::class_<hpmc::PairPotential, std::shared_ptr<hpmc::PairPotential>>(m, "PairPotential")
-        .def(pybind11::init<std::shared_ptr<SystemDefinition>>());
+        .def(pybind11::init<std::shared_ptr<SystemDefinition>>())
+        .def("setParent", &hpmc::PairPotential::setParent);
 
     pybind11::bind_vector<std::vector<std::shared_ptr<PairPotential>>>(m, "PairPotentialList");
     }

hoomd/hpmc/PairPotential.h

@@ -33,8 +33,19 @@ namespace hpmc
 
     Of the two, the additive r_cut provides better performance when particle types have different
     sizes. The non-additive r_cut allows more freedom to implement models. Subclasses can use
-    one of the two or both at the same time. r_cut_additive and r_cut_non_additive default to all
-    0's - subclasses must call setRCutNonAdditive and/or setRCutAdditive to update the arrays.
+    one of the two or both at the same time.
+
+    A PairPotential may be compose another PairPotential (e.g. angular patches modulating a square
+    well). Each PairPotential maintains a weak pointer to its parent. Only the r_cut of the top most
+    parent is considered by IntegratorHPMC when computing interactions. Each PairPotential subclass
+    should override the default computeRCutNonAdditive and computeRCutAdditive methods as needed
+    to compute the various r_cut values as a function of parameters and child potentials. The base
+    implementations return 0.
+
+    The top level potential maintains cached values of the total and maximum r_cut values.
+    Subclasses must call notifyRCutChanged whenever they would change the value of their computed
+    r_cut values. The cached values ensure that we can use non-virtual inlined calls in the
+    inner loops where the total r_cut values are checked.
 */
 class PairPotential
     {
@@ -67,16 +78,31 @@ class PairPotential
         }
 
     /// Returns the maximum non-additive r_cut beteween any two types.
-    inline LongReal getMaxRCutNonAdditive()
+    inline LongReal getMaxRCutNonAdditive() const
         {
         return m_max_r_cut_non_additive;
         }
 
+    /// Set the parent potential
+    void setParent(std::shared_ptr<PairPotential> parent)
+        {
+        m_parent = parent;
+
+        if (auto parent = m_parent.lock())
+            {
+            parent->notifyRCutChanged();
+            }
+        }
+
     /*** Evaluate the energy of the pair interaction
 
         energy is given a pre-computed r_squared as many potentials use this parameter and the
         caller has already computed it.
 
+        To avoid repeated evaluations of r_square < r_cut_squared, the *caller* must perform the
+        check before calling energy. Implementations of energy are free to compute non-zero values
+        beyond r_cut when convenient.
+
         @param r_squared Pre-computed dot(r_ij, r_ij).
         @param r_ij Vector pointing from particle i to j.
         @param type_i Integer type index of particle i.
@@ -99,32 +125,34 @@ class PairPotential
         return 0;
         }
 
+    /// Compute the non-additive cuttoff radius
+    virtual LongReal computeRCutNonAdditive(unsigned int type_i, unsigned int type_j) const
+        {
+        return 0;
+        }
+
+    /// Returns the additive part of the cutoff distance for a given type.
+    virtual LongReal computeRCutAdditive(unsigned int type) const
+        {
+        return m_r_cut_additive[type];
+        }
+
     protected:
     /// The system definition.
     std::shared_ptr<SystemDefinition> m_sysdef;
 
     /// Indexer to access arrays by pairs of type parameters
     Index2D m_type_param_index;
 
-    /// Set the non-additive cutoff radius for both i,j and j,i pairs.
-    void setRCutNonAdditive(unsigned int type_i, unsigned int type_j, LongReal r)
+    /// Notify all parents that r_cut has changed.
+    void notifyRCutChanged()
         {
-        m_r_cut_non_additive[m_type_param_index(type_i, type_j)] = r;
-        m_r_cut_non_additive[m_type_param_index(type_j, type_i)] = r;
-
-        updateRCutTotal(type_i, type_j);
-        updateMaxRCutNonAdditive();
-        }
-
-    /// Set the additive part of the cutoff distance for a given type.
-    void setRCutAdditive(unsigned int type, LongReal r)
-        {
-        m_r_cut_additive[type] = r;
-
-        for (unsigned int type_j = 0; type_j < m_sysdef->getParticleData()->getNTypes(); type_j++)
+        if (auto parent = m_parent.lock())
             {
-            updateRCutTotal(type, type_j);
+            parent->notifyRCutChanged();
             }
+
+        updateRCutCache();
         }
 
     private:
@@ -140,31 +168,38 @@ class PairPotential
     /// Pre-computed maximum additive r_cut.
     LongReal m_max_r_cut_non_additive = 0;
 
-    /// Update the pre-computed r_cut_total value.
-    void updateRCutTotal(unsigned int type_i, unsigned int type_j)
+    /// Parent potential
+    std::weak_ptr<PairPotential> m_parent;
+
+    /// Update r_cut cache
+    void updateRCutCache()
         {
-        unsigned int type_pair = m_type_param_index(type_i, type_j);
-        LongReal r_cut = m_r_cut_non_additive[type_pair]
-                         + LongReal(0.5) * (m_r_cut_additive[type_i] + m_r_cut_additive[type_j]);
-        m_r_cut_squared_total[type_pair] = r_cut * r_cut;
+        const unsigned int n_types = m_sysdef->getParticleData()->getNTypes();
 
-        type_pair = m_type_param_index(type_j, type_i);
-        m_r_cut_squared_total[type_pair] = r_cut * r_cut;
-        }
+        for (unsigned int type = 0; type < n_types; type++)
+            {
+            m_r_cut_additive[type] = computeRCutAdditive(type);
+            }
 
-    /// Update the pre-computed m_max_r_cut_non_additive value.
-    void updateMaxRCutNonAdditive()
-        {
         m_max_r_cut_non_additive = 0;
-        const unsigned int n_types = m_sysdef->getParticleData()->getNTypes();
 
         for (unsigned int type_i = 0; type_i < n_types; type_i++)
             {
             for (unsigned int type_j = type_i; type_j < n_types; type_j++)
                 {
-                unsigned int type_pair = m_type_param_index(type_i, type_j);
-                LongReal r_cut_non_additive = m_r_cut_non_additive[type_pair];
-                m_max_r_cut_non_additive = std::max(m_max_r_cut_non_additive, r_cut_non_additive);
+                LongReal r_cut_ij = computeRCutNonAdditive(type_i, type_j);
+                m_max_r_cut_non_additive = std::max(m_max_r_cut_non_additive, r_cut_ij);
+
+                unsigned int param_index_1 = m_type_param_index(type_i, type_j);
+                unsigned int param_index_2 = m_type_param_index(type_j, type_i);
+                m_r_cut_non_additive[param_index_1] = r_cut_ij;
+                m_r_cut_non_additive[param_index_2] = r_cut_ij;
+
+                LongReal r_cut
+                    = m_r_cut_non_additive[param_index_1]
+                      + LongReal(0.5) * (m_r_cut_additive[type_i] + m_r_cut_additive[type_j]);
+                m_r_cut_squared_total[param_index_1] = r_cut * r_cut;
+                m_r_cut_squared_total[param_index_2] = r_cut * r_cut;
                 }
             }
         }

hoomd/hpmc/PairPotentialLennardJones.cc

@@ -65,7 +65,7 @@ void PairPotentialLennardJones::setParamsPython(pybind11::tuple typ, pybind11::d
     unsigned int param_index_2 = m_type_param_index(type_j, type_i);
     m_params[param_index_2] = ParamType(params);
 
-    setRCutNonAdditive(type_i, type_j, params["r_cut"].cast<LongReal>());
+    notifyRCutChanged();
     }
 
 pybind11::dict PairPotentialLennardJones::getParamsPython(pybind11::tuple typ)

hoomd/hpmc/PairPotentialLennardJones.h

@@ -3,7 +3,7 @@
 
 #pragma once
 
-#include "hoomd/hpmc/IntegratorHPMC.h"
+#include "PairPotential.h"
 
 namespace hoomd
     {
@@ -29,6 +29,13 @@ class PairPotentialLennardJones : public hpmc::PairPotential
                             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);