Merge branch 'refactor/flow' into refactor/integrator

Author: mayukh33

.ruff.toml

@@ -32,6 +32,8 @@ ignore = [
 
 "__init__.py" = ["F401", # __init__.py import submodules for use by the package importer.
 ]
+"*/pytest/*.py" = ["PLR2004", # allow "magic numbers" in tests
+]
 
 [lint.pydocstyle]
 convention = "google"

src/BondEvaluatorDoubleWell.h

@@ -2,151 +2,133 @@
 // Copyright (c) 2021-2024, Auburn University
 // Part of azplugins, released under the BSD 3-Clause License.
 
-/*! \file BondEvaluatorDoubleWell.h
-    \brief Defines the bond evaluator class for a double well potential
-*/
-
 #ifndef AZPLUGINS_BOND_EVALUATOR_DOUBLE_WELL_H_
 #define AZPLUGINS_BOND_EVALUATOR_DOUBLE_WELL_H_
 
-#ifndef NVCC
+#ifndef __HIPCC__
 #include <string>
 #endif
 
-#include "hoomd/HOOMDMath.h"
+#include "BondEvaluator.h"
 
-#ifdef NVCC
+#ifdef __HIPCC__
 #define DEVICE __device__
 #else
 #define DEVICE
 #endif
 
+namespace hoomd
+    {
 namespace azplugins
     {
 namespace detail
     {
 
-//! Class for evaluating the double well bond potential
-/*!
-This bond potential follows the functional form
-\f{eqnarray*}
-
-V_{\rm{DW}}(r)  =   \frac{V_{max}-c/2}{b^4} \left[ \left( r - a/2 \right)^2 -b^2 \right]^2 +
-\frac{c}{2b}\left(r-a/2\right)+c/2
-
-\f}
-which has two minima at r = (a/2 +/- b), seperated by a maximum at a/2 of height V_max when c is set
-to zero.
-
-The parameter a tunes the location of the maximal value and the parameter b tunes the distance of
-the two maxima from each other.  This potential is useful to model bonds which can be either
-mechanically or thermally "activated" into a effectively longer state. The value of V_max can be
-used to tune the height of the energy barrier in between the two states.
-
-If c is non zero, the relative energy of the minima can be tuned, where c is the energy of the
-second minima, the first minima value is at zero. This  causes a small shift in the location of the
-minima and the maxima, because of the added linear term.
-
-The parameters are:
-    - \a V_max (params.x) potential difference between the the first minima and maxima
-    - \a a (params.y) shift for the location of the V_max, the maximum is at approx. a/2
-    - \a b (params.z) scaling for the distance of the two minima at approx. (a/2 +/- b)
-    - \a c (params.w) potential difference between the two minima, default is c=0
-
-*/
-class BondEvaluatorDoubleWell
+//! Parameters of double well bond potential
+struct BondParametersDoubleWell : public BondParameters
     {
-    public:
-    //! Define the parameter type used by this bond potential evaluator
-
-    typedef Scalar4 param_type;
+#ifndef __HIPCC__
+    BondParametersDoubleWell() : r_1(0), r_diff(1.0), U_1(0), U_tilt(0) { }
 
-    //! Constructs the pair potential evaluator
-    /*!
-     * \param _rsq Squared distance beteen the particles
-     * \param _params Per type bond parameters of this potential as given above
-     */
-    DEVICE BondEvaluatorDoubleWell(Scalar _rsq, const param_type& _params)
-        : rsq(_rsq), V_max(_params.x), a(_params.y), b(_params.z), c(_params.w)
+    BondParametersDoubleWell(pybind11::dict v)
         {
+        r_1 = v["r_1"].cast<Scalar>();
+        r_diff = r_1 - v["r_0"].cast<Scalar>();
+        U_1 = v["U_1"].cast<Scalar>();
+        U_tilt = v["U_tilt"].cast<Scalar>();
         }
 
-    //! This evaluator doesn't use diameter information
-    DEVICE static bool needsDiameter()
+    pybind11::dict asDict()
         {
-        return false;
+        pybind11::dict v;
+        v["r_0"] = r_1 - r_diff;
+        v["r_1"] = r_1;
+        v["U_1"] = U_1;
+        v["U_tilt"] = U_tilt;
+        return v;
         }
+#endif
+
+    Scalar r_1;    //!< location of the potential local maximum
+    Scalar r_diff; //!< difference between r_1 and r_0 (location of first minimum)
+    Scalar U_1;    //!< Potential Potential maximum energy barrier between minima
+    Scalar U_tilt; //!< tunes the energy offset (tilt) between minima
+    }
+#if HOOMD_LONGREAL_SIZE == 32
+    __attribute__((aligned(16)));
+#else
+    __attribute__((aligned(32)));
+#endif
 
-    //! Accept the optional diameter values
-    /*!
-     * \param da Diameter of particle a
-     * \param db Diameter of particle b
-     */
-    DEVICE void setDiameter(Scalar da, Scalar db) { }
+//! Class for evaluating the double well bond potential
+/*!
+ * This bond potential follows the functional form
+ * \f{eqnarray*}
+ *
+ * U(r)  =  U_1\left[\frac{\left((r-r_1)^2-(r_1-r_0)^2\right)^2}{\left(r_1-r_0\right)^4}\right]
+ *  +
+ * U_{\rm{tilt}}\left[1+\frac{r-r_1}{r_1-r_0}-\frac{\left((r-r_1)^2-(r_1-r_0)^2\right)^2}{\left(r_1-r_0\right)^4}\right]
+ *
+ * \f}
+ * which has two minima at r = r_0 and r = 2 * r_1 - r_0, seperated by a maximum
+ * at r_1 of height U_1 when U_tilt is set to zero.
+ *
+ * The parameter r_1 tunes the location of the maximal value and the parameter r_0 tunes the
+ * distance of the two minima from each other.  This potential is useful to model bonds which can be
+ * either mechanically or thermally "activated" into a effectively longer state. The value of U_1
+ * can be used to tune the height of the energy barrier in between the two states.
+ *
+ * If U_tilt is non zero, the relative energy of the minima can be tuned, where 2 * U_tilt
+ * is the energy of the second minima, the first minima value is at zero. This causes a
+ * small shift in the location of the minima and the maxima, because of the added linear term.
+ */
+class BondEvaluatorDoubleWell : public BondEvaluator
+    {
+    public:
+    typedef BondParametersDoubleWell param_type;
 
-    //! This evaluator doesn't use charge
-    DEVICE static bool needsCharge()
+    DEVICE BondEvaluatorDoubleWell(Scalar _rsq, const param_type& _params)
+        : BondEvaluator(_rsq), r_1(_params.r_1), r_diff(_params.r_diff), U_1(_params.U_1),
+          U_tilt(_params.U_tilt)
         {
-        return false;
         }
 
-    //! Accept the optional charge values
-    /*!
-     * \param qa Charge of particle a
-     * \param qb Charge of particle b
-     */
-    DEVICE void setCharge(Scalar qa, Scalar qb) { }
-
-    //! Evaluate the force and energy
-    /*!
-     * \param force_divr Output parameter to write the computed force divided by r.
-     * \param bond_eng Output parameter to write the computed bond energy
-     *
-     *  \return True if they are evaluated or false if the bond energy is not defined.
-     */
     DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& bond_eng)
         {
         bond_eng = 0;
         force_divr = 0;
-
-        // check for invalid parameters
-        if (b == Scalar(0.0))
+        // check for invalid parameters r0 = r1
+        if (r_diff == Scalar(0.0))
             return false;
 
-        Scalar c_half = Scalar(0.5) * c;
-        Scalar r = fast::sqrt(rsq);
-        Scalar r_min_half_a = r - Scalar(0.5) * a;
-        Scalar b_sq = b * b;
-        Scalar d = r_min_half_a * r_min_half_a - b_sq;
-
-        bond_eng = ((V_max - c_half) / (b_sq * b_sq)) * d * d + c_half / b * r_min_half_a + c_half;
-        force_divr = -(4 * (V_max - c_half) / (b_sq * b_sq) * d * r_min_half_a + c_half / b) / r;
+        const Scalar r = fast::sqrt(rsq);
+        const Scalar x = (r_1 - r) / r_diff;
+        const Scalar x2 = x * x;
+        const Scalar y = Scalar(1.0) - x2;
+        const Scalar y2 = y * y;
 
+        bond_eng = U_1 * y2 + U_tilt * (Scalar(1.0) - x - y2);
+        force_divr = (Scalar(4.0) * x * y * (U_tilt - U_1) - U_tilt) / (r_diff * r);
         return true;
         }
 
-#ifndef NVCC
-    //! Get the name of this potential
-    /*!
-     * \returns The potential name. Must be short and all lowercase, as this is the name energies
-     * will be logged as via analyze.log.
-     */
+#ifndef __HIPCC__
     static std::string getName()
         {
-        return std::string("doublewell");
+        return std::string("DoubleWell");
         }
 #endif
 
-    protected:
-    Scalar rsq;   //!< Stored rsq from the constructor
-    Scalar V_max; //!< V_max parameter
-    Scalar a;     //!< a parameter
-    Scalar b;     //!< b parameter
-    Scalar c;     //!< c parameter
+    private:
+    Scalar r_1;    //!< r_1 parameter
+    Scalar r_diff; //!< r_diff parameter
+    Scalar U_1;    //!< U_1 parameter
+    Scalar U_tilt; //!< U_tilt parameter
     };
 
     } // end namespace detail
     } // end namespace azplugins
+    } // end namespace hoomd
 
 #undef DEVICE
 

