Evaluation of weighted Voronoi decompositions of physicochemical ensembles

Abstract

Voronoi diagrams are widely used to model disperse systems such as foams, powders, polycrystals and atoms in the classical limit. Voronoi tessellations partition the continuous phase into compartments, or cells, that encompass all space closer to the assigning dispersed object than any other in the system. To account for heterogeneity in object size, weights are applied to avoid unphysical partitioning across non-contacting objects. Power and additive weighting are the most common weighting schemes, wherein power is more computationally tractable but additive weighting correlates more directly with size. In general, the two schemes produce distinct spatial decompositions for any non-monodisperse system. To calibrate the divergent volumetric metrics from the two schemes, and to gain physical insight into their divergence, we compared power and additively weighted Voronoi diagrams of polydisperse ensembles representing physically relevant ranges of polydispersity, density, and overlap. When tested against experimental distributions of gas foams, the results related their divergent power and additively weighted decompositions to the polydispersity of their particle size distributions. Geometric analysis of the Voronoi cells implicated the subpopulation of small objects as the primary contributors to the divergence through their preferential assignment of larger, aspherical power cells relative to their additively weighted counterparts.

---

Introduction

The spatial description of disperse matter is of fundamental interest to the physical and engineering sciences. At nanometer length scales and above, suspensions of materials or entrapped bubbles constitute dispersed phases whose spatial and geometric properties reflect their thermodynamic stability and other colloidal characteristics.^1–4 A physically appropriate description of the volume and surface area of the dispersed objects relative to their embedding medium is essential to understanding and predicting their bulk properties and dynamics such as phase separation and response to mechanical stress.

A standard approach to structuring disperse matter is the Voronoi diagram.⁵ The primitive Voronoi diagram divides the space around a set of points (called generators) into cells such that each cell contains all the space closer to a generator than any other generator. Thus, for a set of generators S = {p, q, …}, the Voronoi cell V associated with p is the locus of points x with the shortest distances to p as follows:

V(p) = {x|d(x,p) ≤ d(x,q) ∀ q ∈ S, q ≠ p} (1)

where d(u,w) ≡ ‖u − w‖ is the Euclidean norm. Other formulations of the distance function (e.g., L^P norms) are also possible. As Voronoi diagrams of points do not consider the sizes of the generators, application to real structures leads to unphysical decompositions in which larger non-overlapping objects protrude out of their assigned Voronoi cells [Fig. 1A]. For spherical objects (balls) s with center c and radius r, the radius is incorporated into various weighted distances that replace the Euclidean norm in eqn (1). Weighting schemes differ in features that make them more or less suited for particular physicochemical applications. For example, multiplicative weighting i.e., d(x,p)/r admit spaces that are not assigned to an object.⁶ Among the space-filling schemes, the additively weighted Voronoi diagram and the power diagram (also known as radical and Laguerre tessellation) are the most commonly used. Power diagrams are defined by the following weighted distance known as the power distance π:⁷

π(x,s) = d²(x,c) − r². (2)

Fig. 1 Illustrative primitive, power, and weighted Voronoi decompositions in the plane. (A) Primitive and weighted dividing lines for two dissimilar disks determined by assignment of randomly placing points (top) or analytically (bottom) according to eqn (1)–(3). In primitive Voronoi diagram, the dividing line is equidistant (green segments) between the centers of the generating objects without regard to size. Unphysical decompositions occur when the dividing line intersects non-overlapping balls, which are avoided by power (P) and additively weighted (AW) diagrams. In the AW scheme, the dividing line is equidistant to the surface of the generators. (B) A network dividing the space housing five disks creating a cell around the center disk. The curved edges in an additively weighted diagram also govern the partitioning of overlap between unequally sized disks. The intersections of edges represent vertices of the Voronoi cells.

In practice, power diagrams are attractive for the speed and convenience of computing flat-faced cells. In contrast, the additively weighted diagram, which is also known as the Voronoi diagram of spheres or the Apollonius diagram, divides space such that each cell contains the points closer to the surface of ball s than any other by way of the weighted distance metric:⁵

δ(x,s) = d(x,c) − r. (3)

Whereas the curved faces and edges resulting from additive weighting have a more direct physical correspondence to the radius, the power distance π in eqn (2) is less straightforwardly interpreted as the capability of s to “influence its neighborhood and claim a larger cell.”⁵

Voronoi diagrams are used to model molecular^8,9 and ionic liquids,¹⁰ colloidal crystals,^11,12 granular materials,^13,14 foams/bubbles,^15,16 and myriad other physicochemical systems. Such systems are commonly considered with Voronoi diagrams of balls because the objects themselves are spherical or can be modeled as assemblages of spheres.^17,18 Molecular systems, whose atomic constituents are approximated as balls in the classical limit, are also most correctly modeled by additively weighted Voronoi diagrams.¹⁹ Compared with power diagrams, additive weighting is considerably complex and computationally intensive,^20–22 especially for balls that overlap [Fig. 1B]. For many years, available code for additively weighted Voronoi diagrams, whether proprietary (e.g., MATLAB, MATHEMATICA) or open source (CGAL library), was limited to the plane. Freely accessible software for three-dimensional work has only recently become available.^23–25 Outside of biomolecular applications,^25–28 power analysis of physicochemical systems (particular continuum matter) has remained the norm.^12,29–35 Given this, and the significant body of literature on weighted Voronoi diagrams, a practical question arises as to how the two weighting schemes differ in their description of physicochemical systems. Although the theoretical linkage among weighting schemes has been made,³⁶ abstract analysis does not readily translate to physical interpretations.

As the only guaranteed condition for convergence of power and additive weighting is radial uniformity of the generators, deviation between the two schemes is the default expectation. Quantitative calibration of weighted diagrams and physical insight into their divergence are useful but remain open knowledge gaps. To make progress in a broadly applicable way, we compared power and additively weighted Voronoi descriptions of spherical ensembles commonly used to model physicochemical systems. These results: (1) provide insight into the fundamental factors that drive divergence between the two weighting schemes, (2) establish a quantitative basis for interpreting volumetric metrics of power diagrams under the additively weighted régime, and (3) extend previous studies on oligo-disperse systems.^37–44 In support of these efforts, we share with the community open-source, python-based software for the integrated computation and analysis of weighted Voronoi diagrams from user-supplied generator coordinates and radii.

Theory and methods

Mathematical background

