A mean-field approach to simulating anisotropic particles

Particle simulation has become a crucial component of the scientific toolbox for investigating a wide range of physical phenomena, with various simulation techniques to choose from based on which ones are best suited to the phenomena of interest.^1,2 A core element for such simulations is the pair potential, which governs the pairwise interactions between any two particles. Pair potentials are ubiquitous in physics and chemistry, describing everything from the gravitational interaction of celestial bodies to the electrostatic forces between atoms and molecules. The straightforward application of pair potentials representing Coulombic and van der Waals forces to point particles representing atoms has been enormously successful in predicting a wide array of behaviors in highly complex systems, from biomolecules to organic compounds and metals.^3–10 However, despite algorithmic developments and improved hardware, the range of applicability of this approach for molecules or particles of complex geometry remains limited because of the complexity in computing all relevant interactions for sufficiently large numbers of particles.^6,11–13

Coarse-graining, the process of collectively handling many atomistic degrees of freedom, is a powerful strategy to increase the time and length scales accessible via simulation.^14,15 Examples of this include the bead-spring polymer model^16,17 and the Martini force field for proteins.¹⁴ However, coarse-grained models still utilize isotropic pair potentials similar in form to those used in all-atom simulations, depending upon bonded networks of isotropic particles to represent nontrivial geometries, such as proteins.^18,19 While such approaches assume that the internal flexibility and interactions within such networks are important, for many phenomena, these models can be simplified by treating each such network as a single rigid body,²⁰ for example, by using sets of overlapping spheres to represent a rod or cube.^21,22 Despite substantially reducing complexity and admitting more efficient simulation techniques,²³ these methods ultimately still require computation of a large number of interactions. Moreover, recent studies have shown that even coarser representations of particle shape are sufficient to predict thermodynamic and kinetic behavior,^24–26 a discovery that motivates the ongoing pursuit of new coarse-graining approaches and simulation methodologies for studying anisotropic particles.

Numerous different methodologies for the simulation of anisotropic particles depending on the phenomena of interest have emerged. Monte Carlo (MC) methods have been extremely successful in elucidating the assembly behavior of shaped particles.^24–26 Since MC methods do not require force fields, potential energies need not be represented by smooth or even continuous functions. This flexibility has been critical to previous studies of shaped particles, which have made extensive use of the hard particle approximation to bypass the need for even an analytical expression to represent interactions between particles.²⁷ However, because MC simulations rely on a statistical sampling rather than a time integration of equations of motion, in principle, they only provide information on equilibrium (thermodynamic), not kinetic, behavior unless we make further assumptions. Furthermore, the hard particle approximation captures all but one aspect of the behavior of attractive shapes since no overlap is permitted, and parallelizing hard particle Monte Carlo (HPMC) simulations to the degree possible in molecular dynamics (MD) remains an unmet challenge. As a result, an MD method for simulating anisotropic particles remains eminently desirable.

Perhaps the best-known method for MD simulations of anisotropic particles is the Gay–Berne (GB) potential, which rescales the distance and energy in the 12-6 Lennard-Jones (LJ) potential by the relative orientations of the pair of particles to simulate ellipsoids.^28–30 The potential is an analytic function that can be computed efficiently, but no known generalization of the Gaussian overlap formulation exists for arbitrary geometries, illustrating the central problem for anisotropic pair potential development. One method for circumventing this difficulty is using event-driven molecular dynamics (EDMD), which can be much faster than force-based MD in some cases^31,32 and has been used to great effect to simulate systems ranging in complexity from simple hard spheres³³ to fibril assembly from proteins.³⁴ However, although EDMD has been used to study the assembly of hard cubes,³⁵ there is no known generalization for arbitrary geometries.

An alternate approach to dynamical simulations of shaped particles is the discrete element method (DEM), in which geometries are decomposed into discrete pieces whose interactions can be independently computed and then summed. This method has been used to great effect in the granular materials community,³⁶ and it has been adapted to conform to the energy and momentum conserving schemes preferred in other particle simulation communities.³⁷ While the DEM can be used to study both dynamical and equilibrium behavior, it, too, suffers from certain limitations. The DEM is generally restricted to working with shapes that admit a decomposition into a discrete collection of elements, and while it readily admits parallelization of numerous calculations, the approach scales poorly to complex geometries where the number of discrete interaction sites becomes prohibitively expensive. As a result, we seek a method that can ameliorate these performance concerns while still accurately representing both smooth (e.g., ellipsoidal) and discontinuous (e.g., polytope) geometries.

