Pair hpmc angular step potential

hoomd/hpmc/CMakeLists.txt

@@ -41,6 +41,7 @@ set(_hpmc_sources   module.cc
                     PairPotentialLennardJones.cc
                     PairPotentialStep.cc
                     PairPotentialUnion.cc
+                    PairPotentialAngularStep.cc
                     ShapeUtils.cc
                     UpdaterBoxMC.cc
                     UpdaterQuickCompress.cc
@@ -82,6 +83,7 @@ set(_hpmc_headers
     PairPotentialLennardJones.h
     PairPotentialStep.h
     PairPotentialUnion.h
+    PairPotentialAngularStep.h
     ShapeConvexPolygon.h
     ShapeConvexPolyhedron.h
     ShapeEllipsoid.h

hoomd/hpmc/PairPotentialAngularStep.cc

@@ -0,0 +1,130 @@
+// Copyright (c) 2009-2024 The Regents of the University of Michigan.
+// Part of HOOMD-blue, released under the BSD 3-Clause License.
+
+#include "PairPotentialAngularStep.h"
+
+namespace hoomd
+    {
+namespace hpmc
+    {
+
+PairPotentialAngularStep::PairPotentialAngularStep(
+    std::shared_ptr<SystemDefinition> sysdef,
+    std::shared_ptr<PairPotential> isotropic_potential)
+    : PairPotential(sysdef), m_isotropic_potential(isotropic_potential)
+    {
+    unsigned int ntypes = m_sysdef->getParticleData()->getNTypes();
+    m_directors.resize(ntypes);
+    m_cos_deltas.resize(ntypes);
+
+    if (!m_isotropic_potential)
+        {
+        throw std::runtime_error("no isotropic potential given.");
+        }
+    }
+
+void PairPotentialAngularStep::setMask(std::string particle_type, pybind11::object v)
+    {
+    unsigned int particle_type_id = m_sysdef->getParticleData()->getTypeByName(particle_type);
+
+    if (v.is_none())
+        {
+        m_directors[particle_type_id].clear();
+        m_cos_deltas[particle_type_id].clear();
+        return;
+        }
+    pybind11::list directors = v["directors"];
+    pybind11::list deltas = v["deltas"];
+
+    auto N = pybind11::len(directors);
+
+    if (pybind11::len(deltas) != N)
+        {
+        throw std::runtime_error("the length of the delta list should match the length"
+                                 "of the director list.");
+        }
+
+    m_directors[particle_type_id].resize(N);
+    m_cos_deltas[particle_type_id].resize(N);
+
+    for (unsigned int i = 0; i < N; i++)
+        {
+        pybind11::tuple director_python = directors[i];
+        if (pybind11::len(director_python) != 3)
+            {
+            throw std::length_error("director must be a list of 3-tuples.");
+            }
+
+        m_directors[particle_type_id][i] = vec3<LongReal>(director_python[0].cast<LongReal>(),
+                                                          director_python[1].cast<LongReal>(),
+                                                          director_python[2].cast<LongReal>());
+
+        // normalize the directional vector
+        m_directors[particle_type_id][i] /= fast::sqrt(dot(m_directors[particle_type_id][i], 
+        m_directors[particle_type_id][i]));
+
+        pybind11::handle delta_python = deltas[i];
+        m_cos_deltas[particle_type_id][i] = cos(delta_python.cast<LongReal>());
+        }
+    }
+
+pybind11::object PairPotentialAngularStep::getMask(std::string particle_type)
+    {
+    unsigned int particle_type_id = m_sysdef->getParticleData()->getTypeByName(particle_type);
+    size_t N = m_directors[particle_type_id].size();
+
+    if (N == 0)
+        {
+        return pybind11::none();
+        }
+
+    pybind11::list directors;
+    pybind11::list deltas;
+
+    for (unsigned int i = 0; i < N; i++)
+        {
+        directors.append(pybind11::make_tuple(m_directors[particle_type_id][i].x,
+                                              m_directors[particle_type_id][i].y,
+                                              m_directors[particle_type_id][i].z));
+        deltas.append(acos(m_cos_deltas[particle_type_id][i]));
+        }
+
+    pybind11::dict v;
+    v["directors"] = directors;
+    v["deltas"] = deltas;
+    return v;
+    }
+
+// protected
+LongReal PairPotentialAngularStep::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
+    {
+    if (maskingFunction(r_squared, r_ij, type_i, q_i, type_j, q_j))
+        {
+        return m_isotropic_potential
+            ->energy(r_squared, r_ij, type_i, q_i, charge_i, type_j, q_j, charge_j);
+        }
+    return 0;
+    }
+
+namespace detail
+    {
+void exportPairPotentialAngularStep(pybind11::module& m)
+    {
+    pybind11::class_<PairPotentialAngularStep,
+                     PairPotential,
+                     std::shared_ptr<PairPotentialAngularStep>>(m, "PairPotentialAngularStep")
+        .def(pybind11::init<std::shared_ptr<SystemDefinition>, std::shared_ptr<PairPotential>>())
+        .def("setMask", &PairPotentialAngularStep::setMask)
+        .def("getMask", &PairPotentialAngularStep::getMask);
+    }
+    } // end namespace detail
+
+    } // end namespace hpmc
+    } // end namespace hoomd

hoomd/hpmc/PairPotentialAngularStep.h

@@ -0,0 +1,100 @@
+// 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 angular step potential between two particles.
+
+For use with HPMC simulations.
+*/
+class PairPotentialAngularStep : public hpmc::PairPotential
+    {
+    public:
+    PairPotentialAngularStep(std::shared_ptr<SystemDefinition> sysdef,
+                             std::shared_ptr<PairPotential> isotropic_potential);
+    virtual ~PairPotentialAngularStep() { }
+
+    /// Set the patch parameters for a given particle type.
+    void setMask(std::string particle_type, pybind11::object v);
+
+    /// Get the patch parameters for a given particle type.
+    pybind11::object getMask(std::string particle_type);
+
+    /*** Evaluate the energy of the pair interaction
+        @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.
+        @param charge_i Charge of particle i.
+        @param q_i Orientation quaternion of particle i.
+        @param type_j Integer type index of particle j.
+        @param q_j Orientation quaternion of particle j.
+        @param charge_j Charge of particle j.
+        @returns Energy of the pair interaction.
+    */
+    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;
+
+    virtual LongReal computeRCutNonAdditive(unsigned int type_i, unsigned int type_j) const
+        {
+        // Pass on the non-additive r_cut from the isotropic (parent) potential.
+        return m_isotropic_potential->computeRCutNonAdditive(type_i, type_j);
+        }
+
+    virtual LongReal computeRCutAdditive(unsigned int type) const
+        {
+        // Returns the additive part of the cutoff distance for a given type.
+        return m_isotropic_potential->computeRCutAdditive(type);
+        }
+
+    protected:
+    bool maskingFunction(LongReal r_squared,
+                         const vec3<LongReal>& r_ij,
+                         const unsigned int type_i,
+                         const quat<LongReal>& q_i,
+                         const unsigned int type_j,
+                         const quat<LongReal>& q_j) const
+        {
+        vec3<LongReal> rhat_ij = r_ij / fast::sqrt(r_squared);
+
+        for (size_t m = 0; m < m_directors[type_i].size(); m++)
+            {
+            LongReal cos_delta_m = m_cos_deltas[type_i][m];
+            vec3<LongReal> ehat_m = rotate(q_i, m_directors[type_i][m]);
+
+            for (size_t n = 0; n < m_directors[type_j].size(); n++)
+                {
+                LongReal cos_delta_n = m_cos_deltas[type_j][n];
+                vec3<LongReal> ehat_n = rotate(q_j, m_directors[type_j][n]);
+
+                if (dot(ehat_m, rhat_ij) >= cos_delta_m && dot(ehat_n, -rhat_ij) >= cos_delta_n)
+                    {
+                    return true;
+                    }
+                }
+            }
+        return false;
+        }
+
+    /// Create a pointer to get the isotropic pair potential
+    std::shared_ptr<PairPotential> m_isotropic_potential;
+
+    /// Type pair parameters of potential
+    std::vector<std::vector<vec3<LongReal>>> m_directors;
+    std::vector<std::vector<LongReal>> m_cos_deltas;
+    };
+
+    } // end namespace hpmc
+    } // end namespace hoomd

hoomd/hpmc/module.cc

@@ -46,6 +46,8 @@ void exportPairPotential(pybind11::module& m);
 
 void exportPairPotentialLennardJones(pybind11::module& m);
 
+void exportPairPotentialAngularStep(pybind11::module& m);
+
 void exportPairPotentialStep(pybind11::module& m);
 
 void exportPairPotentialUnion(pybind11::module& m);
@@ -146,6 +148,7 @@ PYBIND11_MODULE(_hpmc, m)
 
     exportPairPotential(m);
     exportPairPotentialLennardJones(m);
+    exportPairPotentialAngularStep(m);
     exportPairPotentialStep(m);
     exportPairPotentialUnion(m);
     }

hoomd/hpmc/pair/CMakeLists.txt

@@ -3,6 +3,7 @@ set(files __init__.py
         pair.py
         step.py
         union.py
+        angular_step.py
         user.py
  )
 

hoomd/hpmc/pair/__init__.py

@@ -31,4 +31,5 @@
 from .pair import Pair
 from .lennard_jones import LennardJones
 from .union import Union
+from .angular_step import AngularStep
 from .step import Step