The Voronoi diagram contains all the space closer to the generator ball than any other ball. The face of a Voronoi cell is therefore the subset of points equidistant (bisector) from the generator ball and the adjacent ball. Equivalently, each point on the bisector is the center of an sphere which is tangent to both balls. The bisecting surface is a sheet of hyperboloid when the balls differ in radius and a plane when they have the same radius (Fig. 1A). As each face is a subset of the center points of balls tangential to two balls, the intersection of two faces must be a subset of the center points of balls tangential to three balls. Furthermore, the intersection of three edges on the Voronoi cell is a subset of the center points of balls tangential to four balls.

For an additively weighted decomposition space, we can re-write eqn (3) in coordinate form as the distance r from a point (x, y, z) to the surface of a ball s_i with center c = (x_i, y_i, z_i) and radius r_i as:

where the indexed coordinates refer to the center of the ball s_i of radius r_i. Setting r_i = 0 returns the primitive (unweighted) Voronoi diagram. For the corresponding power distance, cf.eqn (2),

(x − x_i)² + (y − y_i)² + (z − z_i)² − r_i² = r (5)

The properties of power diagrams have been extensively described by Aurenhammer.⁷ For additively weighted diagrams, eqn (4) can be considered equivalently as a ball with radius r that is tangent to the ball s_i. The surfaces (s), edges (e) and vertices (v) can be found by rearranging eqn (4) and solving a system of equations between two, three, and four balls respectively:

When r_i ≠ r_j, the surface is a hyperboloid of two sheets, only one of which is the true bisector of the two balls. This property leads to complications when computing the Voronoi edges and the Voronoi vertices because intersections of artefactual sheets will also be solutions to the above systems of equations. Nevertheless, it is possible to recover the point on the true Voronoi cell by finding the weighted distance from the point to each participating ball. The point on the true surface of the Voronoi cell will have the same weighted distance to each of the balls.²²

Construction of statistical ensembles

Statistical ensembles of balls that model random packing were generated using an in-house Python code, FOAM_GEN, based on user-specified distribution type, mean ball radius, density, number of balls, polydispersity, overlap allowance, and periodic boundary condition (PBC). The ensembles are constructed by randomly packing balls with radii sampled from the target statistical distribution in a cubic bounding box whose dimensions correspond to the specified density. Implementation details, particularly the handling of ball overlap in predicting the bounding box dimensions by a machine learning model,⁴⁵ are provided in Supplementary methods (ESI†). Completed ensembles are exported in Protein Data Bank (PDB) format for analysis.

Computation of Voronoi diagrams

Many computationally efficient (hardware-accelerated) techniques^{22,23,46–49} for constructing weighted ordinary Voronoi diagrams of spheres and other parametric geometries such as ellipses,⁵⁰ as well as the more intensive task of general Voronoi diagrams of arbitrary objects^48,51–53 have been reported. Given our priorities for the analytical solution of the Voronoi network and surface triangulation, we based our implementation of the additively weighted Voronoi diagram, named VORPY and coded in Python, on the edge-tracing Voronoi S-network reported by Medvedev and coworkers.⁵⁴ Algorithmic improvements to the original descriptions and their implementations will be the subject of a separate work. VORPY takes as input PDB format with radius replacing the B-factor column. In addition to additively weighted Voronoi diagrams, VORPY computes power and unweighted diagrams for direct comparison with their additively weighted counterparts. Diagrams consist of triangular meshes bounded by polygonal lines for neighboring balls. To calculate the surface area for any given cell, the areas of the triangles in each mesh are summed to the total. The volume of the cell is calculated by summing the tetrahedrons formed by the same triangles and the center of the ball assigned to a cell. VORPY outputs a comprehensive suite of information items, including a log file describing the network (vertices and their linkages to edges and faces) and their visualization in PyMOL. Volumes, surface areas, and other geometric properties are also printed in the output log.

Results and discussion

Effect of size difference on weighted Voronoi diagrams

For an intuitive evaluation of how unequal size differentially perturbs additively weighted and power diagrams in space, we first considered a toy system in which a probe ball is embedded within a fixed lattice of 52 uniform balls [Fig. 2A]. Size difference is rendered through the radius of the embedding balls only, without changing the coordinates of any ball, over a range in which the radius of the probe is 0.45 to 1.38-fold of the embedding balls. The probe was inserted without regard to the coordinates of the embedding balls and overlapped with one or more embedding balls at all radial ratios tested. Overlap among the embedding balls was present only at radial ratios above unity.

Fig. 2 Illustrative comparison of additively weighted Voronoi and power diagrams. (A) A toy model in which a probe ball (blue) is embedded within a fixed lattice of 52 uniform balls of radius scaled with the probe as indicated. The colored cells represent the additively weighted cells for the probe ball. White cells with straight edges and flat surfaces represent the power cell for the same probe. The probe ball envelopes the cell entirely at the highest radial ratio and is not shown. The rendering artifact at unit radial ratio reflects the coincidence of the weighted cells when the probe and embedded cells are identical in radius. (B) The scaling of weighted cell properties with the radial ratio. (C) The volume and surface area ratios of the additively weighted and power cells as a function of radial ratio.

Across the range of radial ratio tested, both the additively weighted and power cell decreased in volume and surface area as the probe became small relative to the embedding balls [Fig. 2B]. The edges of the additively weighted cell transitioned from being concave to flat, to convex as the radial ratio increased past unity. The curved edges are convex for additively weighted cells generated by smaller balls. Relative to power cells, additively weighted cells showed the same number of facets except at the lowest radial ratio, where the sculpted additively weighted cell showed additional small facets absent in the power cell. At high radial ratios, the probe cell became less featured for both weighting schemes. For this model, the additively weighted probe cell was larger than the power cell above and below unit radius [Fig. 2C]. The convergence of the cells as the size tends to uniformity is apparent at unit radial ratio.

Generation of statistical ensembles for Voronoi decomposition

To meet its illustrative purpose, the toy model is a bidisperse lattice that is controlled through a single parameter (the radial ratio) to uniformly inflate or deflate the probe-embedding balls at fixed positions. Recognizing the size and positional heterogeneity of physical systems that are widely modeled as spheres, we considered statistical ensembles used to model physicochemical systems and whose structures have been analyzed with power diagrams.^12,29–34 Experimental studies of diverse physicochemical systems including gas-based foams,⁵⁵ ceramics,⁵⁶ nano-sized polycrystals,⁵⁷ and powdered food⁵⁸ have established size distributions that typically fall into gamma and lognormal distributions. Given their similarity³³ and theoretical preference⁵⁹ for gamma over lognormal distributions for crystalline systems, we focused on the gamma distribution as a general framework for analyzing weighted Voronoi decompositions of polydisperse balls.