In this paper, we introduce a framework for deriving anisotropic pair potentials. The resulting anisotropic pair potentials generalize existing isotropic pair potentials in a way that largely combines the strengths of the previously discussed methods. Since these pair potentials can be used in either MD or MC simulations, the method allows the study of both dynamical systems and systems out of equilibrium. Moreover, simulations using these pair potentials can take advantage of both the parallelization schemes used in MD and the cluster moves and advanced sampling techniques used to overcome kinetic barriers in MC. Moreover, by choosing an appropriate pair potential, the method can be used to study both attractive and purely repulsive particles. For discrete geometries such as polytopes, the non-smoothness associated with parallel faces and sharp edges and vertices introduces discontinuities into the potential energy landscape with respect to the orientational coordinates, but we show that these discontinuities can be ameliorated with a local averaging that achieves energy conservation comparable to the DEM while retaining advantageous performance characteristics.

We organize this paper as follows: In Sec. II, we describe the derivation of anisotropic pair potentials by applying a mean-field approximation to a form of Poisson’s equation that describes a wide range of interactions. The resulting pair potentials are applicable to a wide range of shapes, including ellipsoids and polytopes, which are our principal focus in this work. As an example, we derive an anisotropic generalization of the LJ potential – the anisotropic Lennard-Jones (ALJ) potential – along with a purely repulsive anisotropic Weeks–Chandler–Anderson (AWCA) variant. Section III discusses a practical implementation of this ALJ potential for use in MD simulations on both central processing units (CPUs) and graphical processing units (GPUs). In Sec. IV, we validate this implementation and measure its performance, focusing on the AWCA potential since there are directly comparable methodologies (primarily the DEM) that have been used for simulations of purely repulsive shapes in the past. We also demonstrate the utility of this method for studying particle self-assembly by demonstrating the assembly of some common colloidal crystal structures. Finally, we summarize and consider avenues for future work using this potential in Sec. V. This work makes extensive use of NumPy³⁸ and rowan,³⁹ and figures in this paper were generated using Matplotlib⁴⁰ unless otherwise noted. Data and workflows were managed with the signac framework.^41,42

Consider a system of N anisotropic particles occupying a volume V. The force on a reference particle with position $r→$ and orientation α from the potential energy field generated by all the other particles in some fixed configuration can be written as

where ∇U is the vector force field (the derivative of the underlying potential of interest), ρ defines the distribution of particles within the volume V, and s and β are the positions and orientations of the particles relative to the reference position and orientation, respectively.

For simplicity, we restrict ourselves to shapes that can be constructed by the anisotropic dilation of a sphere contained in that shape. Mathematically, this motivates defining a new coordinate system (r, θ, ϕ) such that x = rf_x(θ, ϕ), y = rf_y(θ, ϕ), and z = rf_z(θ, ϕ), where r is the radial distance from the sphere center and f_x, f_y, and f_z are some functions encoding the angular dependence of x, y, and z, respectively. This parameterization exists if and only if there is a sphere of radius r (not necessarily unique) for which a unique mapping (θ, ϕ) → (x, y, z) exists. Locally, the magnitude of the dilation of the sphere to the shape at any point (r, θ, ϕ) is given by $Ω(θ,ϕ)=fx2+fy2+fz21/2$⁠. This parameterization exists for all convex shapes, and more generally for any shape for which a sphere (not necessarily unique) may be chosen such that Ω(θ, ϕ) is unique relative to the center of that sphere. An example of a nonconvex shape that admits such parameterization is a star; conversely, an example for which this parameterization is impossible is a U-shaped particle (see Sec. III for how to simulate such a particle in practice).

In this work, we define the insphere of a shape as the largest sphere whose center coincides with the centroid of the shape that is also contained inside the shape. In the subsequent derivation, the coordinate system is always defined relative to this insphere. As a clarifying example, for a cube whose insphere radius is 1, Ω = 1, $2$⁠, and $3$ at the face centers, edge centers, and corners, respectively. Such an approach has been successfully employed both to develop an equation of state for hard polyhedra⁴³ and to predict ligand shell morphology about an anisotropic particle.⁴⁴ Transforming to these new coordinates, we find

