admin – Miloslav Torda, PhD

March 31, 2025April 2, 2025

Ionic Molecular Crystallisation Continued: Tetrahedra, Symmetry, and Salt-Like Molecular Architectures

In my previous post on Ionic Molecular Crystallization, I sketched a geometric analysis of a two-dimensional ionic molecular framework, guided by charged particle interactions, symmetry, and regular polygon packings. Since our real interest lies in three-dimensional frameworks, I wondered how this idea would translate into three dimensions. While still theoretical, the model feels promising—and surprisingly beautiful.

To help me visualise the constructions, I discovered a wonderful CAD modelling software, Shapr3D, which is both intuitive to use and very powerful. I had zero experience with CAD before, but within a few hours, I was able to build all kinds of three dimensional frameworks.

Let’s start with the simplest regular polyhedron—the regular tetrahedron—and build a periodic framework by connecting the vertices of tetrahedra, as shown in the picture below.

Unfortunately, regular tetrahedra are not space-filling. In fact, the densest known packing of tetrahedra is currently a dimer packing with $P\bar{1}$ space group symmetry and a density of $\frac{4000}{4671}$ (Chen, E.R., Engel, M., & Glotzer, S.C. Dense Crystalline Dimer Packings of Regular Tetrahedra. Discrete Comput. Geom. 44, 253–280 (2010). https://doi.org/10.1007/s00454-010-9273-0). However, I suspect that this packing might turn out to be a $P2_{1}/c$ or $C2/c$ tetrahedra monomer packing.

The question is: how many regular tetrahedra can share one vertex subject to crystallographic restrictions? The answer is eight, as shown in the image below. This is a simple consequence of the sphere kissing number being 12. It forms a uniform star polyhedron called an octahemioctahedron.

The next question is: what might a local cluster with low net electrostatic energy look like? Here is one such arrangement of blue ( $n_{+} = 8$ ) and green ( $n_{\text{-}}=6$ ) charged particles that satisfies this criterion under the Coulomb potential.

Additionally, I used the GEOMAG construction toy to build a physical model of this configuration. The green bars represent the repulsive forces in the green anion octahedral configuration. The white bars represent the attractive forces between green anions and blue cations; however, this is not entirely correct. These should not be represented as rigid rods but as tendons that can change length.

Geometrically, the GEOMAG structure represents the 1-skeleton of a stellated octahedron. The projection of the two tetrahedra, which together form the stellated octahedron, becomes visible in the shadow cast by the setting sun.

We set the configuration such that the local cluster’s Coulomb energy—introduced in the Ionic Molecular Crystallisation—is zero.

$\begin{equation*} E^{C}(C_0) = \sum_{\substack{\{i,j\} \in C_0 \\ i \neq j}} \phi(i,j) = 0 \end{equation*}$

In terms of an optimisation problem, the configuration is a solution to the following:

$\begin{equation*} C_0 = \text{argmin}_{\iota:C_0 \rightarrow \mathbb{R}^{2k}} \left[ E^{C}(C_0^{\iota}) \right]^2 \end{equation*}$

As in the two-dimensional case, the positively charged particles (blue) are located in the tetrahedral face-centred cubic (FCC) interstitial voids.

The local configuration of positively charged tetrahedral molecules and single-atom negatively charged molecules looks like this:

Once we assemble these octahemioctahedra, the resulting framework exhibits $Fm\bar{3}m$ symmetry (space group 225)—the same as fluorite ( $\text{CaF}_{2}$ ). The similarity is not merely mathematical; chemically, the tetrahedral framework model could approximate a real organic fluorite structure, provided we can identify the appropriate molecular analogue.

One can view this as the tetrahedral-octahedral honeycomb, where the octahedral building blocks have been removed. I positioned the camera to highlight the six-fold and four-fold roto-inversions of the framework.

Following the two-dimensional triangular tiling case, I constructed the following framework by removing tetrahedra from the $Fm\bar{3}m$ framework while preserving crystallographic symmetry:

The framework’s structure becomes clearer when in motion. See the animation below:

This is the three-dimensional equivalent of the lowest-density triangular configuration (Kagome lattice) from the Ionic Molecular Crystallisation post—the quarter cubic honeycomb. Its space group is $Fd\bar{3}m$ (space group 227).

As in the case of triangular tiling, there is a group-subgroup relationship between the tetrahedral-octahedral and the quarter cubic honeycombs via the space group $Pn\bar{3}m$ (space group 224). $Pn\bar{3}m$ is the maximal index $4$ $k$ -subgroup of $Fm\bar{3}m$ , and $Fd\bar{3}m$ is the maximal index $2$ $k$ -subgroup of $Pn\bar{3}m$ . This should not come as a surprise given how the framework was constructed.

If the observations from the two-dimensional framework analysis translate to three dimensions, this framework should represent a local Coulomb energy minimum, and a global minimum if the symmetries of the molecular system are constrained to the $Fd\bar{3}m$ space group symmetry isomorphism class.

Similar to triangular frameworks, one can associate a sphere packing with the $Fm\bar{3}m$ tetrahedral framework, and further with a regular $4$ -polytope, using these to enumerate local energy minima.

This suggests a potential geometric design principle: Higher symmetry often corresponds to greater mechanical stability, whereas its sub-frameworks may offer larger pores or other advantageous properties at the expense of rigidity.

In other words, if our Coulomb cluster-tetrahedral structure minimises energy under $Fm\bar{3}m$ symmetry, then the lower-symmetry frameworks derived from it (e.g., $Fd\bar{3}m$ ) may also represent local minima—stable enough to exist, but metastable within the full energy landscape.

Moreover, we can define a generative method for exploring neighbourhoods of local energy optima by sampling an entire family of ionic molecular frameworks through the simple identification of symmetry-breaking subgroups of a highly symmetrical parent structure. That is, we begin with the framework exhibiting the highest symmetry and then systematically investigate its symmetry-reduced subframeworks.

Interestingly, some coarse-grained molecular simulations already approximate molecules as rigid polyhedra. Our goal is to create stable framework models with large cavities, which may serve as blueprints for the computational design of organic frameworks, applicable to areas such as water and $\text{CO}_{2}$ capture.

In summary, referring back to previous posts (Chiral Interaction Ground States and Ionic Molecular Crystallisation), these configurations are all related through the face-centred cubic lattice packing of spheres, with a density of $\frac{\pi}{\sqrt{18}}$ ,

and the associated absolutely symmetric quadratic form:

$\begin{equation*} X^2 + Y^2 + Z^2 + XZ + YZ \end{equation*}$

courtesy of K. L. Fields (K. L. Fields, 1979. Stable, fragile, and absolutely symmetric quadratic forms. Mathematika, 26(1), 76–79).

This brings us full circle to the symmetries of FCC sphere packing explored in these blog posts:
$\text{Densest } P1 , P\bar{1} \text{, and} P2_{1}/c \text{ packings of spheres}$

$\text{Sphere packings continued - space groups} P2_{1} \text{ and } C2/c$

Based on these, A. I. Kitaigorodsky (link to German Wikipedia page as no English version exists) formulated his Close Packing Principle:

“The mutual arrangement of the molecules in a crystal is always such that the ‘projections’ of one molecule fit into the ‘hollows’ of adjacent molecules.”

— From Molecular Crystals and Molecules by A. I. Kitaigorodsky.

So, what’s next? I’m now looking into how kinematics of the cuboctahedron under space group symmetry constraints can help identify which polyhedral frameworks are mechanically stable—and which ones are just aesthetically pleasing. After my visit to ICERM and a conversation with Robert Connelly, I’m convinced that the mathematics of tensegrities is the key.

March 17, 2025March 29, 2025

What Do Geodesic Polyhedra and X-ray Crystallography Have in Common?

During her visit to Liverpool, Zuzka—my fiancée Lina’s mom—gifted us these glass geodesic polyhedra. When sunlight shines through them, it scatters across our living room walls, creating miniature rainbows.

While preparing to travel to ICERM’s workshop, Matroids, Rigidity, and Algebraic Statistics, I noticed that the geodesic polyhedra formed a dilated triangular lattice pattern on the wall.

This is, in fact, the basic idea behind X-ray diffraction, discovered by Max von Laue in 1912. When electromagnetic waves-such as X-rays-pass through a crystal, the symmetrical arrangement of atoms causes the waves to diffract. When these waves hit an impenetrable surface, we can photograph the resulting wave interference pattern, composed of both constructive and destructive interactions.

One key challenge in X-ray crystallography is reconstructing a crystal’s atomic structure from its diffraction pattern—an example of an inverse problem. For instance, the mini-rainbow triangular lattice in the photo arises from the six-fold rotational symmetry of our two glass tetrahedral geodesic polyhedra.

One can perceive the situation at hand as a stereographic projection of a spherical model (as in Magnus Wenninger‘s book, Spherical Models) onto the two-dimensional plane, thereby relating the infinite triangular lattice to the “sphere of reflection”-i.e., Ewald’s sphere.

The idea of hanging transparent geodesic shapes in front of the window—creating rainbows across the room—came from Lina’s dad and my teacher, Oleg Šuk. He experimented with different materials, including an acrylic geodesic heart (with “acrylic” referring here to PMMA photonic crystal). I can’t help but wonder if the Šuk family is secretly guiding my work. See my previous posts, Chiral Interaction Energy Ground States and Yayoi Kusama’s Chandelier of Grief: Symmetry in Art and Science, for additional context.

March 5, 2025March 31, 2025

Ionic Molecular Crystallization

A few months ago, I encountered an interesting problem related to molecular crystallization. Consider a two-dimensional system with two types of molecules: the first type consists of rigid triangles with positive charges localized at the vertices (shown as blue triangles in the image below), and the second type is a single-atom molecule with a negative charge (shown as a green dot).

Our key questions are:

Does this molecular system crystallize in the thermodynamic limit?
If so, what are the characteristics of its ground states?

Let’s assume the system’s intermolecular interactions follow a normalized Coulomb/electrostatic potential:

$\begin{equation*} \phi_{(i,j)} = \frac{q_i q_j}{r_{(i,j)}}, i \in I, j \in J \end{equation*}$

where $I$ and $J$ denote molecules of either type (green negative or blue positive triangular), $r_{(i,j)}$ is the Euclidean distance between charges $q_i$ and $q_j$ . For green particles, $q_i=-1$ , while for vertices of blue triangles, $q_i=1$ . This creates attractive forces between opposite charges and repulsive forces between like charges, proportional to their relative distance.

The intermolecular energy is defined as the sum of all pairwise interactions $\phi_{(i,j)}$ between molecules $I$ and $J$ :

$\begin{equation*} E(I,J) = \sum_{i \in I,j \in J} \phi_{(i,j)} \end{equation*}$

This quantity generally lacks an infimum, so we need additional constraints.

We aim to find a sequence of molecular configurations $X_1 \subset X_2 \subset \ldots \subset X_N \subset \ldots$ where the following limit converges to a constant:

$\begin{equation*} \lim_{ \ |X_N|\rightarrow \infty} \frac{1}{| X_N|} \min_{\iota: X_N \rightarrow \mathbb{R}^{(2N + n_{+})} } \sum_{\substack{ \{I,J\}| I \neq J \\ I,J \in X_N }} E(I,J) = \epsilon \end{equation*}$

Here, $X_N$ represents a collection of $N$ molecules of both types ( $N=n_{-} + n_{+}$ ).

Let’s focus on local clusters, as shown below:

For any $k$ where $(n_{-} + n_{+}) = 2k$ , an alternating configuration of charges in a local cluster can found for any $\epsilon$ . Moreover, due to the rigidity of triangular molecules and the form of the potential $\phi$ , the maximal cluster contains $n_{-}=n_{+} = 6$ charges.

In our example cluster, $C$ , we set the energy to:

$\begin{equation*} E^{C}(C_0)= \sum_{\substack{\{i,j\} \in C_0 \\ i\neq j}}\phi(i,j) = -1 \end{equation}$

for demonstration purposes.

Alternatively, the cluster configuration can be expressed in the form of the optimization problem:

$\begin{equation*} C_0 = \text{argmin}_{\iota:C_0 \rightarrow \mathbb{R}^{2k}} \left[ E^{C}(C_0^{\iota}) +1 \right]^2 \end{equation}$

In the $C_0$ configuration, the blue positive charges occupy the interstitial voids of the hexagonal close-packed arrangement of negative green charges:

Following these local optimizations, we define a quantity we call the energy per cluster as:

$\begin{equation*} E^C_{\mathcal{A}}(\xi)=E^{C}(C_0)+\sum_{\substack{\{i,j\} \\ i\in C_0 \\ j \in C_\alpha | \alpha \in \mathcal{A}}}\phi(i,j) \end{equation}$

where $\mathcal{A}$ is an index set . Minimizing $E^C_{\mathcal{A}}(\xi)$ for such a cluster collection becomes a well-defined problem:

$\begin{equation*} \min_{\mathcal{\xi}:\mathcal{A} \rightarrow \mathbb{R}^{2|\mathcal{A}|}} E^C_{\mathcal{A}}(\xi) \end{equation}$

For cardinality $|\mathcal{A}| = 6$ , the solution is illustrated below.

Note that the interaction potential $\phi(i,j)$ respects molecular memberships—for $i \in I$ and $j \in J$ , we require $I \neq J$ .

We can expand $E_C^\Alpha$ to $|\mathccal{A}|=18$ :

or any $|\mathcal{A}|=6(2n-1)$ .

What we did is that we transformed the problem of finding molecular system minimizers into finding minimizers of a single-particle system with potential $\frac{1}{6}E^C_{\mathcal{A}}$ for $|\mathcal{A}|=6$ , and we already know that the global minimizer is a dilated, translated, or rotated version of the triangular lattice.

Geometrically, this transformation effectively contracts local clusters into single points, resulting in a regular triangular tiling:

Having resolved the crystallization question, let’s examine the ground states’ characteristics, starting with a natural question: Are there other possible crystal phases? (From the materials science perspective, this question is important because of Crystal polymorphism)

On one side we can view our ground state configuration as a triangular tiling because of the rigidity of the blue molecules. On the other side, because of the way we have reframed the intermolecular interactions in terms of pairwise cluster interactions we may view the ground state configuration as the triangular lattice. The connecting point of these two views are the Archimedean circle packings and their associated semi-regular tessellations.

The densest packing of regular triangles is the $p3m1$ configuration with a density of $1$ . It is, in fact, a tiling of triangles as shown below.

The maximal non-isomorphic subgroup of $p3m1$ is the $p31m$ plane group. The densest packing of regular triangles when the packing configurations are restricted to the $p31m$ isomorphism class has a density of $\frac{3}{4}$ and is shown in the image below.

This is, in fact, a consequence of removing one of the mirror symmetries from the $p3m1$ plane group.

Another way to get here is to look at $p3m1$ , the Archimedean framework

associated with the $p3m1$ disc packing with packing density of $\frac{\sqrt{3}}{9}\pi$

, and the Archimedean framework associated with the hexagonal close packed configuration of disc with density $\frac{\sqrt{3}}{6}\pi$

and its associated Archimedean framework

By interlacing the hexagonal and triangular frameworks

we construct new framework. The common (blue) intersection point of the hexagona and triangular frameworks in fact define a triangular sub-tiling

and the p $31m$ is constructed by removing triangle tiles containing a vertex at its center of mass

Alternatively we can create a complementary triangle packing with density of $\frac{1}{4}$ by removing tiles not containing an interior vertex

The resulting structure is sometimes referred to as the Kagome Lattice, although it is not a lattice per se.

Starting from the $cm$ triangle packing which is in fact a tiling

and decomposing it into two complementary packings with $pg$ symmetry with density of $\frac{1}{2}$

This should come as no surprise since $cm$ is the a maximal t-subgroup of $p3m1$ of index 3.

To demonstrate this point, we start again from the $p3m1$ disc packing framework

and interlace it with the triangular framework associated with $p1$ packing of discs

in such a way that disc centers are are subsets of the hexagon vertices

and we remove one $pg$ symmetry from the tiling, the one not having a hexagon vertex at its center.

In fact, there are two possibilities to achieve the same result by using the triangle tiling above and rotating it by $\pi$ with respect to zero (center of the middle pentagon).

and do exactly the same thing as previously to get the second $pg$ triangle tiling

Combining these two triangle packings,

the $cm$ symmetry of $p3m1$ should be clearly visible.

This tiling is a isomorphic to the one of the Archimedean (semi-regular vertex transitive) tilings associated with the densest $cm$ packing of discs with density of $(2-\sqrt{3})\pi$

with it’s $3^3.4^2$ semi-regular tiling

These constructions are nothing new. For a comprehensive treatment of regular and semi-regular tilings, I recommend Regular Polytopes by H. S. M. Coxeter, first published in 1947. (Coxeter was the geometer who inspired M. C. Escher, which resulted in the Circle Limit I–IV series.)

Returning to our original problem of characterizing and enumerating possible Coulomb potential ground states of our molecular system, all these triangle packings represent viable local energy minima. This conclusion follows from our transformation of the intermolecular interactions into a nearly spherically symmetric “energy per cluster” potential, which preserves all symmetries of the hexagonal close-packing configuration of discs: p2, p2gg, pg, p3, and p1, as well as the rigidity of the triangular molecule.

We have classified the packings by associating each configuration with Archimedean disc packings and tessellations. The global energy minimizer corresponds to the triangular tiling and lattice, as demonstrated by Theil in “A Proof of Crystallization in Two Dimensions” (Commun. Math. Phys. 262, 209–236, 2006). https://doi.org/10.1007/s00220-005-1458-7. Since all our configurations are subsets of the triangular tessellation relative to some tiling symmetry subgroup, they are energy minimizers within their respective symmetry groups. (See Bétermin, L. “Effect of Periodic Arrays of Defects on Lattice Energy Minimizers,” Ann. Henri Poincaré 22, 2995–3023, 2021. https://doi.org/10.1007/s00023-021-01045-0.)

A property of Archimedean tessellations is that they are vertex-transitive, meaning the same number, $q$ , of regular $p$ -gons meet at each vertex. In triangular tessellations, six regular triangles ( $3$ -gons) meet at each vertex. For regular tessellations like the triangular one, this follows a simple equation:

$\begin{equation*} q\left( 1- \frac{2}{p} \right)\pi = 2\pi \end{equation}$

This is only for one type of polygon, but we have more types. However, we can generalize this formula to n-types of polygons:

$\begin{equation*} q_1\left( 1- \frac{2}{p_1} \right)\pi + q_2\left( 1- \frac{2}{p_2} \right)\pi + \ldots + q_n\left( 1- \frac{2}{p_n} \right)\pi = 2\pi \end{equation}$

For example, the tessellation associated with the $cm$ packing of discs shown above has every vertex surrounded by $q_1 = 3$ regular $p_1 = 3$ -gons (triangles) and $q_2 = 2$ regular $p_2 = 4$ -gons (squares).

Ultimately, the equation above is simply a disguised form of the Euler characteristic for plane-connected graphs induced by the vertex figure of a semi-regular tessellation.

What does this mean in the context of our molecular system?

The vertex figure of each ground state is determined by the leading contributions to its total energy.
In training Graph neural networks for Machine-learned interatomic potentials, learning these interactions alone suffices for ground state prediction in such molecular systems, significantly improving computational efficiency.
Combined with the Crystallographic restriction theorem, this formula allows us to explore possible ground state configurations based on molecular shape. In a periodic configuration, for every $p_i$ , there exists at least one $k \in \{2,3,4,6\}$ such that the expression $p_i \ \pmod{k} = 0$ holds. Thus, for any such $p_i$ -gon shaped molecule, the vertex figure equation has only a finite number of integer solutions. Unfortunately, in cases like the pentagon, we might need to delve into the realm of aperiodic tilings.

December 24, 2024February 12, 2025

Chiral Interaction Energy Ground States

A few days ago, my beloved fiancé surprised me with a GEOMAG construction toy for my birthday. GEOMAG consists of small rods with embedded magnets and metallic spheres that can be assembled into various structures. It’s truly amazing—I had no idea something like this existed!

This toy can effectively model crystallography problems. For example, below is a geometric proof of the Crystallographic restriction theorem

It can also demonstrate structural rigidity problems, as shown in this framework

which becomes flexible after removing just one link.

The stable framework is actually the snub trihexagonal tiling – a regular tessellation associated with the p6 packing of discs.

This framework is created by drawing links between the centers of touching disks. The tessellation on the right is called a floret pentagonal tiling. We can create this tiling from the p6 packing of hexagons by connecting each hexagon’s edge vertices to their nearest lattice points.

From the crystallographic perspective, the vertices of the dual tessellation lie in the interstitial sites of the regular tessellation. Thus, the circles of the p6 packing are incircles of the pentagons in the floret pentagonal tiling and the hexagons in the hexagonal p6 packing.

This serves as a physical model of the Lennard-Jones system I described in my blog entry Lennard-Jones hexagonal molecular system. Through this model, I discovered that there are two chiral ground states connected by a continuous path in the configuration space.

Take a look at the animation below to see how it all works!

This is closely related to the vibrational part of lattice energy, degrees of freedom and structural stability of molecular crystals. Let’s have closer looks at this. Below is a visualization illustrating the transition between the chiral ground states of a 3-molecule cluster. Each molecule in the cluster contains six atoms. The intermolecular energy is given by the Lennard-Jones potential:

$\begin{equation*} E_{LJ}=\sum_{\substack{\{I,J\} \\ I \cap J =0}} \sum_{\substack{(i,j ) \\ i \in I, j \in J}} \left[\left(\frac{1}{r_{(i,j)}}\right)^{12} - \left(\frac{1}{r_{(i,j)}}\right)^{6}\right] \end{equation*}$

where $r_{( i,j ) } = \| x_i - x_j \|$ is the Euclidean distance between the $i$ -th atom of molecule $I$ and the $j$ -th atom of molecule $J$ , as shown in the image below.

The chiral ground states are in clear energy potential wells and represent two different global solutions of the energy cluster $E_{LJ}$ minimization problem.

From the GEOMAG model animation, one can imagine this cluster being held together only by the bonds between nearest neighbors. By allowing three of the total bonds to stretch, one introduces one degree of freedom into the otherwise rigid cluster, making it possible for the configuration to move to its neighboring potential well.

This however requires a short range interaction potential, since the pairwise interaction is computed among all atoms of different molecules not only the first neighbors. The Lennard-Jones potential is one such example where the leading contribution to the overall energy is the three atom cluster in the middle and the transition from one ground state to the other is animated as revolving around the centre of mass of this cluster.

It is also by this Lennard-Jones potential property that the ground state of this two-dimensional molecular system coincides with the densest $\mathbf{p6}$ packing of discs. This means that there are also two chiral densest $p6$ packings of all $n$ -gons with six-fold rotational symmetry. This observation connects directly to my earlier work on plane group packings during my PhD, as it shows an aspect these packings I overlooked. The results of this study are published in our Densest plane group packings of regular polygons manuscript.

Moreover, this explains the crystal defects in the Lennard-Jones hexagonal molecular system. Since mirror symmetry is not permitted in this two-dimensional system, we observe incompatible local ground state patches.

I was wondering if a real material with this crystal structure exists. In fact, it does, and it was published only a few years ago. See the image and reference below.

Rusek, M., Kwaśna, K., Budzianowski, A., & Katrusiak, A. (2019). Fluorine··· fluorine interactions in a high-pressure layered phase of perfluorobenzene. The Journal of Physical Chemistry C, 124(1), 99-106.

The above is a space-filling visualization of this high-pressure phase of perfluorobenzene, which is in excellent agreement with our highly simplified model. The sphere radii are set to the van der Waals radii. The left image shows one layer, while the right image displays the same layer with the carbons removed. The space group of this crystal is C2/c. C2/c symmetry consists of a two-fold rotational symmetry (here, within one layer of perfluorobenzene) and a perpendicular mirror symmetry. This effectively means that the layers alternate as two chiral two-dimensional ground states, interchanging chirality layer by layer.

Why is this important beyond just creating appealing pictures and animations? One reason is that it can significantly reduce computational time complexity by eliminating unnecessary calculations in Crystal Structure Prediction (CSP). For example, if we know beforehand that our target compound satisfies a few assumptions, the degrees of freedom to explore can be reduced to a maximum of 12, sometimes even less, as demonstrated in our Close-Packing approach to CSP (Alternatively, we can think of the hard sphere model as a pairwise interaction potential. However, in this blog, we are considering soft core interactions).

And we know time is money, as illustrated by the recent developments around the DeepSeek model and its ensemble learning approach to large language model training. The concept of ensemble here mirrors the Thermodynamic ensemble that Gibbs introduced in 1902.

This also shows how chemistry can inspire the development of novel machine learning algorithms beyond Machine-learned interatomic potentials.

November 8, 2024January 30, 2025

LRC Symposium Supplementary Materials

Visualisation of the MCYPDE packing maximisation run for the space group P2₁ (Left) Maximum packing density; (middle) average packing density of N-best packings at i-th iteration; (right) distribution of the packing generation at the 100-th iteration.

Download handout

September 5, 2024January 29, 2025

Ellsworth Kelly’s Yellow Curve: Axiomatic Art

Below is a photo I took at the exhibition of works by Ellsworth Kelly at the Louis Vuitton Foundation, which I visited on my way back from the CaLISTA Workshop on Geometry-Informed Machine Learning at Mines Paris.

He truly distilled art to its essentials—I’m sure Euclid would have approved.

August 19, 2024August 20, 2024

Sphere packings continued – space groups $P2_1$ and $C2/c$

In this post, we’re diving back into the world of sphere packings, focusing on the space groups $P2_1$ and $C2/c$ . Feel free to check out an earlier discussion on Densest P1, P-1 and P21/c packings of spheres. Both of these groups reach the packing density of $\frac{\pi}{\sqrt{18}}$ .

The densest $P2_1$ packing almost identical to the $P\overline{1}$ group. Even though they belong to different crystal systems—triclinic for $P\overline{1}$ and monoclinic for $P2_1$ —they share the same unit cell parameters: $a = c = 2$ , $b = 2\sqrt{2}$ , and $\beta = \frac{\pi}{2}$ . Here’s a snapshot of what this packing looks like:

Next up is the $C2/c$ space group, which also falls under the monoclinic category. This group is a bit more complex, with eight symmetry operations. Its densest packing has unit cell parameters of $a = b = 4$ , $c = 2\sqrt{3}$ , and $\beta = \pi - \text{arcsin} \left( \frac{\sqrt{2}}{\sqrt{3}} \right)$ . Check out the illustration below of the densest packing configuration:

So far, we’ve seen that the optimal packing density for the space groups $P1$ , $P\overline{1}$ , $P21/c$ , $P2_1$ , and $C2/c$ matches the general optimal sphere packing density. This isn’t too surprising since these groups are all related to the Hexagonal Closed Packed structure, via group – subgroup relations. Interestingly, these space groups are among the top ten of the most frequently occurring in the Cambridge Structural Database, making up 74% of its entries.

July 3, 2024August 20, 2024

Yayoi Kusama’s Chandelier of Grief: Symmetry in Art and Science

This photo captures my fiancé Lina and myself at Yayoi Kusama’s “Infinity Mirror Rooms” exhibit at Tate Modern in London, taken in March 2024. The installation, called “Chandelier of Grief,” uses mirrors and light to explore ideas about symmetry.

The centrepiece of the artwork is a chandelier surrounded by six mirrors arranged in a regular hexagonal prism. This creates an interesting visual effect—an infinite lattice of chandeliers extending in all directions, with each reflected chandelier appearing at the centre of symmetry within this lattice.

What we’re experiencing is a visual representation of a hexagonal lattice, one of the five two-dimensional Bravais lattices. Mathematically, it represents a discrete group of isometries in two-dimensional Euclidean space. The point group of this hexagonal lattice is isomorphic to the dihedral group $D_6$ , highlighting its sixfold rotational symmetry and mirror reflection properties.

Physically, the hexagonal lattice emerges as the ground state configuration for systems of particles interacting via certain potentials, such as the Lennard-Jones potential, in the thermodynamic limit. Moreover, the hexagonal lattice is realised in materials like graphene, where carbon atoms arrange themselves in this highly symmetric pattern, leading to its unique properties.

Kusama’s installation is not only aesthetically pleasing but also serves as a visual representation of abstract concepts in mathematics and physics. The infinite reflections illustrate how artistic expression can intersect with scientific principles, making complex ideas accessible and engaging to a wider audience.

October 29, 2023January 30, 2025

Exercise in Monte-Carlo packing density computation.

I needed a way to calculate the packing densities of structures modeled using the van der Waals sphere model, generated through CSP and DFT computations. Specifically, for periodic configurations, this comes down to computing:: $\rho = \frac{\text{vol}(\mathbf{O})}{\text{vol}(\mathbf{U})}$ . Here, $\text{vol}(\mathbf{U})$ represents the volume of the unit cell $\mathbf{U}$ , and $\text{vol}(\mathbf{O})$ is the volume of the subset of $\mathbf{U}$ occupied by the van der Waals spheres denoted as $\mathbf{O}$ .

Calculating the volume of $\mathbf{U}$ is is straightforward; it’s simply the determinant of the unit cell’s generators: $\text{vol}(\mathbf{U})}=\text{det}(\mathbf{U})$ . However, analytically computing the volume of $\mathbf{O}$ involves solving numerous intricate integrals of intersections of spheres with various radii, making it cumbersome.