To proceed, we simulated ensembles consisting of 1000 balls with radii r sampled from the gamma distribution at specified level of polydispersity, density, and overlap [Fig. 3A]. The gamma distribution is a two-parameter function:

where . To incorporate polydispersity explicitly into the distribution, we re-parameterized eqn (7) in terms of the coefficient of variation CV ≡ σ/μ (σ, standard deviation; μ, mean):where k = θ/μ and θ = μ·CV². For simplicity, we set μ = 1 since Voronoi diagrams are scale-invariant.

Fig. 3 Statistical ensembles of balls: construction and characteristics. (A) Outline of the ensemble construction process as implemented in FOAM_GEN. See Theory and methods for details. (B) Illustrative ball packing in the absence and presence of periodic boundary condition (PBC). With PBC, balls are allowed to protrude from the retaining box and are checked against the balls in the corresponding adjacent unit(s) for overlap. (C) Density adjustment for overlapping systems. Given the set parameters (number of balls, density, CV, PBC, distribution, overlap), a linear model was trained to predict the required adjustment on the nominal non-overlapping density (see ESI†). Twelve illustrative realizations of the model are shown (number of balls, overlap, CV): 1 – (100, 1.5, 0.5), 2 – (100, 1.5, 1.0), 3 – (1000, 1.5, 0.5), 4 – (1000, 1.5, 1.0), 5 – (100, 1.0, 0.5), 6 – (100, 1.0, 1.0), 7 – (1000, 1.0, 0.5), 8 – (1000, 1.0, 1.0), 9 – (100, 0.5, 0.5), 10 – (100, 0.5, 1.0), 11 – (1000, 0.5, 0.5), 12 – (1000, 0.5, 1.0). (D) Representative ensembles across density, CV, and overlap. (E)–(G) Representative positional, radial, and overlap distributions of PBC ensembles. Distributions for the full dataset are given in Fig. S1, ESI.† (E) x, y, and z relative locations of the balls within the bounding cube at (CV, density) = (0.5, 0.25) with and without overlap. The “static” appearance reflects the homogeneity of the locations across systems. (F) Histograms of ball radii to the gamma distribution sampled at the indicated CV values and density = 0.25. (G) Overlap (relative to the smaller ball) of PBC ensembles.

Given a distribution of radii, a retaining box was constructed for packing the balls given the total volume of the set of balls and the desired density. To assess the impact of boundary effects, ensembles were generated with and without periodic boundary condition (PBC) to approximate continuous domains or as an isolated system [Fig. 3B]. In general, the total volume of an ensemble is the summed volume of the constituent balls (the non-overlapped volume) minus any volume involved in overlap. For definitional purposes, we expressed overlap between two balls as the maximum fractional penetration into the smaller ball. For example, overlap = 0.5 means that the larger of any two neighboring balls was allowed to penetrate up to the center of the smaller ball. To account for reduction in density of the ensemble due to overlap, a machine learning model was trained to predict the density adjustment given the set parameters (number of balls, density, CV, PBC, distribution, overlap) [Fig. 3C]. Once sorted from largest to smallest, the balls were randomly placed in the retaining box and checked for permitted overlap with the other balls. Completed ensembles were then analyzed by Voronoi decomposition (vide infra).

Ensembles consisting of 1000 balls covering a wide range of polydispersity (CV from 0.05 to 1.0), density (0.05 to 0.5), and overlap (0.0 to 1.0) were generated. A packing density (fraction) up to 0.5, a range representative of solid⁶⁰ and liquid foams,⁶¹ was sampled. For each composition, twenty independent replicates were realized. Representative ensembles are shown in Fig. 3D. As quality checks, we examined the ensembles for positional homogeneity [Fig. 3E], fidelity to the underlying gamma distribution of size [Fig. 3F], and the prevalence of overlapping balls [Fig. 3G]. The number fraction of overlapping balls scales sharply with density for ensembles with low polydispersity, in accord with the increasing restrictiveness with which similarly sized balls may be packed without overlap by the static packing scheme implemented in FOAM_GEN.

Comparison of power and additively weighted Voronoi decompositions of statistical ensembles

Our implementation of the additively weighted Voronoi diagram (VORPY) is based on an edge-tracing approach.⁵⁴ To test VORPY for correctness, we compared the additively weighted decomposition of 21 ensembles (an independent test set representing the extremes of densities, CV, and overlap of the foregoing ensembles), with VORONOTA,²³ a solver that implements an alternative “gift-wrapping” strategy. The two implementations share >98% matched vertices (over 6 to 9 × 10³ per ensemble), which are close to within 0.5% in position [Fig. S2, ESI†]. These results support the correctness of VORPY for additively weighted diagrams.

For each ensemble, a power and additively weighted decomposition was computed. Since both types of weighted cells map 1 : 1 to the generator set of balls, the assignment of Voronoi cells should be complete in principle under PBC. In practice, a small fraction of cells at very high CV is not assigned due to vertices located beyond a preset cutoff from the retaining box. For non-PBC ensembles, the outer layer of balls form unbounded Voronoi cells which are removed from analysis. The completeness of weighted Voronoi cells due to these factors are shown in Fig. S3 (ESI†).

For each ball with an assigned power and additively weighted cell, the properties of the two cells were compared over the full ensemble. Since deviations in volume or surface area between neighboring cells mutually compensate, the average difference over all balls would be zero under PBC. For non-PBC ensembles, the interior balls would have no effect on the final value and the only remaining difference would come from the exterior balls partitioned values. Given these features, we define as summary metric the average absolute differences in cell volume:

where the subscripts “P” and “AW” denote the power and additively weighted cell of the i-th ball. An average absolute difference in surface area, avg(|Δ_a|) was defined analogously.

Inspecting avg(|Δ_v|) and avg(|Δ_a|) as a function of density over a range of CV values shows that, controlling for density, polydispersity increases the deviation between corresponding power and additively weighted Voronoi diagrams. At sufficient polydispersity, even low-density ensembles deviate significantly in volume between power and additively weighted cells e.g., avg(|Δ_v|) > 50% at CV = 0.6. As the density dependence of avg(|Δ_a|) is qualitatively identical to avg(|Δ_v|) except smaller in magnitude [Fig. S4, ESI†], we will henceforth focus on avg(|Δ_v|). The divergent effect of CV is manifest in both the overlapping and non-overlapping systems, with or without PBC [Fig. 4A]. Non-PBC ensembles exhibit different avg(|Δ_v|)–density profiles than PBC ensembles, suggesting finite-boundary effects, but nevertheless follow the same trend with polydispersity at equivalent density. Parenthetically, a comparison set of 100-ball ensembles shows identical trends with wider variances [Fig. S5, ESI†], indicating that finite-boundary effects do not significantly modulate the intrinsic statistics of the weighted Voronoi decompositions and that the 1000-ball ensembles balance converged statistics and computational efficiency.