src/CMakeLists.txt

@@ -3,8 +3,10 @@ set(COMPONENT_NAME azplugins)
 
 # TODO: List all host C++ source code files in _${COMPONENT_NAME}_sources.
 set(_${COMPONENT_NAME}_sources
+    ConstantFlow.cc
     GroupVelocityCompute.cc
     module.cc
+    ParabolicFlow.cc
     )
 
 # TODO: List all GPU C++ source code files in _${COMPONENT_NAME}_cu_sources.
@@ -14,16 +16,22 @@ set(_${COMPONENT_NAME}_cu_sources
 # TODO: List all Python modules in python_files.
 set(python_files
     __init__.py
+    conftest.py
     bond.py
     flow.py
     pair.py
     )
 
 # TODO: Add names of all bond evaluators
-set(_bond_evaluators )
+set(_bond_evaluators
+    DoubleWell
+    )
 
 # 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 )

src/ConstantFlow.cc

@@ -0,0 +1,34 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2024, Auburn University
+// Part of azplugins, released under the BSD 3-Clause License.
+
+#include "ConstantFlow.h"
+
+namespace hoomd
+    {
+namespace azplugins
+    {
+namespace detail
+    {
+void export_ConstantFlow(pybind11::module& m)
+    {
+    namespace py = pybind11;
+    py::class_<ConstantFlow, std::shared_ptr<ConstantFlow>>(m, "ConstantFlow")
+        .def(py::init<Scalar3>())
+        .def_property(
+            "velocity",
+            [](const ConstantFlow& U)
+            {
+                const auto field = U.getVelocity();
+                return pybind11::make_tuple(field.x, field.y, field.z);
+            },
+            [](ConstantFlow& U, const pybind11::tuple& field)
+            {
+                U.setVelocity(make_scalar3(pybind11::cast<Scalar>(field[0]),
+                                           pybind11::cast<Scalar>(field[1]),
+                                           pybind11::cast<Scalar>(field[2])));
+            });
+    }
+    } // end namespace detail
+    } // end namespace azplugins
+    } // end namespace hoomd

src/ConstantFlow.h

@@ -5,28 +5,40 @@
 #ifndef AZPLUGINS_CONSTANT_FLOW_H_
 #define AZPLUGINS_CONSTANT_FLOW_H_
 
+#ifndef __HIPCC__
+#include <pybind11/pybind11.h>
+#endif
+
 #include "hoomd/HOOMDMath.h"
 
-#ifdef __HIPCC__
+#ifndef __HIPCC__
 #define HOSTDEVICE __host__ __device__
 #else
 #define HOSTDEVICE
-#include <pybind11/pybind11.h>
+#endif // __HIPCC__
+
+#ifndef PYBIND11_EXPORT
+#define PYBIND11_EXPORT __attribute__((visibility("default")))
 #endif
+
 namespace hoomd
     {
 namespace azplugins
     {
 
 //! Position-independent flow along a vector
-class ConstantFlow
+class PYBIND11_EXPORT ConstantFlow
     {
     public:
     //! Constructor
     /*!
      *\param U_ Flow field
      */
-    ConstantFlow(Scalar3 U_) : U(U_) { }
+    ConstantFlow(Scalar3 velocity)
+        {
+        setVelocity(velocity);
+        }
+
     //! Evaluate the flow field
     /*!
      * \param r position to evaluate flow
@@ -52,32 +64,10 @@ class ConstantFlow
     Scalar3 U; //!< Flow field
     };
 
-namespace detail
-    {
-void export_ConstantFlow(pybind11::module& m)
-    {
-    namespace py = pybind11;
-    py::class_<ConstantFlow, std::shared_ptr<ConstantFlow>>(m, "ConstantFlow")
-        .def(py::init<Scalar3>())
-        .def_property(
-            "mean_velocity",
-            [](const ConstantFlow& U)
-            {
-                const auto field = U.getVelocity();
-                return pybind11::make_tuple(field.x, field.y, field.z);
-            },
-            [](ConstantFlow& U, const pybind11::tuple& field)
-            {
-                U.setVelocity(make_scalar3(pybind11::cast<Scalar>(field[0]),
-                                           pybind11::cast<Scalar>(field[1]),
-                                           pybind11::cast<Scalar>(field[2])));
-            });
-    }
-    } // end namespace detail
-
     } // namespace azplugins
     } // namespace hoomd
 
 #undef HOSTDEVICE
+#undef PYBIND11_EXPORT
 
 #endif // AZPLUGINS_CONSTANT_FLOW_H_