Adding Chebyshev anisotropic pair potential by mzf0069 · Pull Request #121 · mphowardlab/azplugins

mzf0069 · 2026-01-26T22:36:43Z

No description provided.

mphoward

This is a good start! Here are some comments, let me know what questions you have.

mphoward

This is going in a good direction. Let's finalize the interface, especially the constructor, then start to implement that.

mphoward · 2026-02-10T02:55:31Z

+                                            const std::vector<Scalar>& r0_data,
+                                            const std::vector<unsigned int>& terms,
+                                            const std::vector<Scalar>& coeffs);
+


r0_data should be a raw pointer to Scalar so that we can copy the data in-place from a NumPy array. You will need a second raw pointer or a std::vector that stores the dimensionality of r0_data along each of the 5 dimensions.

terms and coeffs also need to be raw pointers, and you'll need one more variable that is the number of terms.

Last, domain should probably also be a raw pointer to a regular Scalar so you can accept a NumPy array that is (6,2) in shape. Otherwise, it can be a std::vector<Scalar2>, and you can copy into it.

mphoward · 2026-02-10T02:57:25Z

+        }
+
+    /// r0 data (theta, phi, alpha, beta, gamma) (N x 6)
+    const GPUArray<Scalar>& getR0Data() const


Do you want to be expose the R0Data data directly after you construct the class? Same questions go for some of the members below that return GPUArray. You have to decide whether you are OK with the user modifying their contents or not.

mphoward · 2026-02-10T02:57:57Z

+    /// Allocate storage for term list and coefficients.
+    void resizeTerms(unsigned int Nterms);


Is this something we want the user to be able to do? It seems like you wouldn't want them to do that themselves, most likely. If it's a method the class needs to trigger, though, it can be protected.

mphoward · 2026-02-10T02:58:19Z

+    protected:
+    void computeForces(uint64_t timestep) override;
+
+    private:


Most of the members in this section will need to be protected not private if the GPU version of the class is going to access them.

mphoward

Looks good to me, let's start on the implementation!

mphoward · 2026-02-17T15:12:59Z

+                                      unsigned int Nterms,
+                                      const unsigned int* terms,
+                                      const Scalar* coeffs);


I like to follow the convention that the number of items comes after the pointers

Suggested change

unsigned int Nterms,

const unsigned int* terms,

const Scalar* coeffs);

const unsigned int* terms,

const Scalar* coeffs,

unsigned int Nterms);

mphoward · 2026-02-17T15:14:32Z

+                                      const Scalar* r0_data,
+                                      const unsigned int* r0_shape,


Logically, should these come before the things related to the potential approximation or after? It doesn't really matter, beyond which way it makes more sense for you to read them.

mphoward · 2026-02-17T15:15:43Z

+    // neighbor list object
+    std::shared_ptr<hoomd::md::NeighborList> m_nlist;


Member variables should be documented like this in doxygen style

Suggested change

// neighbor list object

std::shared_ptr<hoomd::md::NeighborList> m_nlist;

std::shared_ptr<hoomd::md::NeighborList> m_nlist; //!< Neighbor list

mphoward · 2026-02-17T15:17:54Z

+    // intenal r0 linear interpolator
+    std::unique_ptr<LinearInterpolator5D> m_r0_interp;
+
+    // r0_data
+    GPUArray<Scalar> m_r0_data;
+
+    std::array<unsigned int, 5> m_r0_shape;


We should discuss how the linear interpolator will work. It may not be necessary to store a separate variable for it, and instead just construct it each time it is needed because it will be a lightweight wrapper around the data and the shape.

Noting: if the shape is stored in an array, it will need to be copied into the interpolator to pass it to the GPU.

mphoward · 2026-02-17T15:18:17Z

+    GPUArray<Scalar> m_coeffs;
+
+    // number of terms
+    unsigned int m_Nterms = 0;


This doesn't need a default value assigned to it

Suggested change

unsigned int m_Nterms = 0;

unsigned int m_Nterms;

mphoward

This is moving in a good direction! Please see my comments and let me know what questions you have.

mphoward · 2026-03-06T16:57:50Z

+    domain = numpy.zeros((5, 2), dtype=numpy.float64)
+    terms = numpy.zeros((2, 6), dtype=numpy.uint32)


I would suggest giving the domain and terms some real values, not just zeros, as it may cause problems later once we start implementing things.

mphoward · 2026-03-06T16:59:13Z