Fig. 4 Divergence of power and additively weighted Voronoi decomposition of statistical ensembles. (A) Deviations between power and additively weighted volumes as measured between equivalent balls across three different overlap allowances, 11 different CV values, 10 different densities, across periodic (top row) and non-periodic (bottom) ensembles of balls. The shaded bands indicate ±1 SD of each set of 20 replicates. (B) Illustrative ensembles sampled at the indicated configurations marked in panel (A). Balls with unbounded cells from non-PBC ensembles are rendered semi-transparent. Balls are colored according to a map of avg(Δ_v) on the same scale across the ensembles. Balls whose power volumes are greater than their additively weighted volumes appear red and balls whose power volume is less than their additively weighted volume appear blue. Overlapped ensembles are labeled in magenta. The scale for avg(Δ_v) in the ensembles shown in 9* and 10* are scaled up by 10× to render the narrower deviations more apparent. (C) Distribution of signed deviation between power and additively weighted cells over ball radius. The average radius is set at unity. Curves represent hyperbolic fits to the points intended only to guide the eye. (D) Representative balls and power (P) and additively weighted (AW) cells are red and purple respectively from the three size régimes relative to the average radius (set to unity). Zoom levels are normalized along columns to balance visualization of geometric details and roughly equal sizing of the additively weighted cells across rows.

In contrast to polydispersity, density (ρ) effects on the deviation between corresponding power and additively weighted cells are more complex. At high CV values e.g., CV > 0.8 under PBC, avg(|Δ_v|) increases monotonically with density within the range tested. Although packing of highly polydisperse balls increases the divergence in volume between power and additively weighted cells, the positive slope moderates with density. At lower polydispersity e.g., CV < 0.4, the slope with density is negative . At intermediate polydispersity, a reversal in the density-slope of avg(|Δ_v|) is observed. The second derivative of avg(|Δ_v|) as a function of density is, therefore, overall negative i.e., , and indicates compensation between a volume-divergent effect by polydispersity on the one hand, and convergent effects on avg(|Δ_v|) by packing on the other.

Some features of the avg(|Δ_v|)–density profiles reflect expected properties of power and additively weighted Voronoi cells. As elaborated in the Appendix, the power cells defined by a pair of unequally sized, non-overlapping balls intersect the ball–ball axis closer to the larger ball than the corresponding additively weighted cells. Moreover, the flat-faced power cell over-captures the volume for smaller balls than the curved additively weighted cells, and under-capture volume for larger balls. The reverse is true for a pair of overlapping balls. The translation to polydisperse ensembles is illustrated in Fig. 4B in which balls are colored by the (signed) Δ_v of their assigned weighted Voronoi cells. Mapped to a common scale, highly polydisperse ensembles exhibit a wider range in avg(Δ_v). Small balls are associated with positive avg(Δ_v) i.e., inflated power cells while larger balls are associated with negative avg(Δ_v). The absence of the inflated power cells in the toy model (Fig. 2) at low radial ratio (i.e., small probe), where the embedding balls are uniform, underscores polydispersity as the primary driver of divergence between the two schemes.

In the presence of ball overlap, the ensembles at equivalent polydispersity and density are progressively less volume-divergent i.e., lower avg(Δ_v) than their non-overlapping counterparts. This may be rationalized by the curved additively weighted faces between unequal neighbors, which are concave with respect to the smaller ball (cf., Fig. 1B). Within an overlap, this concavity increases the share for the smaller ball (versus a power cell) in opposition with the effect without overlap. Thus, overlap uniformly opposes the volume-divergent effect of polydispersity.

With respect to density, it may be rationalized by noting that power and additively weighted edges are locally the least divergent where neighbors are the closest (cf., Fig. 1B). This would suggest that closer packing (higher density) leads to progressively less-curved additively weighted cells. Such an explanation is incomplete, however, given that the vertices of power and additively weighted cells do not generally coincide i.e., power cells are not flat-faced additively weighted cells, and that for polydisperse ensembles at high density. Polydispersity thus divides the density dependence of avg(|Δ_v|) into two distinct régimes.

To dissect the non-uniform effects of polydispersity and packing on power and additively weighted cells, we analyzed avg(Δ_v) as a function of radius [Fig. 4C]. By definition, the sign of avg(Δ_v) for any ensemble is balanced at the radial average which was set at unity. With increasing density, balls below the average radius at low polydispersity differentially inflate power cells, while balls above the average radius are differentially deflationary. This behavior is reinforced by overlap. By contrast, there is no reversal with density in the role of balls above and below the average at high polydispersity. Instead, balls below the average radius strongly inflate power cells with increasing density. This contribution is moderated but not mitigated by overlap. The analysis thus specifically implicates smaller balls as the drivers of avg(Δ_v) and avg(|Δ_v|) via their strong inflation of their power cells at high polydispersity, several examples of which are shown in Fig. 4D.

Application to experimental ensembles

Using the divergence between power and additive weighted volumes of individual cells and their distributions across the ensemble as a framework, we examined the representative experimental distributions such as reported by DeVries,⁶² Gal-Or and Hoelscher (GO–H),⁶³ and Ranadive and Lemlich (R–L)⁶⁴ for aerated gas bubbles. These models had also been evaluated for compressed air foams by Magrabi and coworkers.⁵⁵ The probability distribution functions (PDFs) are:

Although only the R–L distribution is formally a gamma distribution, the three models are similar in line shape [Fig. 5A].

Fig. 5 Weighted Voronoi diagrams of experimental foam distributions. Three different experimental distributions for gas-based foams were reported by DeVries (purple), Gal-Or and Hoelsher (GO–H, green), and Ranadive and Lemlich (R–L, yellow). (A) Probability distribution functions (PDFs) for the three distributions and histograms of randomly sampled radii. (B) The absolute % difference in power volume from additively weighted volume. (C) Distributions of Δ_v at high-density, 0.5 by radius for the three distributions. Curves represent hyperbolic fits to the points intended only to guide the eye. The dotted lines show the avg(|Δ_v|) for each ensemble.