where the Jacobian determinant of the coordinate transformation is simply r²Ω(β) and Ω is a function of the relative orientation of each particle. Taking the divergence of both sides of Eq. (2) results in a form of Poisson’s equation,

where U_in represents the radial component of the potential U. This component can be separated due to two important features of Eq. (3). First, because we assume that U is conservative, all non-radial components of the divergence vanish by symmetry. Second, because we restricted ourselves to shapes that can be modeled by anisotropic dilations of a sphere, U is separable into its radial and angular components, and the θ and ϕ terms in the radial derivative cancel. A crucial consequence of these two features is that all orientational dependence on the right-hand side of the equation is confined to the Ω term, allowing us to consider the derivative of only the radial part of the potential on the right side.

Therefore, we can consider a self-consistent field solution to Eq. (3) using an isotropic pair potential of interest as the input potential U_in, where U_in could be any of the isotropic pair potentials commonly used in simulations. There is no known analytical solution to Eq. (3) because no general expression exists for Ω across the class of anisotropic particles. However, in simulations, we only need to evaluate U for a specific configuration, for which we can employ iterative approaches to determine Ω for a given pair of particles even if no analytical solution is known for computing the distance between that particular pair of shapes. As a result, we can treat Ω as a constant for a given configuration and solve Eq. (4) for the effective pairwise interaction between anisotropic particles, given an input potential U_in.

With a constant Ω, Eq. (3) may be solved numerically given an expression for the density, so we now make three further approximations to derive such an expression. First, rather than solving for the potential in a particular fixed configuration, we instead perform an ensemble average to replace U with the potential of mean force U_m, allowing us to approximate the system density distribution as $ρ=ρoe−βUm$⁠. Second, although U_m is a configurational integral over all positions and orientations, we assume that the field strength at any point in space is dominated by the closest point (θ_c, ϕ_c) and perform a first-order Taylor expansion of U_m about this point. Third, we Taylor expand the exponential itself, yielding the following differential equation:

where U₀ is the first term in the Taylor expansion. For a system of interacting anisotropic particles, the Taylor expansion of U_m amounts to assuming that the field strength at any point in space is dominated by the closest point (θ_c, ϕ_c) on each particle to that point in space. The corresponding distance r_c, which we term the “contact” distance, is a function of the orientation α, but the formulation of Eq. (4) allows us to instead view this equation in terms of the parameters r and r_c, both of which are easily computed from simulations. We are now in a position to solve Eq. (4) for a given isotropic pair potential. For our purposes, we choose the traditional LJ potential. When U_in = U_LJ, the numerical solution to Eq. (4) is – not surprisingly – similar in shape to an LJ potential, shifted along the x axis by Ω [Fig. 1(a)]. This suggests employing the LJ functional form as a basis function for defining approximate analytical solutions for U_m. Using this basis function and a change of variables allow us to solve for U_m. Requiring that U₀ reflects the same anisotropic shape as U_m provides a condition that can be used to solve for U₀, resulting in the following approximate form for the ALJ potential U_ALJ = U₀ + U_m:

where we call the U₀ term the central interaction (because it depends only on r) and the U_m term the contact interaction (because it depends on r_c). σ_c is a free parameter that controls the relative strength of the contact interaction and determines how well the potential represents the shape of interest. In the limit σ_c → 0, the shape will be sharply preserved, but the resulting potential will be extremely hard and require very small time steps to simulate. Conversely, large values of σ_c will lead to an effectively softer (but much more rounded) representation of the shape. The requirement that U₀ must preserve the anisotropic core shape shown in Eq. (5) [Figs. 1(b) and 1(c)] results in defining $ϵs=ϵL1/L2$⁠, with

where σ_ij is the average of the insphere diameters of particles i and j and M_ij = (Ω_i(θ_c, ϕ_c) − 0.5σ_i) + (Ω_j(θ_c, ϕ_c) − 0.5σ_j). Figures 1(d) and 1(e) show the isocontours of U_ALJ if this rescaling is ignored, and we simply set ϵ_s = ϵ, in which case these potential energy surfaces bulge the faces of the shape in a way that does not reflect the underlying shape.