Let’s try to do this the Monte-Carlo way, instead. The volume of $\mathbf{O}$ can be expressed as an integral over the unit cell $\mathbf{U}$ , i.e.

$\begin{equation*}\text{vol}(\mathbf{O})=\int_{\mathbf{O}} d\mathbf{V}\end{equation*}$

Here, $d\mathbf{V}$ is the natural volume form on $\mathbf{R}^3$ . Since this integral is zero everywhere except in the unit cell, it can be redefined as an integral of the indicator function over the unit cell:

$\begin{equation*}\text{vol}(\mathbf{O})=\int \bm{1}_{\mathbf{O}} d\mathbf{U}\end{equation*}$

By changing coordinates to an integral over the unit cube $\textbf{C}$ , we get:

$\begin{equation*}\text{vol}(\mathbf{O})=\int \bm{1}_{\mathbf{U}^{-1}\mathbf{O}} \text{det}(\mathbf{U}) d\mathbf{C}\end{equation*}$

Plugging this into the density equation yields:

$\begin{equation*}$\rho=\int \bm{1}_{\mathbf{U}^{-1}\mathbf{O}} d\mathbf{C}\end{equation*}$

Considering $\mathbf{C}$ as a uniform random vector on the unit cube, this integral represents the expected value of the indicator function with respect to the $3$ -variate uniform distribution. This, in turn, equals the probability of $X \sim \mathbf{C}$ lying in the collection of the van der Waals spheres $\mathbf{U}^{-1}\mathbf{O}$ , scaled by the inverse of the unit cell:

$\begin{equation*}$\rho=\textbf{E}\left[ \bm{1}_{\mathbf{U}^{-1}\mathbf{O}} \right]=\text{P}\left( X \in \mathbf{U}^{-1}\mathbf{O} \right)\end{equation*}$

This means that the crystal’s density can be estimated by sampling $X_1, \ldots X_N \sim \mathbf{C}$ , counting how many fall inside $\mathbf{U}^{-1}\mathbf{O}$ , and dividing by $N$ :

$\begin{equation*}\hat{\rho} = \frac{\# \left \{ X_i \in \mathbf{U}^{-1}\mathbf{O} \right \}}{N}\end{equation*}$

Now, say I want to determine the density of this $P21C$ anthracene structure

Equivalently, I can sample uniformly from the parallelepiped defined by the unit cell, e.g. $100 000$ realisations,

and count how many fall inside at least one of the van der Waals spheres

which is $73056$ . Therefore, the density estimate is $\hat{\rho}=0.73056$ .

But how accurate is this estimate? Or even better, how large does $N$ need to be so that $| \hat{\rho} - \rho | < c$ , where $c$ is some constant?

To answer this, consider each $X_i$ as a Bernoulli trial. Then $\mathbm{X}=\sum_{i=1}^N X_i$ is a binomially distributed random variable with the expected value $N\rho$ and variance $N\rho\left(1-\rho \right)$ . According to the central limit theorem, the random variable $\mathbm{Z}=\frac{\sqrt{N}\left( \hat{\rho} - \rho \right)}{\sqrt{\rho\left(1-\rho \right)}}$ converges in distribution to a normally distributed random variable with a mean of $0$ and a variance of $1$ . Therefore, I can request $\mathbm{Z}$ being in the interval $(-z, z)$ , and I want the size of this interval to be less than some constant with probability $p$ :

$\begin{equation*} P \left( | \hat{\rho} - \rho | < 2 z \sqrt{\frac{\hat{\rho}\left(1-\hat{\rho} \right)}{N}} \right)=p \end{equation*}$

where $z$ is the $\frac{1+p}{2}$ quantile of the standard normal distribution.

For the $P21C$ anthracene structure mentioned earlier, to achieve $| \hat{\rho} - \rho | < 0.001$ with probability $p=0.999$ , I need $N=9,000,000$ . The density estimate in this case is then:

$\begin{equation*} \hat{\rho} \approx 0.731501 \end{equation*}$

As a point of comparison, CCDC’s Mercury reports the packing coefficient of this structure as $0.733255$ .

Closing note: This method could also be used to quantify uncertainty about the true packing during ETRPA‘s optimization run in the tradition of Sequential analysis.

October 15, 2023January 30, 2025

Repairing Unfeasible Solutions in ETRPA.

The ETRPA is designed to optimize a multi-objective function by:

1. Maximizing packing density.
2. Minimizing overlap between shapes in a configuration.

This is achieved through the use of a penalty function, which generally works quite effectively. However, there are some downsides to this approach. Many generated samples are flagged as unfeasible, leading to inefficiency. Additionally, some potentially good configurations with minimal overlaps are discarded. A more effective approach would involve mapping unfeasible solutions to feasible ones, essentially repairing the overlaps.

I’ve designed a relatively straightforward method to accomplish this almost analytically. It boils down to finding roots of a second-order polynomial, which is essentially a high school algebra problem. Admittedly, you’ll need to solve quite a few of these second-order polynomials, but they are relatively easy to handle.

So, with this approach, you can now map a configuration with overlap/s

to one without overlaps

Now, a couple of new issues have cropped up. First, there are some seriously bad solutions that are being sampled. For example, consider this one:

which is then mapped to this configuration

This is not helpful at all. Second, the parameters of this configuration are way above the optimisation boundaries that were used to create all configurations. As a result, it falls outside the configuration space where my search is taking place.

To address this issue, after repairing the configuration, I conduct a check to determine if the repaired configuration falls within my optimization boundaries. Any configurations that fall outside of these bounds are labeled as unfeasible solutions.

It might seems this puts me back in the same situation of multi-objective optimization, but that’s not the case. In my scenario, the Euclidean gradient at time $t$ takes this form:

$\nabla_{\bm{\theta}} J(\bm{\theta}^t) \ = \ \bm{\mu}^t_{\textbf{F}_{1- \frac{1}{q^{*}}}}-\bm{\mu}^t$ .

Here, $\bm{\mu}^t$ represents the mean of all samples in the current batch, and I only consider the configurations without overlap or the repaired ones that lie within the optimization boundaries for $\bm{\mu}^t_{\textbf{F}_{1- \frac{1}{q^{*}}}}$ . This effectively means I’ve eliminated the need for the penalty function and am optimizing only based on the objective function given by the packing density.

This is a promising, but it does introduce another problem. When I repair some of the unfeasible configurations, I end up making slight alterations to the underlying probability distribution, shifting it slightly from $\bm{\theta}^t$

to $\bm{\theta}_{\text{rep}}^t$

The issue here is that if I were to compute the Natural gradient as I used to, I’d obtain $\widetilde{\nabla} J(\bm{\theta}_{\text{rep}}^t) =\mathcal{I}_{\bm{\theta}_{\text{rep} }^t}^{-1} \left({\bm{\mu}_{\text{rep}}^t}_{\textbf{F}_{1- \frac{1}{q^{*}}}}-\bm{\mu_\text{rep}}^t\right)$ in the tangent space $\mathcal{T}_{\bm{\theta}_\text{rep}^t}$ of the statistical manifold, where everything is unfolding, but what I actually need is $\widetilde{\nabla} J(\bm{\theta}^t)$ in the tangent space at the point $\bm{\theta}^t$ .

Fortunately, the statistical manifold of the exponential family is a flat Riemannian manifold. So, first, I move from $\bm{\theta}^t$ to $\bm{\theta}_{\text{rep} }^t$ using this vector $\mathcal{I}_{\bm{\theta}^t}^{-1} \left(\bm{\mu_\text{rep}}^t-\bm{\mu}^t\right)$ (where $\bm{\mu}^t$ corresponds to the dual Legendre parametrization of $\bm{\theta}^t$ , and $\bm{\mu}_\text{rep}^t$ to $\bm{\theta}_\text{rep}^t$ ), and then I compute the natural gradient at the point $\bm{\theta}_\text{rep}^t$ .

In addition, to increase the competition within my packing population, I introduce a sort of tournament by retaining the N-best configurations I’ve found so far and incorporating them with the repaired configurations from the current generation. This leads to yet another slightly modified distribution, $\bm{\theta}_\text{rep + N-best}^t$

and I have to perform the same Riemannian gymnastics once more. In the end, the approximation to the natural gradient at the point $\bm{\theta}^t$ is as follows

$\widetilde{\nabla} J(\bm{\theta}^t)$ $=$ $\mathcal{I}_{\bm{\theta}_{\text{rep + N-best} }^t}^{-1} \left({\bm{\mu}_{\text{rep + N-best}}^t}_{\textbf{F}_{1- \frac{1}{q^{*}}}}-\bm{\mu}_{\text{rep + N-best} }^t\right)$ $+$ $\mathcal{I}_{\bm{\theta}_{\text{rep} }^t}^{-1} \left(\bm{\mu}_{\text{rep + N-best}}^t-\bm{\mu_\text{rep}}^t\right)$ $+$ $\mathcal{I}_{\bm{\theta}^t}^{-1}\left(\bm{\mu_\text{rep}}^t-\bm{\mu}^t\right)$

So far, experiments show that these modifications to the ETRPA are far more effective than the original ETRPA. Here’s a visualization from one of the test optimization runs:

In this particular run, the best packing I achieved attains the density of $0.7845$ . Previously, reaching such a level of density required multiple rounds of the refinement process, each consisting of $8000$ iterations. Now, a single run with just $2000$ iterations is sufficient. It’s worth noting that even after only $200$ iterations, the algorithm already reaches areas with densities exceeding $0.7$ .

It’s important to mention that currently, I’ve only implemented the repair process for unfeasible configurations for the $P1$ space group. However, similar repairs can be done for configurations in any space group. Furthermore, the algorithm’s convergence speed can be further enhanced by fine-tuning its hyperparameters. Perhaps it’s time to make use of the BARKLA cluster for the hyperparameter tuning.