We generated ensembles of balls using radii sampled from the three experimental distributions and analyzed the resultant density dependence of avg(|Δ_v|). The parameterization by DeVries showed larger deviation between the power and additively weighted decompositions over the R–L and GO–H distributions [Fig. 5B]. Based on the developed framework, the greater magnitude of and density dependence of avg(|Δ_v|) for the DeVries distribution is expected to have its basis in a higher polydispersity over the other two distributions. This hypothesis is borne out by a direct estimate of the CV for the three distributions, which is largest for the DeVries distribution (CV = 0.600). The DeVries avg(|Δ_v|)–density profile compares most closely with the non-overlapping ensembles drawn gamma distribution sampled at CV = 0.7 (Fig. 4A), supporting the general applicability of the gamma distribution to physicochemical systems.

Both the R–L and GO–H distributions showed very similar avg(|Δ_v|) over the entire range of density tested, even though their respective PDFs are not themselves superimposable. The R–L distribution exhibits a wider skew with a measured standard deviation of 0.422 (versus 0.317 for GO–H). However, the GO–H distribution is more skewed towards small balls with a mean radius of 0.751. The two statistics compensate to yield an identical CV value of 0.422. The superimposable avg(|Δ_v|) across the full density range tested, thus provides further and orthogonal support for the foregoing results as a framework for interpreting experimental data.

Packing moderates power-cell inflation by promoting more spherical Voronoi cells

Inspection of weighted Voronoi cells such as those in Fig. 4D, suggested that inflated power cells assigned to small balls exhibit elongated or otherwise non-spherical morphology compared with their additively weighted counterparts. This observation prompted us to consider non-spherical morphology as a factor on the divergence between power and additively weighted cells as indicated by avg(|Δ_v|). We therefore examined the combined results from the most polydisperse ensembles (CV = 1.0) at low and high densities, with and without overlap, and quantified their weighted Voronoi cells in terms of sphericity Ψ as follows:where v and a are the cell volume and surface area, respectively. By this definition, Ψ = (0, 1], with Ψ = 1 corresponding to a perfect sphere. To compare the sphericity of corresponding power and additively weighted cells, we define:

ΔΨ ≡ Ψ_P − Ψ_AW (12)

A ball whose power cell is more spherical than its additively weighted counterpart has a positive ΔΨ value and vice versa. The results are presented as two-dimensional histograms (biplots) of avg(Δ_v) and ΔΨ (Fig. 6A). On the one hand, there is a clear correlation between less spherical power cells (negative ΔΨ) with inflated power volumes (positive Δ_v; upper left quadrant). Indeed, the sign of Δ_v also divides the subpopulations of inflated power cells by their relative non-sphericity (left and right halves of the biplots). On the other hand, ΔΨ is strongly determined by ball radius, with the power cells of small balls overwhelmingly exhibiting less spherical geometry. Comparison of the distributions by overlap shows that the bias among the four quadrants is less pronounced when overlap occurs. Taken together, the results implicate small balls in producing less spherical and more volume-inflated power cells that dominate the overall divergence in avg(|Δ_v|) relative to additively weighted Voronoi cells.

Fig. 6 Asphericity of power cells drives volume inflation and divergence from additive weighting. (A) and (B) Representative distributions of Δ_v and sphericity difference (ΔΨ) between power and additively weighted cells across PBC and non-PBC ensembles. Each plot on the left shows the pooled ensembles corresponding to the marked position in the avg(|Δ_v|)–density profile, reproduced from Fig. 4A, on the right. (A) Plots of the PBC sphericity difference of each cell from ensembles with CV = 1.0 across overlap values of 0.0 and 1.0 and density values of 0.1 and 0.5. (B) Corresponding distributions from the non-PBC ensembles with CV = 1.0 and densities of 0.1, 0.4, and 0.5 and overlap values of 0.0 and 1.0.

To probe the generality of the interplay between sphericity and distribution of ball radius further, we scrutinized the non-PBC ensembles under the same conditions of high polydispersity (CV = 1.0) because the case with overlap = 1.0 presents a clear-cut example of the volume-convergent effects of packing [Fig. 6B]. We scrutinized ensembles sampled at extreme dispersities as well as the maximum divergence (marked 14, 15, and 16 in Fig. 4A). As a “control”, the non-PBC ensembles without overlap (points 11, 12, and 13) recapitulate the monotonic trends noted in their PBC counterparts in Fig. 5A. However, the reversal in avg(|Δ_v|) in the case of overlap (points 14 to 16) distinctly tracks with the relative abundance of the subpopulation in the upper left quadrant i.e., aspherical, volume-inflated power cells. The statistics thus strongly complicate the less spherical power cells of small balls as the primary contributor to the volumetric divergence from additively weighted cells.

Rationale for the Voronoi cell sphericity and volume divergence

To gain insight into the mechanism underlying cell sphericity and volume divergence, we considered the local geometry of weighted Voronoi cells. For a pair of non-overlapping balls, radial disparity uniformly places the power boundary closer to the larger ball along the axis defined by the ball centers (Appendix). This property translates, in the context of whole cells, to the two weighted schemes defining neighbors that are common (shared) or exclusive to one or the other scheme. Exclusive neighbors profoundly influence the sphericity of the weighted cell, as illustrated by two simple systems in Fig. 7A. A small ball with more exclusive power neighbors seeds a less spherical power cell than additive weighting. In contrast, a larger seed ball with fewer exclusive power neighbors seeds a more spherical power cell than additive weighting. Further translating to ensembles, Fig. 7B shows the conservation of excess exclusive neighbors for the power cells of small balls and excess exclusive neighbors for the additively weighted cells of larger balls. Crucially, exclusive neighbors are further away from the seed balls than shared neighbors [Fig. 7C]. Integrating these features, more exclusive neighbors drive less spherical cells, for both weighting schemes, as a function of polydispersity and density [Fig. 7D]. The association with more exclusive power neighbors thus explains the less spherical and larger power cells of small balls in Fig. 6. More generally, deviation in sphericity acts as a metric of the local structure around a ball within the ensemble.