The core insights of our theory lie in Eqs. (2) and (5). The success of the coordinate transformation Eq. (2) in reducing Eq. (1) to the tractable form of Eq. (5) demonstrates that for a wide class of anisotropic particles (those for which the coordinate transformation is valid) and a wide class of potential forms (those for which the first-order Taylor expansion about the contact point is a good approximation), the total interaction reduces from a volume integral to a sum of two point interactions. While our primary usage of this construct in this work is to develop a new framework for computational simulations, it could also be used in future theoretical treatments of the interaction of anisotropic particles.

Although this potential includes both attractive and repulsive interactions, as discussed in Sec. I, the most commonly used methods for simulations of shapes at present have largely focused on hard shapes. For comparison with these existing methods, we define a short-ranged repulsive variant of our potential analogous to the Weeks–Chandler–Anderson (WCA) potential. This AWCA potential is simply a truncated and shifted ALJ potential, with the additional note that central and contact components of the potential are truncated at their respective minima (note that this also applies to the L₁ and L₂ terms defining ϵ_s). The results in all future sections use the AWCA potential rather than the ALJ potential unless noted otherwise.

We have implemented the ALJ pair potential in the simulation engine HOOMD-blue.^45,46 As described in Eq. (5), the ALJ potential is the sum of a central potential U₀ and a contact potential U_m evaluated between the closest points on a pair of particles (A, B). To find the closest contact points {a ∈ A, b ∈ B}, we employ a modified version of the Gilbert–Johnson–Keerthi (GJK) distance algorithm. We define torques $τi→=r→i×Fm→$⁠, where F_m = −∇U_m and the moment arms are given by $r→a=a−qa$ and $r→b=b−qb$ (where q_i is the centroid of particle i). As noted in Sec. II, this pair potential is only defined for shapes that are amenable to the coordinate transformations described there. However, in practice, arbitrary nonconvex particles may be simulated by using rigid bodies composed of convex ALJ particles.

Finding center–center interparticle distances is handled by the existing machinery in HOOMD-blue used by other pair potentials. The most complex part of implementing this anisotropic pair potential is, therefore, developing a sufficiently performant and robust implementation of the GJK algorithm. In order to discuss our implementation and some of its particular features, we provide a brief introduction to the GJK distance algorithm here that covers the crucial details. We refer the interested reader to Refs. 47 and 48 for more extensive descriptions of the GJK distance algorithm (the following sections closely follow the notation of Ref. 47).

The GJK algorithm^49,50 for computing the minimum distance between two arbitrary convex sets is based on two crucial geometric concepts: Minkowski sums and support functions. The central idea is that given two convex sets A and B, the vector joining the closest points on the two shapes is equal to the vector defining the closest point to the origin on the Minkowski sum of A and −B (defined as A − B = {a − b: a ∈ A, b ∈ B}); consequently, when two shapes overlap, the origin must be contained in A − B. Support functions are the critical ingredients that provide an implicit representation of A − B via their supporting hyperplanes, bypassing the otherwise prohibitively expensive explicit computation of the Minkowski sum. The support function $hK(x→)=sup{x→⋅y:y∈K}$ of a convex set K provides the distances between its supporting hyperplanes and the origin, and the corresponding support mapping $sK(x→)$ defined by the relation $hK(x→)=sK(x→)⋅x$ gives the corresponding point on the hyperplane.

In essence, the GJK algorithm is a descent algorithm that constructs a sequence of simplices (W₀, W₁, W₂, …) contained in A − B and a corresponding sequence of vectors $(v→0,v→1,v→2,…)$ that is guaranteed to converge to the point of minimum norm $v→*$⁠. At each iteration, the point $w→k=sA−B(vk→)$ is added to W_k, and $v→k+1$ is set to the closest point to the origin in the augmented W_k. The smallest subset of this set that contains $v→k+1$ then becomes W_k+1. Finding this subset efficiently is a critical part of the GJK algorithm and is discussed in more detail below.

We implemented a version of the GJK algorithm with a few noteworthy features that stem from our specific requirements. In addition to the vector $v→*$⁠, we also require the points $a→*∈A$ and $b→*∈B$ that correspond to the endpoints of $v→*$ to apply the torques. To find these points, we note that

such that ∑_iλ_i = 1 and ∀i (λ_i > = 0). These equations imply that $v→k+1$ is a convex combination of ${sA−B(v→k)}$⁠. Defining the index set $Ik={i:v→i∈Wk}$⁠, we can find the necessary points using