+    assert pot._cpp_obj.n_terms == 2
+    assert numpy.isclose(pot._cpp_obj.r_cut, r_cut)


_cpp_obj is an implementation detail, so we should not test it in detail. The check above that it exists and is not None is OK though.

Suggested change

assert pot._cpp_obj.n_terms == 2

assert numpy.isclose(pot._cpp_obj.r_cut, r_cut)

mphoward · 2026-03-06T17:14:58Z

+    @property
+    def r_cut(self):
+        """Cut-off distance in approximation domain"""
+        return self._r_cut


You shouldn't need a property for the ones coming from C++, but you need to use the ParameterDict correctly for those.

Suggested change

@property

def r_cut(self):

"""Cut-off distance in approximation domain"""

return self._r_cut

mphoward · 2026-03-06T17:15:15Z

+    @property
+    def n_terms(self):
+        """Number of terms."""
+        return int(self._terms.shape[0])


Same here, this should come from C++ automatically.

Suggested change

@property

def n_terms(self):

"""Number of terms."""

return int(self._terms.shape[0])

mphoward · 2026-03-06T17:16:06Z

+    @property
+    def r0_shape(self):
+        """r0 table shape."""
+        return tuple(int(x) for x in self._r0_data.shape)


We will want to discuss how to expose the array quantities later on. However, I note for now that we don't want to have a special shape property for r0. It should be a NumPy inside this class, so it should work like r0.shape.

mphoward

This is moving in a very good direction!

mphoward · 2026-03-24T14:49:03Z