Fig. 7 The geometry of power cell asphericity and volume inflation. (A) Illustrative 2-D example of a small (left) and large (right) disk embedded among disks with additively weighted and power boundaries (solid) and cutoff boundaries (dashed) by the closer disks. (B) Number of exclusive neighbors by radius. Each ball is represented twice, once for its power cell and another for its additively weighted cell. (C) Average distance of shared and exclusive neighbors, segregated by exclusivity to one or the other weighted cell, as a function of ball radius. The ensembles correspond to those in panel (B). (D) Correlation between difference in sphericity (ΔΨ) and neighbor number difference. The ensembles correspond to those in panel (B). (E) 2-D diagram of overlap showing the area contributions from regions defined by the additively weighted and power boundaries. As one moves away from center axis (positive D), the contribution from weighting to the smaller disk increases. (F) Construction of the local sector volume, projected from a fixed test point, for the additively weighted (red) and power (blue) cells of two overlapping balls. (G) Difference in local (sector) volume as a function of deviation from the center axis (D) for two ratios (1 : 2, 1 : 10) of different radii and three overlap values. (H) Absolute volume difference avg(|Δ_v|) by cell as a function of neighbor numbers at three overlap allowances, colored by % overlap. Histograms indicate the number distribution of neighbors in each ensemble. The ensembles have values of CV = 1.0 and density 0.5.

If balls intersect, the partition within the overlap produces a contribution that favors the additively weighted cell of the small ball (cf., Fig. 2B), in opposition to the situation outside the overlap or in the disjoint case. As derived in the Appendix and illustrated in Fig. 7E for a pair of overlapping balls, this contribution arises from a reversal of the relative positions of the power and additively weighted boundary along the center axis. To calibrate the compensation afforded by the overlap against the inflated power cell of the smaller ball, we analyzed the local volumes of the sectors enclosed by the projection from a test point on the center axis to its power and additively weighted boundaries as a function of deviation (denoted by D) from the axis [Fig. 7F]. This analysis models the compensatory effect of density on avg(Δ_v) in overlapping ensembles. With increasing distance from the center axis, the loss of volume compensation is primarily sensitive to the extent of overlap, but not difference in ball sizes [Fig. 7G]. The translation to ensembles [Fig. 7H] recapitulates this mechanism by which the overlapping balls comprise the least volume-divergent constituents across the full range of local packing scenarios.

Conclusion

The studies reported here indicate that the divergence between a power and an additively weighted Voronoi description of a statistical ensemble is the balance between polydispersity (CV) on the one hand and packing (density and overlap) on the other. Indeed, at sufficiently low CV, dense packing leads to a convergence under the two schemes. The underlying mechanism is traced to the strong propensity of small balls (relative to the average radius) to generate power cells that are less spherical than, and more volumetrically divergent from, their additively weighted counterparts. The generality of this behavior is evidenced in the fully overlapping non-PBC ensembles. In these ensembles, the center of every ball must be at least one radial distance from the boundary, resulting in preferential placement of small balls around the boundary and removal from the final ensembles as their Voronoi cells are incomplete. In this biased population, the volume-convergent effects of packing and overlap overcome the volume-divergent effect of polydispersity (Fig. 6B). On this basis, we expect that PBC ensembles at high CV will exhibit volume convergence at higher packing density than our static algorithm affords. Volume fractions up to ∼0.9 are achievable using dynamic and/or constructive packing algorithms that simulate polydisperse powders and other granular materials.^65–69 Similarly, one expects Voronoi descriptions of ensembles modeled by physics-based schemes, such as Plateau's laws (surface minimization) for soap films,⁷⁰ to also be covered under this framework.