Therefore, to find $a→*$ and $b→*$⁠, we track the support evaluations for the two shapes s_A and s_B corresponding to the current $v→k$⁠. When the algorithm converges, we use these evaluations and the computed λ_i values to recover the contact points $a→*$ and $b→*$⁠, and the difference gives us $v→*$⁠.

Accurately estimating $a→*$ and $b→*$ in addition to $v→*$ using the classical GJK algorithm is infeasible due to numerical inaccuracies in the original method for determining W_k+1 from W_k, which is known as the Johnson (JH) subalgorithm. This solves a linear system of equations to find the minimal subset W_k+1 from W_k, but near-degenerate configurations are known to cause severe robustness issues for the method, and we found that even the most robust termination conditions⁴⁷ for the JH subalgorithm led to solutions that, in the worst cases, deviated from $v→*$ by up to 5%. More critically, in many cases, even much smaller errors in $v→*$ led to dramatically incorrect approximations of $a→*$ and $b→*$⁠, frequently leading to torques in the wrong direction and causing unacceptably severe violations of energy and momentum conservation. To solve this problem, the original GJK algorithm paper proposes an expensive backup procedure that essentially amounts to exhaustive testing of all possible points in the set,⁴⁹ but this procedure is prohibitively expensive for our purposes, particularly on GPUs, where it introduces severe warp divergence.

To ameliorate this issue, we implemented the recently developed signed volumes (SV) subalgorithm.⁵¹ This subalgorithm is substantially more robust and avoids the incorrect torques resulting from the JH one. It also significantly improves performance, both because it accelerates convergence and because it requires far less caching of intermediate calculations, reducing memory requirements in otherwise compute-bound GPU kernels. We observed speedups of nearly 30% on the GPU even for simple polyhedral geometries where the total number of iterations is tightly bounded by the sum of the number of vertices of both shapes. We expect even larger improvements for curved geometries where the SV subalgorithm improves convergence; however, due to the stability issues encountered in simulations using the JH one, we did not conduct this performance comparison.

To improve energy and momentum conservation, we enforce a strict ordering on evaluations of the GJK algorithm. On the GPU, HOOMD-blue computes all interactions twice to avoid communicating between threads, but due to limitations of finite precision arithmetic, differences between gjk(i, j) and gjk(j, i) are enough to lead to violations of momentum conservation even with convergence tolerances on the order of machine precision. To avoid this issue, we enforce that for any pair, U(i, j) = U(j, i) by evaluating gjk(i, j) if i < j and otherwise by evaluating gjk(j, i). This ensures that Newton’s third law is obeyed [F(i, j) = −F(j, i)] and momentum is conserved. We tested alternative estimators for the true distance, including min(gjk(i, j), gjk(j, i)) and mean(gjk(i, j), gjk(j, i)), but we found no measurable difference in the energy conservation of the system over the tested timescales. We therefore chose the current minimum index heuristic for computational efficiency since it only requires a single evaluation of the distance algorithm.

The contact point approximation breaks down for polytopes because there can be multiple pairs of points between two polytopes whose pairwise distances are arbitrarily similar. Although forces computed remain reasonable because the contact distance is a smooth function, changes in the contact point can be arbitrarily large relative to the simulation timestep, introducing a discontinuity in the energy with respect to the orientational coordinates of the system. For a simple case illustrating this discontinuity, consider the interaction between a pair of edge-aligned squares in 2D. In this configuration, an arbitrarily small rotational move of either particle about its respective centroid would result in a shift of the contact point from one vertex to the other. Since this shift is not a function of the integration timestep, the change in torque can become arbitrarily large relative to the move size, leading to significant violations of energy conservation.

To mitigate this issue, instead of only computing interactions about the contact point, we compute a sum of all interacting features a la the DEM.³⁷ In 2D, we compute all vertex-edge interactions between polygon edges, whereas in 3D, we compute both vertex–face and edge–edge interactions between polyhedra. To avoid unnecessary computation of negligible interactions, we construct the vector joining the centroid of one shape to the contact point on the other shape (found using the GJK algorithm) and identify the face this vector passes through. In all relevant cases, this face accounts for all interactions that are necessary to consider. This method, therefore, leads to energy conservation comparable to the DEM while ensuring that performance scales with the number of vertices in the largest face of the polytope rather than scaling with the total number of vertices (which is the principal reason for poor performance scaling in the DEM).