+        }
+
+        {
+        GPUArray<Scalar> r0_arr(n_r0, m_exec_conf);


Is GPUArray OK with a size_t number (no warnings)?

mphoward · 2026-03-24T14:56:08Z

+    // need to start from a zero force and torque
+    m_force.zeroFill();
+    m_torque.zeroFill();
+
+    ArrayHandle<Scalar4> h_force(m_force, access_location::host, access_mode::readwrite);
+    ArrayHandle<Scalar4> h_torque(m_torque, access_location::host, access_mode::readwrite);


Is this the same pattern used in PotentialPair? A more efficient memory access could be to grab these in overwrite, then memset to zero. But, Josh changed alot of those calls to use this zeroFill method, so we should follow the same pattern as him. Just double check.

zeroFill is used in PotentialPair. I kept it the same.

mphoward

This is really nicely done! Here are some suggestions, but they don't all necessarily need to be done now. Take a look and let me know what questions you have, then I think we should try the GPU implementation if you are happy with the CPU one.

mphoward · 2026-04-27T16:54:43Z

+    Scalar m_nlist_r_cut; //!< Neighbor-list cutoff = ceil(max(r0) + r_cut)
+
+    std::shared_ptr<GPUArray<Scalar>> m_r_cut_nlist; //!< r_cut matrix shared with nlist
+    bool m_attached = true;                          //!< Whether attached to the simulation


Why is this being declared here? I am not used to seeing this in HOOMD code for force computes, but maybe it was added later to code that uses neighbor lists.

I think PotentialPair is using the same pattern.

mphoward · 2026-04-27T17:03:00Z

+    Scalar max_r0 = *std::max_element(r0_data, r0_data + n_r0);
+    m_nlist_r_cut = std::ceil(max_r0 + m_r_cut);
+
+    m_r_cut_nlist = std::make_shared<GPUArray<Scalar>>(1, m_exec_conf);


Doesn't this need to be sized large enough to be all the types in the system? Probably you need to do that here and fill them all in.

Longer term, we should make a plan for how we deal with which particles this potential applies to.

mphoward · 2026-04-27T17:13:31Z

+    // Determine the maximum Chebyshev degree needed for each of the 6 coordinates
+    unsigned int max_deg[6] = {0, 0, 0, 0, 0, 0};
+    for (unsigned int t = 0; t < m_Nterms; ++t)
+        {
+        for (unsigned int c = 0; c < 6; ++c)
+            {
+            const unsigned int deg = h_terms.data[t * 6 + c];
+            if (deg > max_deg[c])
+                max_deg[c] = deg;
+            }
+        }
+
+    // Chain-rule scale factors: d(x_scaled)/d(x) = 2 / (hi - lo)
+    Scalar cheb_scale[6];
+    cheb_scale[0] = Scalar(2);
+    for (unsigned int d = 0; d < 5; ++d)
+        {
+        cheb_scale[d + 1] = Scalar(2) / (h_domain.data[d].y - h_domain.data[d].x);
+        }
+
+    // Flat 1D Chebyshev storage
+    unsigned int max_deg_global = 0;
+    for (unsigned int c = 0; c < 6; ++c)
+        {
+        if (max_deg[c] > max_deg_global)
+            max_deg_global = max_deg[c];
+        }
+
+    const Index2D cheb_idx(max_deg_global + 1, 6);
+    std::vector<Scalar> cheb_T_flat(cheb_idx.getNumElements());
+    std::vector<Scalar> cheb_dT_flat(cheb_idx.getNumElements());


This should be done during setup since it does not change. You can add a separate method if it is needed to keep the code cleaner.

Also, noting that the Index2D solution is fine for now, but we will probably want to use a different structure later to save memory on the GPU. These storage requirements will be much larger than necessary if you have an anisotropic grid.

mphoward · 2026-04-27T17:13:54Z

+    //! beta threshold for the Jacobian (avoids 1/sin(beta) singulrity).
+    const Scalar beta_tol = Scalar(1e-5);


This is coded into the lower bounds of the domain, no? Can probably remove and use that.

Or, maybe not—should the lower bounds of the domain just be the true 0s (for sampling purposes), and then the nonzero values are used here for evaluating?

Since we are representing orientation with ZXZ Euler angles and reducing configurations in the same coordinates, the force and torque expressions in these variables become singular at beta = 0 and beta = pi, where sin(beta) = 0. Because the interpolator is also constructed on a domain that excludes those singular endpoints, I think it is safest to keep the current nonzero lower bound in the symmetry class and use that same bound consistently during evaluation.

Add header for Chebyshev anisotropic pair potential

a631d5a

mzf0069 requested a review from mphoward January 26, 2026 22:37

mphoward reviewed Jan 30, 2026

View reviewed changes

Comment thread src/ChebyshevTensorAnisotropicPairPotential.h Outdated

Comment thread src/ChebyshevTensorAnisotropicPairPotential.h Outdated

Comment thread src/ChebyshevTensorAnisotropicPairPotential.h Outdated

Update header interface

c0323f1

mzf0069 requested a review from mphoward February 5, 2026 23:04

mphoward reviewed Feb 10, 2026

View reviewed changes

Update header interface

189f25f

mzf0069 requested a review from mphoward February 16, 2026 15:20

mphoward reviewed Feb 17, 2026

View reviewed changes

Update Chebyshev anisotropic pair potential and add linear interpolant

67e1bdc

mzf0069 requested a review from mphoward February 27, 2026 16:05

mphoward reviewed Mar 6, 2026

View reviewed changes

Update linear interpolation, and initiate compute force

45a1740

mzf0069 requested a review from mphoward March 19, 2026 16:29

mphoward reviewed Mar 24, 2026

View reviewed changes

mzf0069 added 2 commits April 10, 2026 13:54

Implement force and torque calculation

3c07536

Update linear interpolant

c1fc34a

mzf0069 requested a review from mphoward April 10, 2026 19:02

mzf0069 added 2 commits April 16, 2026 19:13

Update force and torque calculations and add symmetry evaluator.

e95b0dd

Update documentation

0d92038

mzf0069 requested review from mphoward and removed request for mphoward April 17, 2026 00:26

mzf0069 added 3 commits April 19, 2026 10:25

Add the symmetry class of cube

926b57f

Refactor Chebyshev anisotropic pair potential for symmetry support

c00f843

Remove domain as user-defined variable

a466061

mphoward reviewed Apr 27, 2026

View reviewed changes

Implement suggestions and start GPU implementation

0f7ec02

		/// Allocate storage for term list and coefficients.
		void resizeTerms(unsigned int Nterms);

		// neighbor list object
		std::shared_ptr<hoomd::md::NeighborList> m_nlist;

		domain = numpy.zeros((5, 2), dtype=numpy.float64)
		terms = numpy.zeros((2, 6), dtype=numpy.uint32)

		assert pot._cpp_obj.n_terms == 2
		assert numpy.isclose(pot._cpp_obj.r_cut, r_cut)

		//! beta threshold for the Jacobian (avoids 1/sin(beta) singulrity).
		const Scalar beta_tol = Scalar(1e-5);

Conversation

mzf0069 commented Jan 26, 2026

Uh oh!

mphoward left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

mphoward left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mphoward left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mphoward left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

mphoward left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

mphoward left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!