Cluster Expansion

The cluster expansion method is a method to construct a complete basis that can fully describe any property that depends on the particular arrangement of the species onto a lattice. A system composed of a lattice of $N$ sites and $K$ different species has $K^N$ possible arrangements.

Formalism

The formalism is mainly attributed to the work of Sanchez and collaborators.

Site-basis Functions

Given $K$ site-basis functions $\varphi(\sigma)$ that assigns a value to the input specie, one can map the species sitting on the $N$ lattice sites of any arrangement $\vec{\sigma}$ as:

\[\begin{split}\sigma_i \longrightarrow \vec{\varphi}(\sigma_i) = \begin{bmatrix} \varphi_1(\sigma_i) \\ \varphi_2(\sigma_i) \\ \vdots \\ \varphi_K(\sigma_i) \\ \end{bmatrix} \in \mathbb{R}^K \quad \forall i \in 1, 2, \ldots, N\end{split}\]

A common basis for site-basis functions is typically the so-called Chebyshev basis. In a system of 3 species, the basis reads:

\[\begin{split}B_{chebyshev} = \begin{bmatrix} 1 & 1 & 1 \\ -\sqrt{3/2} & 0 & \sqrt{3/2} \\ \sqrt{1/2} & -\sqrt{2} & \sqrt{1/2} \end{bmatrix}\end{split}\]

Each row is a site-basis function, while each column is a specie.

Cluster Functions

By taking tensor products of these evaluated site-basis functions across all $N$ lattice sites, one can construct a complete basis of cluster function $\Phi(\vec{\sigma})$ as:

\[\Phi_\vec{\alpha}(\vec{\sigma}) = \prod_{(i,\nu) \in \vec{\alpha}} \varphi_\nu(\sigma_i)\]

where $\vec{\alpha}$ is a size-$N$ list of tuples containing the lattice site index $i$ and the associated index $\nu$ of the site-basis function for this site.

If the basis contains a constant site-basis function, i.e., a site-basis function that returns 1 for all species (e.g., the first row of the Chebyshev basis), then lattice sites associated with this site-basis function do not contribute to the outcome value of the cluster function, regardless of the species sitting on these sites.

The cluster function that associates the constant site-basis function with all $N$ lattice sites always returns the value 1. It is called the empty cluster function.
The cluster functions that associate the constant site-basis function with all lattice sites except one returns the value of the non-constant site-basis function evaluated at the specie sitting on that lattice site. These functions are called point cluster functions.
The cluster functions that associate the constant site-basis function with all lattice sites except two returns the value of the product of the two non-constant site-basis functions evaluated at the species sitting on these two lattice sites. These functions are called pair cluster functions.
In general, the set $C=\{\vec{s}_i\}$ of lattice sites associated with non-constant site-basis functions form a so-called cluster. The number of sites within $C$ defines the body-order of the cluster. It allows for a hierarchical construction of cluster functions.

There are exactly $K^N$ possible different cluster functions. When the site-basis functions are orthogonal to each other, the basis formed by the cluster functions is orthogonal as well. Any property that is a function of the arrangement on the lattice can be expanded as a linear combination of cluster functions (i.e., a cluster expansion) as:

\[f(\vec{\sigma}) = \sum_\vec{\alpha} J_\vec{\alpha} \Phi_\vec{\alpha}(\vec{\sigma}) = \sum_\vec{\alpha} J_\vec{\alpha} \prod_{(i,\nu) \in \vec{\alpha}} \varphi_\nu(\sigma_i)\]

where the $J_\vec{\alpha}$ are the coefficient of expansions. When the property $f(\vec{\sigma})$ corresponds to the energy, these coefficients are generally called effective cluster interactions (ECI).

Correlation Functions

If the property of interest is a scalar (e.g., the energy), it must be invariant under symmetry operations and so must be the cluster expansion. To insure the proper invariance, one can either add conditions on the ECI, or one can modify the cluster function to make them invariant under symmetry operations. In the second alternative, one can apply the Reynolds operator to the cluster functions. The symmeterization of the cluster function $\Phi_\vec{\alpha}$ yields the so-called correlation function $\Theta_A$ defined as:

\[\begin{split}\Theta_A(\vec{\sigma}) &= \frac{1}{|G|} \sum_{g \in G} g \cdot \Phi_\vec{\alpha}(\vec{\sigma}) \\ &= \frac{1}{|G|} \sum_{g \in G} \Phi_\vec{\alpha}(g \cdot \vec{\sigma}) \\ &= \frac{1}{|G|} \sum_{g \in G} \Phi_{g^{-1} \cdot \vec{\alpha}}(\vec{\sigma}) \\ &= \frac{1}{|G|} \sum_{\vec{\alpha'} \in A} \Phi_\vec{\alpha'}(\vec{\sigma}) \\ &= \frac{1}{|G|} \sum_{\vec{\alpha'} \in A} \prod_{(i,\nu) \in \vec{\alpha'}} \varphi_\nu(\sigma_i) \quad \text{where} \quad A = G \cdot \vec{\alpha} \quad \text{(the group action of $\vec{\alpha}$)}\end{split}\]

Correlation functions are simply the average values of all symmetrically equivalent cluster functions. As such, the correlation functions could be equivalently written as $\Theta_A(\vec{\sigma}) = \langle \prod_{\nu \in \vec{\alpha}} \varphi_\nu \rangle_A(\sigma)$.

_images/CE_formalism.png — Illustration of the computation of the correlation functions for a square lattice system composed of 4 sites with 2 species. The Chebyshev basis is used. There is a total of 16 ($2^4$) clusters, but due to the symmetry of the square lattice, there are only 6 different orbits (1 for the empty cluster, 1 for the point cluster, 2 for the pair clusters, 1 for the triplet cluster, and 1 for the quadruplet cluster) and as many correlations functions. This set of evaluated correlation function (i.e., [1 0 -1 1 0 1]) is unique and fully describes the arrangement illustrated (up to a symmetry operation). For each one of the 6 symmetrically unique clusters, the filled atomic sites represent atoms in the cluster, while transparent atomic sites are sites associated with the constant site-basis function.

Any scalar property that is a function of the arrangement on the lattice can be expanded as a linear combination of correlation functions as:

\[f(\vec{\sigma}) = \sum_A J_A \Theta_A(\vec{\sigma}) = \sum_A J_A \sum_{\vec{\alpha} \in A} \prod_{(i,\nu) \in \vec{\alpha}} \varphi_\nu(\sigma_i)\]

In practice, when inferring the energy, one truncates the expansion to consider compact clusters containing only few atoms. The rationalization behind it is that interaction energy depends on a small number of surrounding atoms and quickly decays as the distance of interaction increases.

Implementation

Given a set of $K$ site-basis functions $\{ \varphi_k \}$ and a cluster $C$ of body order $I$, one can construct all correlation functions $\Theta_\alpha$ associated to the cluster $C$ by applying the following algorithm:

List all $K^I$ possible cluster functions $\Phi_\alpha$.

cluster_functions = np.vstack(
         [
            arr.flatten()
            for arr in np.meshgrid(*[range(K) for i in range(|C|)], indexing="ij")
         ],
         dtype=int,
     ).T
Let’s, for instance, consider a triplet cluster $C$ with 2 non-constant site-basis function (the number of species $K=3$), the list of all possible cluster functions would yields:

\[\begin{split}[ \; &\varphi_2(\sigma_1) \cdot \varphi_2(\sigma_2) \cdot \varphi_2(\sigma_3), \\ &\varphi_2(\sigma_1) \cdot \varphi_2(\sigma_2) \cdot \varphi_3(\sigma_3), \\ &\varphi_2(\sigma_1) \cdot \varphi_3(\sigma_2) \cdot \varphi_2(\sigma_3), \\ &\varphi_3(\sigma_1) \cdot \varphi_2(\sigma_2) \cdot \varphi_2(\sigma_3), \\ &\varphi_2(\sigma_1) \cdot \varphi_3(\sigma_2) \cdot \varphi_3(\sigma_3), \\ &\varphi_3(\sigma_1) \cdot \varphi_2(\sigma_2) \cdot \varphi_3(\sigma_3), \\ &\varphi_3(\sigma_1) \cdot \varphi_3(\sigma_2) \cdot \varphi_2(\sigma_3), \\ &\varphi_3(\sigma_1) \cdot \varphi_3(\sigma_2) \cdot \varphi_3(\sigma_3) \; ]\end{split}\]

Form a list $\Gamma_C$ of lists $\gamma_\alpha$ containing all symmetrically-permutated cluster functions for all listed cluster functions $\Phi_\alpha$:

\[\gamma_\alpha = \Big\{ \pi(\Phi_\alpha) = \prod_{(i,\nu) \in \alpha} \varphi_\nu(\sigma_{\pi(i)}) = \Phi_{\alpha'} \; : \; \forall \pi \in \{ \pi \} \Big\}\]
permuted_cluster_functions = cluster_functions[:, permutations]
Let’s, for instance, consider that the sites 2 and 3 in the triplet cluster $C$ above are symmetrically equivalent within $C$, then $\exists g \in G_C$ such that the sites 2 and 3 are swapped. The list $\Gamma_C$ would yields:

\[\begin{split}[ \; &[ \;\varphi_2(\sigma_1) \cdot \varphi_2(\sigma_2) \cdot \varphi_2(\sigma_3), \; \varphi_2(\sigma_1) \cdot \varphi_2(\sigma_3) \cdot \varphi_2(\sigma_2) \;], \\ &[ \;\varphi_2(\sigma_1) \cdot \varphi_2(\sigma_2) \cdot \varphi_3(\sigma_3), \; \varphi_2(\sigma_1) \cdot \varphi_2(\sigma_3) \cdot \varphi_3(\sigma_2) \;], \\ &[ \;\varphi_2(\sigma_1) \cdot \varphi_3(\sigma_2) \cdot \varphi_2(\sigma_3), \; \varphi_2(\sigma_1) \cdot \varphi_3(\sigma_3) \cdot \varphi_2(\sigma_2) \;], \\ &[ \;\varphi_3(\sigma_1) \cdot \varphi_2(\sigma_2) \cdot \varphi_2(\sigma_3), \; \varphi_3(\sigma_1) \cdot \varphi_2(\sigma_3) \cdot \varphi_2(\sigma_2) \;], \\ &[ \;\varphi_2(\sigma_1) \cdot \varphi_3(\sigma_2) \cdot \varphi_3(\sigma_3), \; \varphi_2(\sigma_1) \cdot \varphi_3(\sigma_3) \cdot \varphi_3(\sigma_2) \;], \\ &[ \;\varphi_3(\sigma_1) \cdot \varphi_2(\sigma_2) \cdot \varphi_3(\sigma_3), \; \varphi_3(\sigma_1) \cdot \varphi_2(\sigma_3) \cdot \varphi_3(\sigma_2) \;], \\ &[ \;\varphi_3(\sigma_1) \cdot \varphi_3(\sigma_2) \cdot \varphi_2(\sigma_3), \; \varphi_3(\sigma_1) \cdot \varphi_3(\sigma_3) \cdot \varphi_2(\sigma_2) \;], \\ &[ \;\varphi_3(\sigma_1) \cdot \varphi_3(\sigma_2) \cdot \varphi_3(\sigma_3), \; \varphi_3(\sigma_1) \cdot \varphi_3(\sigma_3) \cdot \varphi_3(\sigma_2) \;] \; ]\end{split}\]

Form a list $\tilde{\Gamma}_C$ of lists $\tilde{\gamma}_\alpha$ of unique symmetrically-permuted cluster functions. The task is performed by finding all unique cluster functions and then group them together if there exists a symmetric permutation of atomic sites from one cluster function to the other.

\[\begin{split}&\text{If} \quad \Phi_\alpha, \Phi_\beta \in \tilde{\gamma}: \\ &\quad \text{then:} \quad \exists \pi \quad \text{such that} \quad \pi(\Phi_\alpha) = \prod_{(i,\nu) \in \alpha} \varphi_\nu(\sigma_{\pi(i)}) = \Phi_{\beta}\end{split}\]
# Remove duplicate cluster functions
(cluster_functions, inverse_indexes) = np.unique(
    permuted_cluster_functions.reshape(-1, |C|),
    return_inverse=True,
    axis=0,
)
# Remove duplicate permutations
inverse_indexes = inverse_indexes.reshape(
    permuted_cluster_functions.shape[:-1]
)
unique_permutations_indexes = np.unique(
    np.sort(inverse_indexes, axis=-1), axis=0
)

# Fill matrices with unique symmetrically-permuted cluster functions
for permutation_indexes in unique_permutations_indexes:
    indexes = np.unique(permutation_indexes)
    unique_cluster_functions = cluster_functions[indexes]
Continuing with the previous example, the list $\tilde{\Gamma}_C$ of $\tilde{\gamma}_\alpha$ would reads:

\[\begin{split}[ \; &[ \;\varphi_2(\sigma_1) \cdot \varphi_2(\sigma_2) \cdot \varphi_2(\sigma_3) \;], \\ &[ \;\varphi_2(\sigma_1) \cdot \varphi_2(\sigma_2) \cdot \varphi_3(\sigma_3), \; \varphi_2(\sigma_1) \cdot \varphi_2(\sigma_3) \cdot \varphi_3(\sigma_2) \;], \\ &[ \;\varphi_3(\sigma_1) \cdot \varphi_2(\sigma_2) \cdot \varphi_2(\sigma_3) \;], \\ &[ \;\varphi_2(\sigma_1) \cdot \varphi_3(\sigma_2) \cdot \varphi_3(\sigma_3) \;], \\ &[ \;\varphi_3(\sigma_1) \cdot \varphi_2(\sigma_2) \cdot \varphi_3(\sigma_3), \; \varphi_3(\sigma_1) \cdot \varphi_2(\sigma_3) \cdot \varphi_3(\sigma_2) \;], \\ &[ \;\varphi_3(\sigma_1) \cdot \varphi_3(\sigma_2) \cdot \varphi_3(\sigma_3) \;] \; ]\end{split}\]

Construct the correlation functions $\Theta_A$ from the unique groups $\tilde{\gamma}_\alpha$.

\[\Theta_A = \frac{1}{|G \cdot C|} \sum_{g \in G \cdot C} \frac{1}{\sqrt{|\tilde{\gamma}_\alpha|}} \sum_{\vec{\alpha} \in \tilde{\gamma}_\alpha} \prod_{(i, \nu) \in g \cdot \vec{\alpha}} \varphi_\nu(\sigma_i) \; , \quad \forall \tilde{\gamma}_\alpha \in \tilde{\Gamma}_C\]