Note that for pathological geometries where multiple faces on a shape are arbitrarily close to parallel, averaging over a single face may not be enough. In this case, a local search of the vertex graph to find and test distances for all facets neighboring the closest facet pair would be necessary to ensure perfect equivalence with the DEM. Since none of the current shapes of interest meet this criterion, we omit this additional smoothing procedure here. Currently, this averaging is enabled by default in the implementation in HOOMD-blue but may be turned off by users in Python scripts in situations with less restrictive energy conservation requirements. All presented results are gathered with the full face averaging unless stated otherwise.

We first validated our method for smooth shapes, where the contact-point formalism should hold exactly. In Fig. 2, we show long-time energy conservation for a system of AWCA ellipsoids; over the timescales observed, we see no measurable drift. It is worth noting that simulations using the GB potential to represent equivalent ellipsoids result in relative fluctuations that are 1–2 orders of magnitude smaller. This difference can largely be attributed to that fact that, since the surface component of the AWCA potential is extremely hard (it is essentially an LJ potential with a very small σ that has been shifted to the surface), small movements result in much larger changes in energy. There may also be small effects due to the use of the GJK algorithm to compute the contact point since any numerical error in this calculation could cause larger shifts in the computed distance from one timestep to the next than the true change. However, our implementation is accurate up to machine precision,⁵¹ and we observed no measurable drift over the tested timescale, suggesting that any such fluctuations would be very minor, if present. Note that the true shape of the ellipsoid is determined by the sum of the underlying principal axes and the rounding radii; for instance, an ellipsoid specified with a = 3, b = 2, c = 1, and σ_c = 0.1 would effectively be an ellipsoid with principal axes a = 3.1, b = 2.1, and c = 1.1. Any desired ellipsoid can still be achieved by setting the principal axes of the internal ellipsoid appropriately.

We next validated the method for AWCA polytopes. All simulations of polyhedra in this paper were conducted using contact spheres of size σ_c = 0.15σ, where σ is defined as the insphere diameter of the shape and is used to evaluate the central potential. For instance, for a unit volume cube, the diameter of the insphere is unity, so σ_c = 0.15. Although this choice makes the potential very hard and necessitates a very small dt, we chose it to retain very accurate representations of our shape. It has previously been shown that a judicious choice of the rounding radius can still result in desired assemblies while also allowing timesteps one to two orders of magnitude larger,³⁷ so such choices can be made to accelerate simulations on a case-by-case basis.

Without the averaging described in Sec. III B, system energies drift significantly even over relatively short times in both 2D and 3D simulations of polytopes. Simulations in HOOMD-blue can be set to completely ignore anisotropic degrees of freedom by setting the parameter aniso = False in the integrator. When doing so, we observe energy conservation on par with LJ sphere simulations, indicating that the discontinuity is indeed entirely in the rotational motion of the system.

In Fig. 3, we show the result of simulations that employ the face averaging described in Sec. III B. In general, the energy conservation of the AWCA is on par with that of the DEM. While proper face averaging is sufficient to avoid energy drift in two dimensions, we note that that in three dimensions, both the DEM and the face-averaged AWCA show significant long-time drift [see Fig. 3(b)]. This drift is the result of the edge–edge interactions, which cause small discontinuities in the torques for near-parallel edges because the closest point jumps from one vertex to the other (essentially the same issue that arises in simulations without face averaging). However, the edge–edge interactions are necessary to account for all possible collision configurations: for instance, a pair of shapes whose closest points lie on perpendicular edges would exhibit incorrect collision dynamics if only vertex–face interactions were considered.

Since any attempt to smooth this interaction imposes an effective distortion or rounding on the shape, we prioritize faithfully reproducing the shape since (as discussed in Ref. 37) the measured energy drift is acceptable for coarse-grained simulations for which this method (like the DEM) is most likely to be used. In fact, the AWCA energies in these plots are nearly identical to the DEM energies, which conforms to our expectations: in repulsive systems, the cutoff of the AWCA potential is short enough that the contact term dominates most pair interactions. As a result, the behavior of the two methods is essentially identical (although the AWCA potential is much faster in many cases as we will show). However, we note that this is not the case when the attractive, long-ranged ALJ potential is considered; in that case, the center–center interaction becomes the dominant contribution because at large separation distances, the sharply decaying contact term quickly becomes negligible.