From a practical point of view, significant rationale exists for colloidal suspensions and engineering foams to enrich fine particles (e.g., by sonication, impelling, extrusion, sparging) for enhanced mass and thermal transfer with their embedding phase.^29,71,72 Such physically important ensembles would be expected to exhibit, at equivalent density, more divergent power approximation of their additively weighted Voronoi volumes. Conversely, the coarsening or ripening of foams and other colloidal systems, where small constituents are preferentially removed from the population, is expected to favor volume convergence. In the literature, hybrid weighting schemes have been proposed that define regions in which additively or power weighting are applied to reduce computational burden,⁷³ or evolve physics-based curvature (e.g., Plateau's laws) in a power framework.⁷⁴ As an alternative, the present results provide a convenient calibration from one weighting scheme to the other over a broad range of polydispersity, density, and overlap when direct computation of the underlying diagram is not needed (or possible). Usefully, as applications to experimental foam systems show, these calibrations give consistent predictions and physical interpretations even for distributions that are not gamma in form.

As part of this work, we developed fully open-source code (VORPY) for constructing and analyzing power and additively weighted Voronoi diagrams. Current programs that compute weighted Voronoi diagrams have catered to biomolecular applications^23,25 in which the constituent balls (atoms) have radii discretized by composition and distance and angular constrained by chemical bonding. The atomic density (or packing coefficient) in protein crystals is typically between 0.5 and 0.6,⁷⁵ and therefore similar in magnitude to the overlapping ensembles investigated here. However, the balls in statistical ensembles may have any radius and geometry permitted by the preset size distribution and overlap. The permissive size specification and organization of balls in statistical ensembles present analytically challenging geometries for additive weighting not likely encountered with structured macromolecules such as proteins. The solutions to these scenarios, which are detailed in a separate work, commend VORPY for a wide scope of applications. VORPY therefore complements the field as an integrated tool for analyzing physicochemical systems by weighted Voronoi diagrams.

Author contributions

J. M. E.: data curation, formal analysis, investigation, software development, validation, visualization; writing – original draft, review & editing; N. W.: methodology, supervision, writing – review & editing; G. M. K. P.: conceptualization, project administration, resources, visualization, supervision, writing – original draft, review & editing, funding acquisition.

Conflicts of interest

The authors declare no competing financial interest.

Data availability

Python-based software for generating statistical ensembles (FOAM_GEN) and analyzing their weighted Voronoi decompositions (VORPY) can be obtained on github as follows: FOAM_GEN (1.0.1) https://github.com/jackericson98/foam_gen.git. VORPY (3.0.5) https://github.com/jackericson98/vorpy.git. The datasets for this article, consisting of ensembles generated by FOAM_GEN and their weighted Voronoi analysis by VORPY, are available at Zenodo: https://doi.org/10.5281/zenodo.14902130.

Appendices

Relative position of the power and additively weighted Voronoi boundary in the locale of two balls

Consider two balls S₁ and S₂. Without loss of generality, let the center of S₁ be at the origin, and the center of S₂ lie on the positive x-axis at x = a > 0 i.e., S₁((0, 0, 0), r₁) and S₂((a, 0, 0), r₂):

We further assume that r₂ ≥ r₁, and only consider pairs of balls S₁ and S₂ such that S₁ is not contained in S₂. This implies:

r₁ ≤ r₂ < a + r₁ (13)

We aim to determine where the Voronoi surface of power and additively weighted Voronoi diagrams (awVD) between S₁ and S₂ intersects the x-axis.

Power	awVD
	|xaw| − r1 = |xaw − a| − r2

To solve the absolute value equation for the awVD with respect to x, we need to distinguish three cases:

(1) x_aw < 0: in this case, from |x_aw| − r₁ = |x_aw − a| − r₂, we get −x_aw − r₁ = −x_aw + a − r₂, which leads to r₂ = a + r₁. This case was excluded by assumption.

(2) x_aw > 0: in this case, from |x_aw| − r₁ = |x_aw − a| − r₂, we obtain x_aw − r₁ = x_aw − a − r₂, which leads to a = r₁ − r₂. Since, by assumption, a > 0 and r₂ ≥ r₁, this case cannot occur either.

(3) 0 ≤ x_aw ≤ a: in this case, from |x_aw| − r₁ = |x_aw − a| − r₂, we obtain x_aw − r₁ = a − x_aw − r₂ and therefore .

Thus, the intersection points are:

Power	awVD

Since r₂ ≥ r₁, the correction terms which are subtracted from are always non-negative. If r₁ = r₂, the intersection point is at in both diagrams. As r₂ increases, the intersection shifts toward the origin and thus toward S₁. This increases the contribution of S₂ to the Voronoi cell of S₁.

From the formulae, one can also deduce when the face from the power diagram or the awVD lies closer to S₁:

(1) If S₁ and S₂ are disjoint, then r₁ + r₂ < a, and thus . In this case, the face from the awVD lies closer to S₁ than the face from the power diagram.

(2) If S₁ and S₂ overlap, then r₁ + r₂ ≥ a, and thus . In this case, the face from the power diagram lies closer to S₁ than the face from the awVD.

Acknowledgements

We thank Grant C. Derdeyn-Blackwell and April M. Cole for preliminary analysis as well as Drs W. David Wilson and Markus G. Germann for helpful suggestions. This work is supported by NSF grant MCB 2028902 to G. M. K. P.

References

1. N. Galvani, M. Pasquet, A. Mukherjee, A. Requier, S. Cohen-Addad, O. Pitois, R. Hohler, E. Rio, A. Salonen, D. J. Durian and D. Langevin, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2306551120.
2. I. Cantat, S. Cohen-Addad, F. Elias, F. Graner, R. Höhler, O. Pitois, F. Rouyer, A. Saint-Jalmes, R. Flatman and S. Cox, Foams, Oxford University Press, 2013.
3. D. Langevin, Emulsions, microemulsions and foams, Springer, 2020.
4. D. L. Weaire and S. Hutzler, The physics of foams, Oxford University Press, 1999.
5. F. Aurenhammer and R. Klein, Voronoi Diagrams, FernUniversität in Hagen, Hagen, 1996.
6. F. Aurenhammer and H. Edelsbrunner, Pattern Recogn., 1984, 17, 251–257.
7. F. Aurenhammer, SIAM J. Comput., 1987, 16, 78–96.
8. E. D. Kadtsyn, V. A. Nichiporenko and N. N. Medvedev, J. Mol. Liq., 2022, 349, 118173.
9. J. S. Hunjan and B. C. Eu, J. Chem. Phys., 2010, 132, 134510.
10. D. Dudariev, V. Koverga, O. Kalugin, F. A. Miannay, K. Polok, T. Takamuku, P. Jedlovszky and A. Idrissi, J. Phys. Chem. B, 2023, 127, 2534–2545.
11. P. K. Kao, B. J. VanSaders, S. C. Glotzer and M. J. Solomon, Sci. Rep., 2021, 11, 11042.
12. D. P. Bourne, P. J. J. Kok, S. M. Roper and W. D. T. Spanjer, Philos. Mag., 2020, 100, 2677–2707.
13. G. Frenning, Comput. Part. Mech., 2023, 10, 427–443.
14. C. H. Zhang, S. W. Zhao, J. D. Zhao and X. W. Zhou, Powder Technol., 2021, 379, 251–264.
15. J. M. Sullivan, in Foams and Emulsions, ed. J. F. Sadoc and N. Rivier, Springer Netherlands, Dordrecht, 1999, ch. 23, pp. 379–402 DOI:10.1007/978-94-015-9157-7_23.
16. M. Z. Cao, T. W. Qiu, B. X. Deng, Y. K. An, Y. D. Xing and E. R. Zhao, Mater. Des., 2024, 238, 112729.
17. C. E. Flores, K. B. Sautter, P. Bucher, A. Cornejo, A. Franci, K. U. Bletzinger and R. Wüchner, Comput. Part. Mech., 2023, 10, 1181–1203.
18. E. A. Lazar, J. Y. Lu and C. H. Rycroft, Am. J. Phys., 2022, 90, 469–480.
19. A. Goede, R. Preissner and C. Frommel, J. Comput. Chem., 1997, 18, 1113–1123.
20. F. Anton, D. Mioc and M. Santos, J. Comput. Sci. Technol., 2013, 28, 255–266.
21. D. S. Kim, Y. Cho and K. Sugihara, Comput.-Aided Des., 2010, 42, 874–888.
22. Z. Hu, X. Li, A. Krishnamurthy, I. Hanniel and S. McMains, Comput.-Aided Des. Appl., 2017, 14, 572–581.
23. K. Olechnovic and C. Venclovas, J. Comput. Chem., 2014, 35, 672–681.
24. D. S. Kim, Y. Cho and D. Kim, Comput.-Aided Des., 2005, 37, 1412–1424.
25. R. Staritzbichler, N. Ristic, A. Goede, R. Preissner and P. W. Hildebrand, Nucleic Acids Res., 2021, 49, W685–W690.
26. F. Persson, P. Soderhjelm and B. Halle, J. Chem. Phys., 2018, 148, 215101.
27. K. Olechnovic, R. Banciul, J. Dapkunas and C. Venclovas, Proteins:Struct., Funct., Bioinf., 2025 DOI:10.1002/prot.26792.
28. J. Dapkunas, A. Timinskas, K. Olechnovic, M. Tomkuviene and C. Venclovas, Nucleic Acids Res., 2024, 52, W264–W271.
29. Z. W. Nie, Y. Y. Lin and Q. B. Tong, Int. J. Heat Mass Transfer, 2017, 113, 819–839.
30. T. Wejrzanowski, J. Skibinski, J. Szumbarski and K. J. Kurzydlowski, Comput. Mater. Sci., 2013, 67, 216–221.
31. Z. W. Nie, Y. Y. Lin and Q. B. Tong, Comput. Mater. Sci., 2017, 131, 160–169.
32. O. Busaryev, T. K. Dey, H. M. Wang and Z. Ren, ACM Trans. Graph., 2012, 31, 1–8.
33. Z. G. Fan, Y. G. Wu, X. H. Zhao and Y. Z. Lu, Comput. Mater. Sci., 2004, 29, 301–308.
34. A. Liebscher, Philos. Mag., 2015, 95, 2777–2792.
35. C. Lautensack, J. Appl. Stat., 2008, 35, 985–995.
36. F. Aurenhammer and H. Imai, Geometriae Dedicata, 1988, 27, 65–75.
37. W. M. Visscher and M. Bolsterli, Nature, 1972, 239, 504–507.
38. P. Richard, L. Oger, J. P. Troadec and A. Gervois, Eur. Phys. J. E:Soft Matter Biol. Phys., 2001, 6, 295–303.
39. P. Richard, L. Oger, J. P. Troadec and A. Gervois, Phys. A, 1998, 259, 205–221.
40. S. M. K. Rassouly, Powder Technol., 1999, 103, 145–150.
41. R. B. Aïm and P. L. Goff, Powder Technol., 1968, 1, 281–290.
42. L. Y. Yi, K. J. Dong, R. P. Zou and A. B. Yu, Powder Technol., 2012, 224, 129–137.
43. R. P. Zou, X. Bian, D. Pinson, R. Y. Yang, A. B. Yu and P. Zulli, Part. Part. Syst. Charact., 2003, 20, 335–341.
44. A. M. Kraynik, D. A. Reinelt and F. van Swol, Colloids Surf., A, 2005, 263, 11–17.
45. F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot and E. Duchesnay, J. Mach. Learn. Res., 2011, 12, 2825–2830.
46. K. E. Hoff, J. Keyser, M. Lin, D. Manocha and T. Culver, presented in part at the Proceedings of the 26th annual conference on Computer graphics and interactive techniques, 1999.
47. G. Rong and T.-S. Tan, presented in part at the Proceedings of the 2006 symposium on Interactive 3D graphics and games – SI3D’06, Redwood City, California, 2006.
48. Z. Yuan, G. Rong, X. Guo and W. Wang, Presented in part at the Eighth International Symposium on Voronoi Diagrams in Science and Engineering, 2011 DOI:10.1109/ISVD.2011.18.
49. N. Lindow, D. Baum and H. C. Hege, IEEE Trans. Visualization Comput. Graphics, 2011, 17, 2025–2034.
50. I. Z. Emiris and G. M. Tzoumas, presented in part at the Proceedings of the 2007 ACM symposium on Solid and physical modeling, Beijing, China, 2007.
51. G. Mollon and J. D. Zhao, Granular Matter, 2012, 14, 621–638.
52. S. Dorn, N. Wolpert and E. Schömer, 2020 Ieee International Conference on Robotics and Automation (Icra), 2020, pp. 137–143 DOI:10.1109/icra40945.2020.9196775.
53. F. M. Schaller, S. C. Kapfer, M. E. Evans, M. J. F. Hoffmann, T. Aste, M. Saadatfar, K. Mecke, G. W. Delaney and G. E. Schröder-Turk, Philos. Mag., 2013, 93, 3993–4017.
54. N. N. Medvedev, V. P. Voloshin, V. A. Luchnikov and M. L. Gavrilova, J. Comput. Chem., 2006, 27, 1676–1692.
55. S. A. Magrabi, B. Z. Dlugogorski and G. J. Jameson, Chem. Eng. Sci., 1999, 54, 4007–4022.
56. C. F. Revelo and H. A. Colorado, Process. Appl. Ceram., 2019, 13, 287–299.
57. R. Pielaszek, J. Alloys Compd., 2004, 382, 128–132.
58. I. Jongenburger, M. W. Reij, E. P. Boer, M. H. Zwietering and L. G. Gorris, Int. J. Food Microbiol., 2012, 157, 35–44.
59. M. F. Vaz and M. Fortes, Scr. Metall., 1988, 22, 35–40.
60. L. Gibson and M. Ashby, Cellular Solids: Structure and Properties, 1997, vol. 2, pp. 175–234.
61. R. S. Farr and R. D. Groot, J. Chem. Phys., 2009, 131, 244104.
62. A. J. DeVries, Rubber Chem. Technol., 1958, 31, 1142–1205.
63. B. Gal-Or and H. Hoelscher, AIChE J., 1966, 12, 499–508.
64. R. Lemlich, Chem. Eng. Commun., 1982, 16, 153–157.
65. K. Bagi, Granular Matter, 2005, 7, 31–43.
66. Y. F. Cheng, S. J. Guo and H. Y. Lai, Powder Technol., 2000, 107, 123–130.
67. R. R. Valera, I. P. Morales, S. Vanmaercke, C. R. Morfa, L. A. Cortés and H. D. G. Casañas, Comput. Part. Mech., 2015, 2, 161–172.
68. G. T. Nolan and P. E. Kavanagh, Powder Technol., 1992, 72, 149–155.
69. E. Lozano, D. Roehl, W. Celes and M. Gattass, Comput. Math. Appl., 2016, 71, 1586–1601.
70. J. E. Taylor, Ann. Math., 1976, 103, 489–539.
71. F. Deindoerfer and A. Humphrey, Ind. Eng. Chem., 2002, 53, 755–759.
72. F. Eslami and J. A. Elliott, J. Phys. Chem. B, 2014, 118, 14675–14686.
73. M. Manak, M. Zemek, J. Szkandera, I. Kolingerova, E. Papaleo and M. Lambrughi, J. Mol. Graphics Modell., 2017, 74, 225–233.
74. A. M. Kraynik, D. A. Reinelt and F. van Swol, Phys. Rev. Lett., 2004, 93, 208301.
75. K. M. Andersson and S. Hovmoller, Acta Crystallogr., Sect. D:Biol. Crystallogr., 2000, 56, 789–790.

---

Footnote

† Electronic supplementary information (ESI) available: Supplemental methods on FOAM_GEN and Fig. S1–S6. See DOI: https://doi.org/10.1039/d5cp00763a

---