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.