We next discuss the performance of our method. Since our method is implemented in HOOMD-blue, it inherits most of the advantageous properties such as parallel execution of the message passing interface (MPI) on multiple nodes and strong scaling on GPUs.⁴⁶ The main distinction between our potential and the typical pair potentials in HOOMD-blue is that the ALJ potential is substantially more compute- and memory-intensive to evaluate. As a result, whereas in typical simulations of LJ particles (for instance) significant portions of the simulation time are taken up by tasks like neighbor list builds, in ALJ simulations, the energy evaluation becomes the bottleneck, so its intrinsic properties become very important to the feasibility of simulations.

For the simulation of polytopes, we again compare with the DEM method since it is the most similar method from this perspective. There are two main differences that we expect to govern the performance differences between the two methods. The first difference is in the operation count of the two methods. The DEM computes interactions between all pairs of features, so in general, there are $O(N2)$ distinct interactions to compute between two N-polytopes. Meanwhile, the AWCA potential uses the GJK algorithm to find relevant interactions to compute, so there are fewer total interactions computed; however, the GJK algorithm is expensive and adds a constant additional cost for every pair potential evaluation. Therefore, in general, we expect the DEM to perform better for small systems of simple shapes where computing all interactions in parallel is cheaper than evaluating with the GJK algorithm, whereas the AWCA potential should scale better to very complex shapes because the number of distinct interactions increases much more slowly than in the DEM (as we see in Fig. 4).

The second major difference is that the two use different parallelization schemes. The AWCA is parallelized like any other pair potential, i.e., it is evaluated in parallel for every pair of particles. In DEM simulations, all possible pairs of interaction sites between each pair of particles are known a priori, so the DEM can take advantage of this information to exploit an additional level of parallelism and compute all possible interaction pairs in parallel (of which there are multiple per particle pair). Since the AWCA potential uses the GJK algorithm to determine the face to average over each time the potential is evaluated for a given pair, this level of parallelization is not possible. Therefore, we expect the AWCA to be especially fast compared to the DEM on the CPU, whereas this advantage will be reduced on the GPU.

Figure 4 shows the investigations of the practical effect of these tradeoffs for various system sizes on both CPU and GPU architectures. In two-dimensional systems, the performance of the two methods is of the same order of magnitude. The DEM is approximately twice as fast for triangles and remains faster than the AWCA potential for all shapes with six or fewer vertices (triangles to hexagons). For shapes with more than six vertices, the AWCA potential becomes faster on the CPU with a relative speedup of ∼2× for most realistic shapes. On the GPU, however, the DEM retains a performance advantage up to higher vertex counts because of its more fine-grained parallelism. As a result, for smaller systems of 1000 particles, the DEM is essentially always faster by about a factor of 2. As system sizes increase and the GPU becomes saturated, however, the AWCA potential becomes faster since it has to do strictly less computation.

The performance difference is much more stark for three-dimensional systems. On the GPU, performance for the two methods is nearly equal for tetrahedra, and both cubes and octahedra experience a speedup of 2–3× when using the AWCA potential. However, more complex shapes run at least 5 times faster and usually more than 10 times faster, regardless of system size. On the CPU, the AWCA potential is essentially always faster, and for complex systems, the performance difference exceeds 1.5 orders of magnitude. The reason for this difference between two and three dimensions is that in three dimensions, there are substantially more interactions to account for. In two dimensions, there are exactly 2N² vertex-edge pairs for which calculations need to be performed, but in three dimensions, the total number of interactions is 2NF + 2EE, where N is the number of vertices, F is the number of faces, and E is the number of edges. Not only does this number grow much more quickly for a given N-polytope in three dimensions but each distance calculation is also more expensive than the corresponding calculations in two dimensions, accounting for the different performance characteristics between 2D and 3D. In both 2D and 3D, once system sizes are sufficiently large to saturate the GPU, we see that the relative performance of the AWCA potential increases as the effect of the DEM’s greater parallelism can no longer keep up with the total number of calculations to be performed. System size is far less impactful on the CPU since the total number of operations scales identically for both methods.

We omit quantitative comparison to the GB potential for ellipsoids in this paper except to note that at present, the AWCA potential is 1–2 orders of magnitude slower. This performance difference can be almost entirely attributed to the cost of using the GJK algorithm to find the contact distance. The GJK algorithm could easily be replaced with an exact solution for uniaxial ellipsoids or more efficient specialized approximate solutions for biaxial ellipsoids.²⁹ Moreover, since the neighbor list constructed by HOOMD-blue is based on the major axis of the ellipsoid, it includes many ellipsoids that in fact should not interact, given their orientations. One possibility for accelerating this would be filtering the neighbor list using, for instance, oriented bounding box (OBB) overlap checks,^55,56 but this method would need to be carefully implemented and profiled to understand the tradeoffs since such OBB overlap checks are also expensive calculations. Since the GJK algorithm maintains a lower bound for the distance as it iterates, the algorithm could also be modified to terminate early if the lower bound is ever larger than the distance cutoff; however, the resulting warp divergence renders this method relatively ineffective on the GPU. Since the immediate applications of our ALJ potential are for studying polyhedral particles, we leave such optimizations for future works specifically focused on studying ellipsoidal particles.

As discussed in Sec. I, some of the most successful methods for studying anisotropic particles are HPMC methods that enable the study of purely entropic systems driven by excluded volume effects alone. Since such methods have been used in prior works to elucidate the assembly propensities of various shapes,²⁴ we validate our method by showing that MD simulations using the AWCA potential can correctly reproduce these behaviors. We ran simulations of 1000 particles of various shapes in the NVT ensemble using a Nose–Hoover thermostat, as described by Martyna et al.⁵⁷ Simulations were conducted at a temperature of unity in units of $kBTϵ$⁠. Particles are modeled using the vertices of a unit volume polyhedron, the σ of the AWCA interaction is set to the insphere radius of the polyhedron, and the contact interaction length is σ_c = 0.15σ. Particle volumes, masses, and inertia tensors were calculated from the Minkowski sum of the core polyhedral particle and the contact sphere. We initialized systems in a dilute gas and compressed them to an effective packing fraction of 0.55. Particles were assumed to be of unit density.

In Fig. 5, we show the results for shapes with well-known assembly propensities, namely, the shapes corresponding to the Voronoi cells of SC, BCC, and FCC crystal structures. The cubes shown in Fig. 5(a) form an SC structure, the rhombic dodecahedra shown in Fig. 5(b) form an FCC structure, and the truncated octahedra shown in Fig. 5(c) form a BCC structure, as expected. The crystal structures are identified using polyhedral template matching (PTM), which is also used to color the particles. RDFs and diffraction patterns (computed using freud⁵⁸) are included as further evidence of structure formation.

Studies of anisotropic particles have appeared frequently in the literature in recent years. The role of shape in various phenomena has been examined by a number of methods, but thus far, researchers have struggled to find a general framework within which such shapes could be studied using classical MD simulations. In this paper, we presented such a framework and showed how it can be used to generalize a wide range of MD pair potentials to account for anisotropic shapes. We showed how this method can be used to study the behavior of nearly hard particle shapes, and we discussed its performance relative to other methods currently in use.

Our method enables the study of mixtures of arbitrary shapes, substantially generalizing existing dynamical methods and providing information that cannot be obtained from existing MC techniques. Concave shapes can be studied by constructing rigid bodies based on the convex decomposition of such shapes, and such rigid bodies would require far fewer particles than would be necessary using rigid bodies made of spheres. We are currently using this new method to study how particle anisotropy modulates dynamical behavior and will publish this research in a future article. Among other uses, we also foresee that this potential will be used to study the assembly behavior of attractive anisotropic particles and the properties of heterogeneous mixtures of anisotropic particles in the near future.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This research was supported in part by the National Science Foundation, Division of Materials Research Award No. DMR 1808342 (implementation and validation), and by the Department of the Navy, Office of Naval Research, under ONR Award No. N00014-18-1-2497 (theory). V.R. acknowledges the 2019–2020 J. Robert Beyster Computational Innovation Graduate Fellowship from the College of Engineering, University of Michigan. This research was supported in part through computational resources and services supported by Advanced Research Computing at the University of Michigan, Ann Arbor. Hardware provided by NVIDIA Corp. is gratefully acknowledged.