# Expanding Density-Correlation Machine Learning Representations for
Anisotropic Coarse-Grained Particles

Arthur Lin  Department of Chemical and Biological Engineering, University of
Wisconsin, Madison, WI, USA  Kevin K. Huguenin-Dumittan  Laboratory of
Computational Science and Modeling, École Polytechnique Fédérale de Lausanne,
Lausanne, Switzerland  Yong-Cheol Cho  Department of Chemical and Biological
Engineering, University of Wisconsin, Madison, WI, USA  Department of Computer
Science and Engineering, University of Wisconsin, Madison, WI, USA  Jigyasa
Nigam  Laboratory of Computational Science and Modeling, École Polytechnique
Fédérale de Lausanne, Lausanne, Switzerland  Rose K. Cersonsky  Department of
Chemical and Biological Engineering, University of Wisconsin, Madison, WI, USA

###### Abstract

Physics-based, atom-centered machine learning (ML) representations have been
instrumental to the effective integration of ML within the atomistic
simulation community. Many of these representations build off the idea of
atoms as having spherical, or isotropic, interactions. In many communities,
there is often a need to represent groups of atoms, either to increase the
computational efficiency of simulation via coarse-graining or to understand
molecular influences on system behavior. In such cases, atom-centered
representations will have limited utility, as groups of atoms may not be well-
approximated as spheres. In this work, we extend the popular Smooth Overlap of
Atomic Positions (SOAP) ML representation for systems consisting of non-
spherical anisotropic particles or clusters of atoms. We show the power of
this anisotropic extension of SOAP, which we deem AniSOAP , in accurately
characterizing liquid crystal systems and predicting the energetics of Gay-
Berne ellipsoids and coarse-grained benzene crystals. With our study of these
prototypical anisotropic systems, we derive fundamental insights on how
molecular shape influences mesoscale behavior and explain how to reincorporate
important atom-atom interactions typically not captured by coarse-grained
models. Moving forward, we propose AniSOAP as a flexible, unified framework
for coarse-graining in complex, multiscale simulation.

††preprint: AIP/123-QED

##  I Introduction

In understanding molecular interactions and modeling their resultant behavior,
it is very often a worthwhile endeavor to group (i.e. coarse-grain) one or
more atoms as a theoretical “unit” or particle. This choice can be practical;
often, the time- and length-scales necessary to simulate molecular processes
limit our ability to simulate with all-atom resolution. Conversely, choosing
particle-based, rather than atom-based computational approaches can also be a
scientific choice; when we selectively limit the degrees of freedom within our
systems, we can identify factors that are causal to molecular behavior or
phenomena. Thus, these simplified coarse-grained simulation approaches serve
as both a tool and a lens with which to study chemical systems.

Similarly, machine learning (ML) methods have emerged as a powerful tool for
scientific inquiry, with the ability to elucidate new patterns within or
relationships between chemical spaces and observed properties, often in order
to predict the properties of unseen systems. So, how do we incorporate the
idea of atom grouping in the context of machine learning? Many approaches
consider our configurational space as a manifold and use a variety of machine
learning architectures (variationalMonroe and Shen (2022a, b) and
hierarchicalCosta  _et al._ (2023) auto-encoders, different forms of neural
networksMajewski  _et al._ (2023); Ruza  _et al._ (2020); Husic  _et al._
(2020); Wang  _et al._ (2019); Zhang  _et al._ (2018), ensemble learning Wang
_et al._ (2020)) to determine a latent space in which this grouping is
implicitly embedded that minimizes a chosen fitness function. Despite high
performance and generalizability, these end-to-end models are often limited by
data requirements or in fundamental analyses by their lack of intrinsic
interpretability.

Another approach within the machine-learning community is so-called _feature-
forward_ modeling, wherein one explicitly transforms the raw chemical data
into numerical “features” that reflect the underlying physics or chemistry of
interest prior to applying ML methods. While both end-to-end and feature-
forward approaches have merits and overlaps, the latter method is often noted
for its interpretability and comparable performance using shallower ML
architecturesDe  _et al._ (2016); Deringer  _et al._ (2021); Musil  _et al._
(2021a). Within this umbrella of approaches, there are many ways to encode the
raw chemical data into features, and the suitable choice depends entirely on
the scientific context. For cheminformatics, wherein we are often looking to
compare different chemistries or identify the role of specific functional
groups, string-based featurizations such as SMILESWeininger (1988) or
SELFIESKrenn  _et al._ (2020) are popular, as they encode important parameters
such as present functional groups and connectivity and can be parsed using
natural language processing (NLP) models and other text-based technologies.
However, in thermodynamic contexts where the chemistry and connectivity remain
largely unchanged, such as in molecular simulation, it is more typical to use
configuration-dependent featuresMusil  _et al._ (2021a); Wigh, Goodman, and
Lapkin (2022), including Behler-Parinello symmetry functionsBehler (2011),
smooth overlap of atomic positions (SOAP)Bartók, Kondor, and Csányi (2013),
and molecular graphs Wu  _et al._ (2018).

So, for coarse-grained systems, how do we apply a feature-forward approach? A
large challenge is that groups of atoms, hereon referred to as “particles”,
are non-isotropic, and so methods based on atomistic ML will fail to capture
the anisotropy of interparticle interactions. Frameworks based on density
expansions (e.g., SOAPBartók, Kondor, and Csányi (2013), NICENigam,
Pozdnyakov, and Ceriotti (2020)) present a compelling avenue for expansion,
given that they can putatively be made flexible to any density expansion, even
anisotropic density fields or hard particle volumes. Furthermore, by extending
these frameworks to isolate molecule-level interactions within atom-atom site
potentials, we gain the ability to combine or compare representations across
multiple scales.

Here, we propose and demonstrate the first such anisotropic expansion of
symmetrized density-based frameworks for ML representations by taking the
popular SOAP (Smooth Overlap of Atomic Positions) formalism and demonstrating
its expansion to simple anisotropic bodies, in a representation we deem
AniSOAP , which can be read as either “The Smooth Overlap of Anisotropic
Particles” or “Anisotropic-SOAP”. While here we demonstrate the expansion for
multivariate Gaussian densities, similar expansions can be made for arbitrary
anisotropic density fields.

In Section II, we provide the underlying mathematical theory behind this
expansion, including how deliberate selection of basis sets enables analytical
evaluation of the expansion coefficients. In Section III, we demonstrate three
performance-defining case studies for its usage across multiple simulation
length scales and materials systems. We start by analyzing two classic cases
of explicitly ellipsoidal particles: liquid crystalline configurations and
those governed by the Gay-Berne interaction potential. In doing so, we show
the conceptual overlap and divergence of AniSOAP from traditional mesoscale
order parameters, as well as its ability to generalize to supervised tasks.
Then, we analyze a set of systems that is only implicitly ellipsoidal –
benzene molecules arranged in stable and unstable crystalline configurations,
and show that AniSOAP can be used to provide molecule-level approximations to
first-principles energetics. We then couple AniSOAP with the traditional SOAP
formalism to demonstrate how to combine representations across multiple
scales. The corresponding open-source code, AniSOAP , is available at
github.com/cersonsky-lab/anisoap.

##  II Theory

General Variables
---
$\bm{r}$ |  a vector in Cartesian space
$r$ |  magnitude of vector $\bm{r}$
$\hat{r}$ |  direction of vector $\bm{r}$
$\bm{r}_{ij}$ |  the vector between point $i$ and point $j$
$\mathcal{R}$ |  a 3x3 rotation matrix
$\bm{x}$ |  a feature vector for one configuration
$\bm{X}$ |  a matrix containing, as rows, feature vectors for multiple configurations
Density Expansions
$g(\bm{r})$ |  a potentially anisotropic density function centered at $\bm{r}$
$\rho_{ij}$ |  the contribution of the density located at the position of point $j$ to the density expansion of point $i$
$n,l,m$ |  indices for radial bases and spherical harmonics
$R_{nl}$ |  a basis function indexed by $n$ and $l$
$\sigma_{\text{GTO}}$ |  the width of a Gaussian-type orbital (GTO) basis
$Y_{l}^{m}$ |  the $m^{th}$ component spherical harmonic of order $l$
$\mathcal{D}_{mm^{\prime}}^{l}$ |  the Wigner matrix used to rotate spherical harmonics
Multivariate Gaussians (MVG)
$\sigma_{1}$, $\sigma_{2}$, $\sigma_{3}$ |  the three principal components of an MVG
$\bm{D}$ |  the principal axis decomposition of an MVG, where $\bm{D}\equiv\mathrm{diag}\left(\frac{1}{\sigma_{1}^{2}},\frac{1}{\sigma_{2}^{2}},\frac{1}{\sigma_{3}^{2}}\right)$
$\mathbf{A}$ |  the rotated principal axis decomposition of an MVG, where $\mathbf{A}\equiv\mathcal{R}\bm{D}\mathcal{R}^{T}$
$\bm{S}$ |  the Gay-Berne analog to $\bm{D}$, where $\bm{S}\equiv\bm{D}^{-2}$
Table 1: Notation Guide. Throughout the text, we adopt the notation typical of
the atom-centered symmetrized representation community as detailed in the
table above.

In the traditional formalism for SOAP and related representations, we treat
each atom as a localized isotropic field in three dimensions and construct an
“atom-density" by summing over all neighboring atoms within a spherical shell.
The atom-density is usually written on a basis of radial and angular
functions, the latter of which are chosen to be spherical harmonics to make
the atom-density more easily amenable to rotational symmetrization.. To
represent the many-body nature of these atoms, we can then compute n-body
correlations of these density expansions or introduce message-passing into the
featurization formalismNigam  _et al._ (2022); Batatia  _et al._ (2022).

###  II.1 The Effect of Symmetry Breaking in Density Frameworks

Using the braket notation prescribed in Ref. 22, the contribution of an atom
$j$ to center atom $i$’s expansion coefficients are given by

$\displaystyle\braket{nlm}{\rho_{ij}}=\int_{\mathbb{R}^{3}}g(\bm{r}-\bm{r}_{ij})R_{nl}(r)Y_{l}^{m}(\hat{r})d^{3}\bm{r}.$
(1)

where $g(\bm{r})$ is a localized function (usually a gaussian), $R_{nl}(r)$ is
a radial basis function, and $Y_{l}^{m}(\hat{r})$ are spherical harmonics. For
a general $g(\bm{r})$ and $R_{nl}(r)$, there is no hope to evaluate this
analytically. There is not even a general way to evaluate a general one-
dimensional integral $\int f(x)\mathrm{d}x,$ which is why, after all, there
are so many books on integral tables.

If we use any arbitrary density $g(\bm{r})$, e.g. to model more closely the
shape of rigid molecules, nanoparticles, or arbitrary bodies, we are forced to
use numerical integration, which could be done using Lebdev grids. A fully
numerical implementation would also provide us with complete freedom regarding
the choice of the radial basis function $R_{nl}$, which would allow us to
choose the basis based on nice mathematical properties like the Laplacian
eigenstate (LE) basisBigi  _et al._ (2022). The main downside of the numerical
approach, at least for sufficiently general densities $g(\bm{r})$, is the
inevitable and potentially severe numerical cost-accuracy trade-off. The
necessary inaccuracies introduced by numerical integration may negate the
fidelity of the anisotropic density field, thus, we would therefore like to
examine the possibility of a fully analytical approach for the evaluation of
these coefficients.

For isotropic representations, it is traditional to do an implicit
transformation to align $\bm{r}_{ij}$ onto the $z-axis$ of our coordinate
system, eliminating any dependence of $g(\bm{r}-\bm{r}_{ij})$ on the angular
integrands. This simplifies (1), even for complicated bases, to

$\displaystyle\braket{nlm}{\rho_{ij}}=F(r_{ij})Y_{l}^{m}(\hat{r}_{ij}),$ (2)

where $F(r_{ij})$ is some function that contains the complete dependence on
the distance between the two atoms. The angular dependence is fully captured
in the spherical harmonic factor, allowing us to precompute $F(\bm{r}_{ij})$
on the full relevant interval for efficient spline evaluationsMusil  _et al._
(2021b).

For non-isotropic densities, however, this is no longer the case. Even taking
the relatively simple case of a multivariate Gaussian, we cannot decouple our
radial (e.g., distance) and angular components, as the Gaussian requires a
transformation into the three principal components
$\sigma_{1},\sigma_{2},\sigma_{3}$ and corresponding axes. In order to
maintain consistency with (1), we replace all instances of $\bm{r}$ with
$\mathcal{R}\mathcal{R}^{T}\bm{r}$, where $\mathcal{R}$ is our rotation matrix
into this component axis. We then simplify by changing our axes of integration
from $\bm{r}$ to $\bm{r}^{\prime}\equiv\mathcal{R}^{T}\bm{r}$ and noting that
$r=\mathcal{R}\mathcal{R}^{T}r$ to obtain

$\braket{nlm}{\rho_{ij}}=\int_{\mathbb{R}^{3}}\mathrm{d}^{3}(\bm{r}^{\prime})~{}g(\mathcal{R}\bm{r}^{\prime}-\mathcal{R}\bm{r}_{ij}^{\prime})R_{nl}(r)Y_{l}^{m}(\mathcal{R}\hat{r}^{\prime}).$

To perform the rotation of our spherical harmonics, one commonly sums over
Wigner-D matricesWigner (1931) to obtain

$\displaystyle\braket{nlm}{\rho_{ij}}$
$\displaystyle=\sum_{m^{\prime}=-l}^{l}\mathcal{D}^{l}_{mm^{\prime}}(\mathcal{R})\int_{\mathbb{R}^{3}}\mathrm{d}^{3}(\bm{r}^{\prime})~{}g(\mathcal{R}\left(\bm{r}^{\prime}-\bm{r}_{ij}^{\prime}\right))R_{nl}(r^{\prime})Y_{l}^{m\prime}(\hat{r}^{\prime})$
(3)

With this step, we can reduce our calculation in the special case in which the
matrix $\mathcal{D}^{l}_{mm^{\prime}}(\mathcal{R})$ is diagonal, which proves
useful for analytical evaluation techniques. We thus focus on the factor

$\displaystyle\int_{\mathbb{R}^{3}}\mathrm{d}^{3}\bm{r}^{\prime}~{}g(\mathcal{R}\left(\bm{r}^{\prime}-\bm{r}_{ij}^{\prime}\right))R_{nl}(r^{\prime})Y_{l}^{m\prime}(\hat{r^{\prime}}).$
(4)

which will be the focus of our discussion from hereon.

###  II.2 Multivariate Gaussian Densities

The simplest anisotropic function that can be put in place of $g(r)$ in (4) is
the multivariate Gaussian density (hereon MVG). The goal of this subsection is
to show that by choosing $g(\bm{r})$ to be an MVG and $R_{nl}(r)$ to be of the
monomial, GTO, or STO form, we can reduce the evaluation of the expansion
coefficient to the evaluation of integrals of the form

$\begin{split}\braket{nlm}{\rho_{ij}}=\int_{\mathbb{R}^{3}}e^{-\frac{1}{2}(\bm{r}-\bm{r}_{ij})^{T}\bm{D}(\bm{r}-\bm{r}_{ij})}\mathrm{p}(\bm{r}),\end{split}$
(5)

where $\mathrm{p}(\bm{r})$ is some polynomial expression in
$\bm{r}\equiv(x,y,z)$ and $\bm{D}$ is some $3\times 3$ diagonal matrix.

Consider the three-dimensional Gaussian defined by

$\displaystyle g:\mathbb{R}^{3}$ $\displaystyle\rightarrow$ (6)
$\displaystyle\bm{r}$ $\displaystyle\rightarrow
g(\bm{r})=\exp\left(-\frac{1}{2}\bm{r}^{T}\mathbf{A}\bm{r}\right),$ (7)

where $A$ is a symmetric and (strictly) positive definite $3\times 3$ matrix.
Any matrix satisfying these two conditions can be orthogonally diagonalized,
also called the principal axis decomposition, allowing us to write

$\displaystyle\mathbf{A}=\mathcal{R}\bm{D}\mathcal{R}^{T},$ (8)

where $\mathcal{R}\in\mathrm{SO}(3)$ is a rotation matrix that specifies the
three principal axes and

$\displaystyle\bm{D}\equiv\mathrm{diag}\left(\frac{1}{\sigma_{1}^{2}},\frac{1}{\sigma_{2}^{2}},\frac{1}{\sigma_{3}^{2}}\right)=\begin{pmatrix}\frac{1}{\sigma_{1}^{2}}&&\\\
&\frac{1}{\sigma_{2}^{2}}&\\\ &&\frac{1}{\sigma_{3}^{2}}\end{pmatrix}$ (9)

is a diagonal matrix containing the widths $\sigma_{j}$ of the Gaussian along
the three principal directions. With this decomposition, we can write

$\displaystyle g(\mathcal{R}\bm{r}^{\prime})$
$\displaystyle=\exp\left(-\frac{1}{2}\bm{r}^{\prime
T}\mathcal{R}^{T}\mathcal{R}D\mathcal{R}^{T}\mathcal{R}\bm{r}^{\prime}\right)$
(10) $\displaystyle=:\exp\left(-\frac{1}{2}\bm{r}^{\prime
T}\bm{D}\bm{r}^{\prime}\right),$ (11)

where $\bm{r}^{\prime}=\mathcal{R}^{T}\bm{r}$ are the coordinates with respect
to the principal axes.

Putting the equation for our MVG into (4), we get

$\displaystyle\int_{\mathbb{R}^{3}}\mathrm{d}^{3}\bm{r}^{\prime}\exp\left(-\frac{1}{2}(\bm{r}^{\prime}-\bm{r}_{ij}^{\prime})^{T}\bm{D}(\bm{r}^{\prime}-\bm{r}_{ij}^{\prime})\right)R_{nl}(r^{\prime})Y_{l}^{m\prime}(\hat{r^{\prime}}).$
(12)

We have now set the stage for the evaluation of the general coefficient. Our
goal is to choose a suitable radial basis that would allow us to evaluate the
coefficients analytically. One possibility is to use certain polynomial bases,
which is motivated by the fact that for any $n=0,1,2,\dots$, we have

$\displaystyle\int\mathrm{d}xx^{n}e^{-\frac{x^{2}}{2\sigma^{2}}}=\begin{cases}\left(2\sigma^{2}\right)^{\frac{n+1}{2}}\Gamma\left(\frac{n+1}{2}\right)&n\quad\mathrm{even}\\\
0&n\quad\mathrm{odd}\end{cases}$ (13)

where $\Gamma$ is the gamma function. If we could reduce the evaluation of the
integral (5) to such (multivariate) polynomials, this might provide us with
analytical expressions for the expansion coefficients. Thus, following (12),
we require that $R_{nl}(r)Y_{l}^{m}(\hat{r})$ is a polynomial.

##### Integrability of Monomial Basis

We will first show that when $R_{nl}(r)=r^{l+2n}$, the expansion coefficients
can be evaluated analytically. First, we separate $R_{nl}(r)$ into two
factors, where

$r^{l+2n}=(r^{2})^{n}\cdot r^{l}$ (14)

Given that $r^{l}Y_{l}^{m}(\hat{r})$ is a polynomial in $(x,y,z)$, as is
$r^{2}=x^{2}+y^{2}+z^{2}$ and thus, so is $(r^{2})^{n}$, the monomial basis
can be analytically integrated.

##### Integrability of GTO Basis

With one extra step, we can also show that a suitable modification of the GTO
basis, $R_{nl}(r)=r^{l+2n}e^{-\frac{r^{2}}{2\sigma^{2}}}$, provides an equally
well-suited basis. . By separating out all the exponential and polynomial
factors, we can write the integrand as

$\begin{split}&e^{-\frac{1}{2}(\bm{r}^{\prime}-\bm{r}_{ij}^{\prime})^{T}D(\bm{r}^{\prime}-\bm{r}_{ij}^{\prime})}R_{nl}(r^{\prime})Y_{l}^{m}(\hat{r^{\prime}})\\\
=&e^{-\frac{1}{2}(\bm{r}^{\prime}-\bm{r}_{ij}^{\prime})^{T}D(\bm{r}^{\prime}-\bm{r}_{ij}^{\prime})}r^{l+2n}e^{-\frac{r^{2}}{2\sigma^{2}}}Y_{l}^{m}(\hat{r})\\\
=&e^{-\frac{1}{2}(\bm{r}^{\prime}-\bm{r}_{ij}^{\prime})^{T}D(\bm{r}^{\prime}-\bm{r}_{ij}^{\prime})-\frac{r^{2}}{2\sigma^{2}}}r^{l+2n}Y_{l}^{m}(\hat{r})\\\
=&e^{-\frac{1}{2}(\bm{r}^{\prime}-\bm{r}_{ij}^{\prime})^{T}D(\bm{r}^{\prime}-\bm{r}_{ij}^{\prime})-\frac{r^{2}}{2\sigma^{2}}}\mathrm{p}(x,y,z)\\\
=&e^{-\mathrm{quadratic}(\bm{r})}\mathrm{p}(x,y,z)\end{split}$ (15)

We can see that this is again an exponential of a (now different) quadratic
form multiplied by a polynomial in $(x,y,z)$, which can be integrated
analytically.

Concluding this section, we can see that for both choices of the basis, the
final expression we need to evaluate is of the form

$\displaystyle\int_{\mathbb{R}^{3}}\mathrm{d}^{3}\bm{r}e^{-\frac{1}{2}(\bm{r}-\bm{r}_{ij})^{T}\bm{D}(\bm{r}-\bm{r}_{ij})}r^{2n}R_{l}^{m}(\hat{r}),$
(16)

with $R_{l}^{m}(\hat{r})=\bm{r}^{l}Y_{l}^{m}(\hat{r})$ and where the diagonal
matrix $\bm{D}$ is either the one directly obtained from the principal
components of the Gaussian density (monomial basis), or the modified version
(GTO basis). Similar arguments can be made for Slater-type orbitals (STO),
which follow a similar form to the GTO basis. From this expansion, we can then
calculate n-body correlations for a given particle, which for the 3-body term
is

$\displaystyle\braket{nl;n^{\prime}l^{\prime}}{\rho_{i}^{\nu=2}}=\sum_{m}\braket{nlm}{\rho_{i}}\braket{n^{\prime}l^{\prime}m}{\rho_{i}}$
(17)

where $\braket{nlm}{\rho_{i}}=\sum_{j}\braket{nlm}{\rho_{ij}}.$ In Appendix.
A, we have included a discussion of the practical steps toward computing the
expansion coefficient, including Löwdin symmetric orthogonalizationLöwdin
(1970) (discussed in Appendix. A.8) to impose orthonormality within the basis
set.

##  III Results and Discussion

Even for general point clouds, determining the optimal anisotropic analog is
non-trivial, and the focus of ongoing workSatorras, Hoogeboom, and Welling
(2022). Thus, we choose to demonstrate the efficacy of AniSOAP on systems
where the choice of an ellipsoidal proxy is trivial, being both historically
founded and well-defined. We ground our expansion in the rich history of
identifying causal mechanisms through anisotropic proxy particlesZhang  _et
al._ (2005), as many complex behaviors within molecular systems can be
explained by analyzing the steric interactions of their molecular
volumesGlotzer and Solomon (2007); Cersonsky  _et al._ (2018a, b); Keys,
Iacovella, and Glotzer (2011); Van Anders  _et al._ (2014). Future studies
will focus on optimizing AniSOAP for less trivial cases and further varieties
of molecular anisotropy.

Figure 1: The feature space of AniSOAP in the context of LC phases. The middle
plot shows the first two components of the structure-average AniSOAP vectors
for 2000 generated liquid crystal configurations showing either no
translational order (nematic phases, blue), planar stacking (smectic phases,
red), and varying degrees of orientation order (color saturation, with white
being little to no order). The inset shows the first two principal components
of the popular Steinhardt Order parameters, which clearly delineate the
smectic (red) phases, but in these and all following components, do not
separate the different nematic (blue) phases. Representative snapshots are
shown to the left and right corresponding to labeled points on the plots.

###  III.1 Ellipsoids in Liquid Crystals

In Sections III.1-III.2, we start by analyzing two archetypal ellipsoidal
systems: particles in different liquid crystalline (LC) phases, and particles
governed by the classical Gay-Berne potential. These case studies demonstrate
the similarities and differences between AniSOAP and traditional LC order
parameters, namely the orientational (or nematic) order parameter and pair-
wise Steinhardt order parameters, and examine the ability of AniSOAP to
perform continuous supervised tasks. In Section III.3, we then move to the
opposite end of the molecular spectrum and analyze benzene crystals –
molecules that are ellipsoidal in shape but whose atomic interactions are
complex and fundamental to their interaction landscape. Of the systems that
ellipsoidal AniSOAP is already well-suited to describe, these cases provide a
relative upper- and lower-bound to model performance from a system where
interactions are, by definition, entirely defined by the ellipsoidal
correlations to those where the energetics are minimally determined by the
shape of the molecule alone.

Unless otherwise specified, when discussing $X_{AniSOAP}~{}$, we are referring
to the power spectrum (the 3-body term) given in Eq. 17, but note the
simplicity of continuing onto higher body-order terms for greater accuracy and
transferabilityNigam, Pozdnyakov, and Ceriotti (2020).

Liquid crystals (LC) are mesoscopic phases that form from rigid ellipsoid-like
particles or molecules and can be controlled to direct the flow of light for a
variety of applications.Hamley (2007) Typically, LCs are characterized as
phases with orientational order but limited positional or translational order.
Molecules with orientational order but no translational order are generally
considered “nematic”, whereas those exhibiting one direction of translational
order (e.g. the particles appear as stacking planes, where within the planes,
there is little translational order) are deemed smectic.

In this analysis, we aim to show how our new AniSOAP representation is similar
to previously established descriptors and how it enables functionality and
analyses beyond the current capabilities. In characterizing LC phases,
scientists have relied on a library of different order parameters, including
different orientational order parameters Ramasubramani  _et al._ (2020) to
classify the orientational alignment of particles, Steinhardt order parameters
(SOPs Steinhardt, Nelson, and Ronchetti (1983)) to characterize the
neighborhood of different particles. We note the similarity of SOPs to Eq.
(2), wherein an SOP for a given particle involves the integration of spherical
harmonics over that particle’s neighbors:

$\braket{lm}{q_{i}}=\frac{1}{N_{j}}\sum_{j}w_{ij}Y_{l}^{m}(r_{ij})$ (18)

where $j$ is a neighbor of $i$, and weights $w_{ij}$ can be introduced based
on neighbor distances or Voronoi tesselationsLechner and Dellago (2008). The
$l^{th}$ SOP is computed by combining these terms such that

$X_{SOP}=\braket{l}{q_{i}}=\sqrt{\frac{4\pi}{2l+1}\sum_{m=-l}^{l}|\braket{lm}{q_{i}}|^{2}}.$
(19)

Thus, for comparison, for all configurations, we also compute the
orientational order parameter $\bm{X}_{OOP}$ (traditionally called the nematic
order parameter) and distance-weighted Steinhardt order parameters
$\bm{X}_{SOP}$ with freudRamasubramani  _et al._ (2020), using similar
$l_{\text{max}}$ and cutoff radius to ensure comparability. All
representations are properly scaled and centered – that is, each matrix of
feature vectors is centered to have zero mean with unit variance, where
$\bm{X}_{SOP}$ and $\bm{X}_{AniSOAP}~{}$ are scaled non-column-wise in order
to retain important relative variance information.

We generate 1000 liquid crystal configurations, 500 exhibiting nematic order,
and 500 exhibiting smectic order. Phases were generated to correspond with a
range of orientational order parameters, from 0 (particles are randomly
oriented) to 1 (particles are all ordered along the same director). We
populate each of these configurations with prolate ellipsoids with a length-
to-diameter ratio (L/D) equal to 3. For each of these configurations, we
compute the AniSOAP radial and power spectrum $\sigma_{GTO}=3.0$,
$n_{\text{max}}=6$, $l_{\text{max}}=6$ (for the power spectrum), and a cutoff
sufficient to include the surrounding neighbors of each ellipsoid, including
in neighboring smectic planes.

We first look at the principal components (PCs) of the AniSOAP vectors for
these phases (Fig. 1). The PCA reflects that the AniSOAP features can be used
to delineate translational order (blue versus red coloring), and orientational
order (saturated versus unsaturated color). These mappings highlight
qualitatively the information contained in the AniSOAP representation and show
that it simultaneously and smoothly represents the translational and
orientational order of the configurations. In Fig. S2, we demonstrate that the
AniSOAP features also smoothly delineate different particle shapes, with a
similar plot provided with data for both L/D=2 and 3 ellipsoids. By combining
this information into one smooth feature space, we are able to better recreate
the “nearsightedness” of molecular interactionsProdan and Kohn (2005), as
these aspects (shape, mutual orientation, and translational order) are often
interrelated and correlated in the ways in which they influence molecular
behavior.

Figure 2: Comparison of Different Featurizations for Liquid Crystal
Configurations. The Global Feature Reconstruction Error (GFRE) encodes the
error using one feature representation $\bm{X}_{\text{input}}$ to reconstruct
another $\bm{X}_{\text{target}}$, and can be used to infer the information
density of one featurization compared to another. Higher values indicate poor
reconstruction, whereas lower indicates better reconstruction. Here we see
that AniSOAP features better reconstruct traditional Steinhardt order
parameters ($q_{0-6}$) and orientational order parameters than vice versa.
While Steinhardt OPs carry large mutual information with the AniSOAP vectors,
they are unable to reconstruct the Orientational Order Parameter.

Quantitatively, we can compare the information density of AniSOAP with other
order parameters for our configurations using the global feature
reconstruction error (GFREGoscinski  _et al._ (2021a), computed using scikit-
matterGoscinski  _et al._ (2023)). For two representations $\bm{X}_{1}$ and
$\bm{X}_{2}$ of the same dataset, $\text{GFRE}(X_{1},X_{2})$, determines how
much information $\bm{X}_{1}$ contains relative to $\bm{X}_{2}$, where 0.0
infers $\bm{X}_{1}$ can perfectly reconstruct $\bm{X}_{2}$, and 1.0
corresponds to poor reconstruction. As shown in Fig. 2, we see that AniSOAP
representations (both the 2-body radial spectrum $X^{(\nu=1)}_{AniSOAP}~{}$
and 3-body power spectrum $X^{(\nu=2)}_{AniSOAP}~{}$) are able to reconstruct
the analogous $\bm{X}_{SOP}$ with lower error than vice versa (6-8% versus
12-14% reconstruction error). Furthermore, the AniSOAP representations are
able to reconstruct the orientational order parameter with low error (6% for
the 3-body, and 24% for the 2-body, compared to 39% for $X_{SOP}$). The
improvement in information density can be attributed to the smooth radial
bases that underlie the AniSOAP construction. It is worth noting that the
Steinhardt order parameter is not well-suited for heterogenous datasets, as
shown in Fig. S3, and its ability to reconstruct the AniSOAP representation
will decrease when a dataset contains multiple particle types.

###  III.2 Gay-Berne Ellipsoids

The Gay-Berne (GB) potential is a Lennard Jones-type potential that contains
additional terms to account for ellipsoidal anisotropies. We use the
generalized formulation of the Gay-Berne potential outlined by Everaers and
Ejtehadi (2003) and first introduced by Berardi, Fava, and Zannoni (1995).
This formulation calculates the pairwise potential between (potentially
dissimilar) ellipsoids, i.e., ellipsoids with three unequal semi-axes:
$a_{i},b_{i},c_{i}$. These three semi-axes define a diagonal structure matrix
for particle i:

$\bm{S}_{i}=\begin{pmatrix}a_{i}&0&0\\\ 0&b_{i}&0\\\ 0&0&c_{i}\\\
\end{pmatrix}$ (20)

where $\bm{S}$ is analogous to $\bm{D}^{-2}$ of Eq. (11), although $(a,b,c)$
are discrete semiaxes, rather than Gaussian widths, like $\sigma$. In our
case-study, we generate 25,000 dimers of $a_{i}=1,b_{i}=1.5,c_{i}=2$
ellipsoids at random offsets and orientations.

The center position and orientation of ellipsoid $i$ is given by position
vector $\bm{r}_{i}$ and a 3$\times$3 rotation matrix $\mathcal{R}_{i}$,
respectively. The relative position is $\bm{r}_{12}=\bm{r}_{2}-\bm{r}_{1}$.
The Gay-Berne potential takes into account interparticle distance dependence
and orientation dependence of the ellipsoid pairs. The potential is given as a
product of three terms:

$\displaystyle
U(\mathcal{R}_{1},\bm{S}_{1},\mathcal{R}_{2},\bm{S}_{2},\vec{r_{12}})=\underbrace{U_{r}(\mathcal{R}_{1},\bm{S}_{1},\mathcal{R}_{2},\bm{S}_{2},\vec{r_{12}})}_{\text{LJ-
like term}}$ (21)
$\displaystyle\cdot\underbrace{\eta_{12}(\mathcal{R}_{1},\bm{S}_{1},\mathcal{R}_{2},\bm{S}_{2})\cdot\chi_{12}(\mathcal{R}_{1},\mathcal{R}_{2},\hat{r_{12}})}_{\text{anisotropy
corrections}}$ (22)

The first term $U_{r}$ is given as follows:

$U_{r}=4\epsilon_{GB}\left[\left(\frac{\sigma_{GB}}{h_{12}+\gamma\sigma_{GB}}\right)^{12}-\left(\frac{\sigma_{GB}}{h_{12}+\gamma\sigma_{GB}}\right)^{6}\right]$
(23)

where the distance $h_{12}$ between two ellipsoids is estimated as the Perram
distance of closest approachPerram and Wertheim (1985)

$h_{12}(\mathcal{R}_{1},\mathcal{R}_{2},\bm{S}_{1},\bm{S}_{2},\bm{r}_{12})=r_{12}-s_{12}(\mathcal{R}_{1},\mathcal{R}_{2},\bm{S}_{1},\bm{S}_{2},\hat{r}_{12})$
(24)

where

$s_{12}(\mathcal{R}_{1},\mathcal{R}_{2},\bm{S}_{1},\bm{S}_{2},\hat{r}_{12})=\left[\frac{1}{2}\hat{r}_{12}^{T}\bm{\Xi}_{12}^{-1}(\mathcal{R}_{1},\mathcal{R}_{2},\bm{S}_{1},\bm{S}_{2})\hat{r}_{12}\right]^{-1/2}$
(25)

and

$\bm{\Xi}_{12}(\mathcal{R}_{1},\mathcal{R}_{2},\bm{S}_{1},\bm{S}_{2})=\mathcal{R}_{1}^{T}\bm{S}_{1}^{2}\mathcal{R}_{1}+\mathcal{R}_{2}^{T}\bm{S}_{2}^{2}\mathcal{R}_{2}$
(26)

The terms $\epsilon_{GB}$ and $\sigma_{GB}$ are the potential well depth and
location analogous to LJ. $\gamma$ is a parameter that shifts the potential
well and is typically set to 1. Since the interparticle distance is defined to
be the “distance of closest approach”, $h_{12}$, rather than center-center
distance, so $U_{r}$ must also contain dependence on the geometry and
orientation of each ellipsoid to calculate $h_{12}$. The remaining pieces of
Eq. (22) account for the ellipsoidal geometry ($\eta_{12}$) and misorientation
($\chi_{12}$); further explanation of these terms are given in Sec. B.1.

One of the less trivial aspects of the Gay-Berne implementation is calculating
the “distance of closest approach” $h_{12}$. The approximation given by Perram
and Wertheim (1985) is computationally efficient, but open to failure in some
cases, like in the case of two bodies with very unequal radii Everaers and
Ejtehadi (2003). It is worth noting that it is the nature of density
expansions to implicitly contain the distance of closest approach, provided
that the chosen cutoff distance and radial resolution are sufficiently large
to include any necessary neighbors, which is evidenced by the perfect
agreement in the computed and predicted $h_{12}$ distances given in Fig.
3111The $h_{12}$ distances here fall within the confidence region of the
Perram and Wertheim (1985) algorithm, so they can be taken as exact. And thus,
in AniSOAP , we alleviate the requirement to explicitly compute $h_{12}$ or
rely on its approximations.

Figure 3: Predictions of ellipsoid interactions based on the Gay-Berne
Potential. The ellipsoid dimer dataset consists of two ellipsoids with varying
inter-center distance d and orientations based on rotations along the x
($\phi_{x}$), y ($\phi_{y}$), and z ($\phi_{z}$) axes. These varying
orientations and distances lead to different characteristic energy wells,
shown in the top right. Parity plots detail the performance in predicting the
distance of closest approach $h_{12}$ (bottom-left, units of characteristic
length scale $\sigma$) and energies (bottom-right, units of potential well
depth $\epsilon$). For clarity, only the test points are shown; the train
points show very similar trends.

From here, learning the Gay-Berne potential is straightforward, similar to how
traditional SOAP vectors can learn the LJ potential. With the AniSOAP power
spectrum vectors, it is possible to interpolate, using solely regularized
linear regression, the interaction potential of ellipsoids of arbitrary
distances from one another, as shown in Fig. 3. Performance decreases for
repulsive, heavily overlapped configurations (right end of the parity plot in
the lower right of Fig. 3), likely due to the fact that, in the repulsive
regime, small feature differences correspond to disproportionately large
energy differences. Our learning exercise converges with as few as 1,000
training points, as evidenced by the learning curve given in Fig. S4.

###  III.3 Benzene Molecules

Figure 4: Representations of the T-Shaped Benzene Dimer. From left to right. A
molecular image of the dimer at approximately 5.5Å separation, generated using
OvitoStukowski (2009). The underlying representation to the SOAP formalism,
wherein each atom is represented by a Gaussian field, and we expand n-body
terms of each atomic neighborhood. The underlying representation for the
AniSOAP formalism, wherein each _molecule_ is represented by an anisotropic
Gaussian field, and we expand n-body terms of each _molecular_ neighborhood.

For the following analysis, we intentionally construct a set of benzene
configurations to demonstrate the possibilities and, simultaneously, the
limitations for AniSOAP .

We start with a set of symmetry-constrained crystals comprised of planar
benzene molecules with the software PyXtalFredericks  _et al._ (2021) across
all 230 space groups. We then augmented this dataset by randomly rotating or
translating the molecules, where many configurations correspond to the same
_positions_ or _orientations_ of the molecules. The resultant dataset contains
roughly 7,000 benzene crystals. We then compute energetic quantities using
QuantumEspresso v7.0 Giannozzi  _et al._ (2009) using Perdew–Burke-Ernzerhof
(PBE) pseudopotentials and cutoff parameters reported by Prandini  _et al._
(2018), a Grimme D3-dispersion correctionGrimme  _et al._ (2016), and a 3x3
Monkhorst-Pack k-point gridMonkhorst and Pack (1976). Computations were
managed with the signac and signac-flow packagesRamasubramani  _et al._
(2018); Adorf  _et al._ (2018), and computed using the ASE computational
front-endHjorth Larsen  _et al._ (2017).

#### Hyperparameter Tuning

With any new featurization, there is always valid concern about hyperparameter
optimization and the sensitivity of the featurization to changes in these
valuesCheng  _et al._ (2020). In many ways, this is the appeal of the SOAP
formalism – the hyperparameters, namely the widths of the Gaussian densities,
basis sets, and cutoff parameters, all have roots in chemical physics and can
be chosen from a large range of “reasonable” values with minimal impact on
model interpretability and performance.

As the primary axes of planar benzenes are well-defined, the only new
hyperparameters to consider are the semiaxes lengths of the MVG. For
simplicity, we will only consider $\sigma_{1}=\sigma_{2}$, and prove that, in
line with conventional knowledge on benzene geometry, best results are
obtained with $\sigma_{1}>\sigma_{3}$ (an oblate ellipsoid).

To tune these parameters, we again use the GFRE of a given AniSOAP
representation and an analogous SOAP representation of our benzene crystal
configurations. By taking $\text{GFRE}(X_{\text{AniSOAP}},X_{\text{SOAP}})$,
we determine how much information is lost by moving from the atomistic to
coarse-grained AniSOAP representation, again where 0.0 infers perfect
reconstruction of the atomic correlations, and higher numbers correspond to
poor reconstruction. Note that this hyperparameter tuning occurs independently
of any traditional supervised learning task.

For our SOAP representation, we compute the structure-averaged representation
using rascalineFraux  _et al._ (2023), with a cutoff radius of $7.0$Å,
Gaussian density width of $0.5$Å, $l_{\text{max}}=10$, and $n_{\text{max}}=4$.
This results in a length $1,188$ vector for each configuration. We compare
against AniSOAP vectors using a similar cutoff radius, number of angular and
radial channels, varying $\sigma_{1}=\sigma_{2}$, $\sigma_{3}$, and
$\sigma_{\text{GTO}}$. Each AniSOAP vector has length $146$.

Figure 5: Results of tuning $\sigma_{1}$ and $\sigma_{3}$ for benzene
crystals. The color of the scatter points and contour levels denote
$\text{GFRE}(X_{\text{AniSOAP}},X_{\text{SOAP}})$, where lower values indicate
a greater fidelity of the AniSOAP representation to the SOAP analog, and
higher values indicate a greater information loss. Our MVG is oblate where
$\sigma_{1}>\sigma_{3}$, prolate where $\sigma_{1}<\sigma_{3}$, and spherical
where $\sigma_{1}=\sigma_{3}$. As expected, oblate representations carry more
mutual information to the atomistic representation than prolate or spherical
ones. The spherical representation $\sigma_{1}=\sigma_{2}=\sigma_{3}$ is
consistent with constructing a traditional SOAP representation from the
molecule centers.

From the results in Fig. 5, we see how the oblate MVG with
$\sigma_{1}=\sigma_{2}=$4.0 and $\sigma_{3}=$0.5 minimizes the information
loss compared to the atomistic SOAP representation with a $\sigma_{GTO}=1.5$.
Prolate ($\sigma_{1}<\sigma_{3}$) and spherical ($\sigma_{1}=\sigma_{3}$)
Gaussian densities perform notably worse. It is interesting, however, that
there is a wide range of $\sigma_{1},\sigma_{3}$ values that obtain similar
results, signifying that the AniSOAP hyperparameters are robust with respect
to small changes, provided the MVG remain oblate, and the semiaxis lengths are
within reasonable proportion. We include in the SI similar analyses for
$\sigma_{\text{GTO}}=\\{0.5,1.0,1.5,2.0,2.5,3.0\\}$, demonstrating that
$\sigma_{\text{GTO}}\in[1.0\text{\r{A}},3.0\text{\r{A}}]$ obtains similar
results to $\sigma_{\text{GTO}}=1.5$Å.

#### Learning of Benzene Energetics

Taking the optimized hyperparameters from the previous section, we now show
how the AniSOAP representation performs in simple supervised tasks,
demonstrating where AniSOAP performs best and where it is incomplete for
atomistic systems. It is worth noting that the intention of AniSOAP is to
provide a molecule-level analog to SOAP, and is, by design,
incompletePozdnyakov  _et al._ (2020) with respect to many atom-atom
correlations.

Figure 6: Parity plots comparing the performance of SOAP (built from the atom
centers), optimized SOAP (built on the molecule centers), and AniSOAP  for the
dataset of benzene crystalline configurations. All regressions were conducted
using 5-fold cross-validated regularized linear regression, and errors are
reported for the same 90/10 training/test split. The SOAP representations was
first reduced via PCA to be of similar rank to AniSOAP . In the parity plots,
color denotes the error in meV/atom, with the corresponding distribution of
errors shown in the lower panel.

To highlight the importance of molecular representation, we will use
regularized ridge models, employing a 90/10 training/test split and five-fold
cross-validation. We first perform three regression tasks – learning the
baselined, per-atom energy (equivalent to learning the per-molecule energy)
using the all-atom SOAP representation, a SOAP representation built solely
from the molecule centers (coinciding with the ridge line at
$\sigma_{1}=\sigma_{2}=\sigma_{3}$ in Fig. 5), and our optimized AniSOAP
representation. We have enforced that all three representations have similar
ranks, as higher-rank features will outperform lower-rank ones based on size
alone. To do so, we perform dimensionality reduction via principal components
analysis222principal covariates regressionHelfrecht  _et al._ (2020), while
ideal for regression-related dimensionality reductions, would again limit our
ability to compare the _relevant information content_ of these
representations..

The results, as shown in Fig. 6, demonstrate how incorporating the anisotropy
of intermolecular correlations can greatly improve our ability to learn
energetic quantities, even in contexts where atom-atom interactions are
important. A SOAP representation on the molecule centers (center panels)
333using the same major semiaxis length as the width of the isotropic Gaussian
field, and keeping all other hyperparameters consistent, has heavily limited
regression performance, with a typical RMSE on the order of $45-50$meV/atom,
and an $R^{2}\leq 0.4$. By introducing intermolecular anisotropy, our
performance jumps $R^{2}\approx 0.9,\text{RMSE}\approx 15-20$meV/atom, much
closer to the all-atom regression in the left panel of Fig. 6.

Figure 7: A parity plot showing the results of baselining our energies on the
AniSOAP approximation, followed by atom-level regression via SOAP. The inset
demonstrates the improvement in the distribution of errors.

Given that, by definition, AniSOAP does not contain the information on
atomistic correlations, it follows naturally that we can learn a “correction”
to the AniSOAP energies from our original SOAP representation, as corroborated
by Fig. 7. This is important for a few reasons – firstly, it shows that
AniSOAP can be used as a coarse-graining baseline to complex energetics.
Simply learn the energies due to body-body interactions, and incorporate an
atom-level correction where necessary. Secondly, it helps identify these
spaces where atom-atom interactions are driving the majority of energetic
quantities, which can be confirmed by leveraging the linear relationship
between the atom-level SOAP vector representation and the errors in the
AniSOAP regression and perform Principal Covariates Regression
(PCovR)Helfrecht  _et al._ (2020); Goscinski  _et al._ (2023). We have
included in the SI such a mapping, explorable via the online tool
chemiscopeFraux, Cersonsky, and Ceriotti (2020). We also posit that one can do
a combined AniSOAP and SOAP representation of a system to optimize the
segmentation of molecular energetics; however, this is beyond the scope of
this current study.

##  IV Conclusions

The tools developed under the umbrella of “machine-learning” are powerful, and
revolutionizing society and the scientific community at an aggressive pace.
One important advancement for chemical sciences is the concept of treating
chemistry as data – choosing specifically the numerical representations with
which to cast chemical questions into statistical and analytical models. These
representations are strongest when grounded in the established physical
principles that we are trying to emulate or predict.

Here, we have demonstrated one such physically-driven approach to machine-
learning representations, aimed at leveraging machine learning for coarse-
grained and mesoscale entities, demonstrating such an approach for particles
and molecules that are well-represented by ellipsoidal bodies. Our results
show that, for both classical and quantum mechanical datasets, AniSOAP
provides a suitable representation for ellipsoidal bodies, as these
representations are able to accurately and _linearly_ map onto complex
energetics without the need for deeper ML infrastructures. We have constructed
AniSOAP in such a way as to retain compatibility with the popular atomistic
SOAP formalism – in practice, we hope that these two technologies are used
hand-in-hand to simplify energetic landscapes and explicitly incorporate many-
body effects. Furthermore, this consistency between formalisms can enable
multiscale machine learning approaches, wherein one can use coarse-grained
representations for long-scale molecular motion and atom-atom representations
for nearsighted interactions. Future work will focus on the efficient
implementation of anisotropic decompositions for many-body effects and their
application to such systems.

##  V Data availability

Data for this paper will be made available via MaterialsCloudTalirz  _et al._
(2020) upon acceptance, including all raw input files and analysis scripts.

##  VI Author contributions

KKHD derived the expansion of multi-variate Gaussian and wrote the
corresponding proof-of-principle code. AL, JN, Y-CC and RKC refined the code,
implementing orthonormalization and additional basis sets. AL and RKC wrote
the manuscript and designed and executed the case studies.

##  VII Conflicts of interest

There are no conflicts to declare.

##  VIII Acknowledgements

This project was funded by the Wisconsin Alumni Research Fund (RKC), by NSF
through the University of Wisconsin Materials Research Science and Engineering
Center (DMR-2309000, AL), and the European Research Council (ERC) under the
research and innovation program (Grant Agreement No. 101001890-FIAMMA, JN,
KKHD).

We extend our un-ending gratitude to Guillaume Fraux and the developers of
rascaline for fielding our many questions during the implementation and
validation of AniSOAP .

## References

  * Monroe and Shen (2022a) J. I. Monroe and V. K. Shen, “Learning Efficient, Collective Monte Carlo Moves with Variational Autoencoders,” Journal of Chemical Theory and Computation 18, 3622–3636 (2022a).
  * Monroe and Shen (2022b) J. I. Monroe and V. K. Shen, “Systematic control of collective variables learned from variational autoencoders,” The Journal of Chemical Physics 157, 094116 (2022b).
  * Costa  _et al._ (2023) A. d. S. Costa, I. Mitnikov, M. Geiger, M. Ponnapati, T. Smidt, and J. Jacobson, “Ophiuchus: Scalable Modeling of Protein Structures through Hierarchical Coarse-graining SO(3)-Equivariant Autoencoders,”  (2023), arXiv:2310.02508 [cs].
  * Majewski  _et al._ (2023) M. Majewski, A. Pérez, P. Thölke, S. Doerr, N. E. Charron, T. Giorgino, B. E. Husic, C. Clementi, F. Noé, and G. De Fabritiis, “Machine learning coarse-grained potentials of protein thermodynamics,” Nature Communications 14, 5739 (2023), number: 1 Publisher: Nature Publishing Group.
  * Ruza  _et al._ (2020) J. Ruza, W. Wang, D. Schwalbe-Koda, S. Axelrod, W. H. Harris, and R. Gómez-Bombarelli, “Temperature-transferable coarse-graining of ionic liquids with dual graph convolutional neural networks,” The Journal of Chemical Physics 153, 164501 (2020).
  * Husic  _et al._ (2020) B. E. Husic, N. E. Charron, D. Lemm, J. Wang, A. Pérez, M. Majewski, A. Krämer, Y. Chen, S. Olsson, G. de Fabritiis, F. Noé, and C. Clementi, “Coarse graining molecular dynamics with graph neural networks,” The Journal of Chemical Physics 153, 194101 (2020).
  * Wang  _et al._ (2019) J. Wang, S. Olsson, C. Wehmeyer, A. Pérez, N. E. Charron, G. de Fabritiis, F. Noé, and C. Clementi, “Machine Learning of Coarse-Grained Molecular Dynamics Force Fields,” ACS Central Science 5, 755–767 (2019), publisher: American Chemical Society.
  * Zhang  _et al._ (2018) L. Zhang, J. Han, H. Wang, R. Car, and W. E, “DeePCG: Constructing coarse-grained models via deep neural networks,” The Journal of Chemical Physics 149, 034101 (2018).
  * Wang  _et al._ (2020) J. Wang, S. Chmiela, K.-R. Müller, F. Noé, and C. Clementi, “Ensemble learning of coarse-grained molecular dynamics force fields with a kernel approach,” The Journal of Chemical Physics 152, 194106 (2020).
  * De  _et al._ (2016) S. De, A. P. Bartók, G. Csányi, and M. Ceriotti, “Comparing molecules and solids across structural and alchemical space,” Physical Chemistry Chemical Physics 18, 13754–13769 (2016).
  * Deringer  _et al._ (2021) V. L. Deringer, A. P. Bartók, N. Bernstein, D. M. Wilkins, M. Ceriotti, and G. Csányi, “Gaussian Process Regression for Materials and Molecules,” Chemical Reviews 121, 10073–10141 (2021).
  * Musil  _et al._ (2021a) F. Musil, A. Grisafi, A. P. Bartók, C. Ortner, G. Csányi, and M. Ceriotti, “Physics-Inspired Structural Representations for Molecules and Materials,” Chemical Reviews 121, 9759–9815 (2021a), publisher: American Chemical Society.
  * Weininger (1988) D. Weininger, “SMILES, a chemical language and information system. 1. Introduction to methodology and encoding rules,” Journal of Chemical Information and Computer Sciences 28, 31–36 (1988), publisher: American Chemical Society.
  * Krenn  _et al._ (2020) M. Krenn, F. Häse, A. Nigam, P. Friederich, and A. Aspuru-Guzik, “Self-referencing embedded strings (SELFIES): A 100% robust molecular string representation,” Machine Learning: Science and Technology 1, 045024 (2020), publisher: IOP Publishing.
  * Wigh, Goodman, and Lapkin (2022) D. S. Wigh, J. M. Goodman, and A. A. Lapkin, “A review of molecular representation in the age of machine learning,” WIREs Computational Molecular Science 12, e1603 (2022), _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1603.
  * Behler (2011) J. Behler, “Atom-centered symmetry functions for constructing high-dimensional neural network potentials,” The Journal of Chemical Physics 134, 074106 (2011).
  * Bartók, Kondor, and Csányi (2013) A. P. Bartók, R. Kondor, and G. Csányi, “On representing chemical environments,” Physical Review B 87, 184115 (2013).
  * Wu  _et al._ (2018) Z. Wu, B. Ramsundar, E. N. Feinberg, J. Gomes, C. Geniesse, A. S. Pappu, K. Leswing, and V. Pande, “MoleculeNet: a benchmark for molecular machine learning,” Chemical Science 9, 513–530 (2018), publisher: Royal Society of Chemistry.
  * Nigam, Pozdnyakov, and Ceriotti (2020) J. Nigam, S. Pozdnyakov, and M. Ceriotti, “Recursive evaluation and iterative contraction of N-body equivariant features,” The Journal of Chemical Physics 153 (2020), publisher: AIP Publishing.
  * Nigam  _et al._ (2022) J. Nigam, S. Pozdnyakov, G. Fraux, and M. Ceriotti, “Unified theory of atom-centered representations and message-passing machine-learning schemes,” The Journal of Chemical Physics 156, 204115 (2022).
  * Batatia  _et al._ (2022) I. Batatia, D. P. Kovacs, G. Simm, C. Ortner, and G. Csanyi, “MACE: Higher Order Equivariant Message Passing Neural Networks for Fast and Accurate Force Fields,” Advances in Neural Information Processing Systems 35, 11423–11436 (2022).
  * Musil  _et al._ (2021b) F. Musil, M. Veit, A. Goscinski, G. Fraux, M. J. Willatt, M. Stricker, T. Junge, and M. Ceriotti, “Efficient implementation of atom-density representations,” The Journal of Chemical Physics 154, 114109 (2021b).
  * Bigi  _et al._ (2022) F. Bigi, K. K. Huguenin-Dumittan, M. Ceriotti, and D. E. Manolopoulos, “A smooth basis for atomistic machine learning,” The Journal of Chemical Physics 157, 234101 (2022).
  * Wigner (1931) E. Wigner,  _Gruppentheorie und ihre Anwendung auf die Quantenmechanik der Atomspektren_ (Vieweg+Teubner Verlag, Wiesbaden, 1931).
  * Löwdin (1970) P.-O. Löwdin, “On the Nonorthogonality Problem**The work reported in this paper has been sponsored in part by the Swedish Natural Science Research Council, in part by the Air Force Office of Scientific Research (OSR) through the European Office of Aerospace Research (OAR), U.S. Air Force under Grant AF-EOAR 67-50 with Uppsala University, and in part by the National Science Foundation under Grant GP-5419 with the University of Florida.” in  _Advances in Quantum Chemistry_, Vol. 5, edited by P.-O. Löwdin (Academic Press, 1970) pp. 185–199.
  * Satorras, Hoogeboom, and Welling (2022) V. G. Satorras, E. Hoogeboom, and M. Welling, “E(n) Equivariant Graph Neural Networks,”  (2022), arXiv:2102.09844 [cs, stat].
  * Zhang  _et al._ (2005) Zhang, A. S. Keys, T. Chen, and S. C. Glotzer, “Self-Assembly of Patchy Particles into Diamond Structures through Molecular Mimicry,” Langmuir 21, 11547–11551 (2005).
  * Glotzer and Solomon (2007) S. C. Glotzer and M. J. Solomon, “Anisotropy of building blocks and their assembly into complex structures,” Nature Materials 6, 557–562 (2007).
  * Cersonsky  _et al._ (2018a) R. K. Cersonsky, J. Dshemuchadse, J. Antonaglia, G. Van Anders, and S. C. Glotzer, “Pressure-tunable photonic band gaps in an entropic colloidal crystal,” Physical Review Materials 2, 125201 (2018a).
  * Cersonsky  _et al._ (2018b) R. K. Cersonsky, G. Van Anders, P. M. Dodd, and S. C. Glotzer, “Relevance of packing to colloidal self-assembly,” Proceedings of the National Academy of Sciences 115, 1439–1444 (2018b).
  * Keys, Iacovella, and Glotzer (2011) A. S. Keys, C. R. Iacovella, and S. C. Glotzer, “Characterizing complex particle morphologies through shape matching: Descriptors, applications, and algorithms,” Journal of Computational Physics 230, 6438–6463 (2011).
  * Van Anders  _et al._ (2014) G. Van Anders, D. Klotsa, N. K. Ahmed, M. Engel, and S. C. Glotzer, “Understanding shape entropy through local dense packing,” Proceedings of the National Academy of Sciences 111 (2014), 10.1073/pnas.1418159111.
  * Hamley (2007) I. W. Hamley,  _Introduction to Soft Matter– Revised Edition_ (John Wiley & Sons, Ltd, Chichester, UK, 2007).
  * Ramasubramani  _et al._ (2020) V. Ramasubramani, B. D. Dice, E. S. Harper, M. P. Spellings, J. A. Anderson, and S. C. Glotzer, “freud: A software suite for high throughput analysis of particle simulation data,” Computer Physics Communications 254, 107275 (2020).
  * Steinhardt, Nelson, and Ronchetti (1983) P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, “Bond-orientational order in liquids and glasses,” Physical Review B 28, 784–805 (1983), publisher: American Physical Society.
  * Lechner and Dellago (2008) W. Lechner and C. Dellago, “Accurate determination of crystal structures based on averaged local bond order parameters,” The Journal of Chemical Physics 129, 114707 (2008).
  * Prodan and Kohn (2005) E. Prodan and W. Kohn, “Nearsightedness of electronic matter,” Proceedings of the National Academy of Sciences 102, 11635–11638 (2005), publisher: Proceedings of the National Academy of Sciences.
  * Goscinski  _et al._ (2021a) A. Goscinski, G. Fraux, G. Imbalzano, and M. Ceriotti, “The role of feature space in atomistic learning,” Machine Learning: Science and Technology 2, 025028 (2021a), publisher: IOP Publishing.
  * Goscinski  _et al._ (2023) A. Goscinski, V. P. Principe, G. Fraux, S. Kliavinek, B. A. Helfrecht, P. Loche, M. Ceriotti, and R. K. Cersonsky, “scikit-matter : A Suite of Generalisable Machine Learning Methods Born out of Chemistry and Materials Science,” Open Research Europe 3, 81 (2023).
  * Everaers and Ejtehadi (2003) R. Everaers and M. R. Ejtehadi, “Interaction potentials for soft and hard ellipsoids,” Physical Review E 67, 041710 (2003), arXiv:cond-mat/0306096.
  * Berardi, Fava, and Zannoni (1995) R. Berardi, C. Fava, and C. Zannoni, “A generalized Gay-Berne intermolecular potential for biaxial particles,” Chemical Physics Letters 236, 462–468 (1995).
  * Perram and Wertheim (1985) J. W. Perram and M. S. Wertheim, “Statistical mechanics of hard ellipsoids. I. Overlap algorithm and the contact function,” Journal of Computational Physics 58, 409–416 (1985).
  * Note (1) The $h_{12}$ distances here fall within the confidence region of the Perram and Wertheim (1985) algorithm, so they can be taken as exact.
  * Stukowski (2009) A. Stukowski, “Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool,” Modelling and Simulation in Materials Science and Engineering 18, 015012 (2009).
  * Fredericks  _et al._ (2021) S. Fredericks, K. Parrish, D. Sayre, and Q. Zhu, “PyXtal: A Python library for crystal structure generation and symmetry analysis,” Computer Physics Communications 261, 107810 (2021).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, “QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials,” Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Prandini  _et al._ (2018) G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, and N. Marzari, “Precision and efficiency in solid-state pseudopotential calculations,” npj Computational Materials 4, 1–13 (2018), number: 1 Publisher: Nature Publishing Group.
  * Grimme  _et al._ (2016) S. Grimme, A. Hansen, J. G. Brandenburg, and C. Bannwarth, “Dispersion-Corrected Mean-Field Electronic Structure Methods,” Chemical Reviews 116, 5105–5154 (2016), publisher: American Chemical Society.
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, “Special points for Brillouin-zone integrations,” Physical Review B 13, 5188–5192 (1976), publisher: American Physical Society.
  * Ramasubramani  _et al._ (2018) V. Ramasubramani, C. Adorf, P. Dodd, B. Dice, and S. Glotzer, “signac: A Python framework for data and workflow management,” (Austin, Texas, 2018) pp. 152–159.
  * Adorf  _et al._ (2018) C. S. Adorf, P. M. Dodd, V. Ramasubramani, and S. C. Glotzer, “Simple data and workflow management with the signac framework,” Computational Materials Science 146, 220–229 (2018).
  * Hjorth Larsen  _et al._ (2017) A. Hjorth Larsen, J. Jørgen Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. Bjerre Jensen, J. Kermode, J. R. Kitchin, E. Leonhard Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. Bergmann Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, and K. W. Jacobsen, “The atomic simulation environment—a Python library for working with atoms,” Journal of Physics: Condensed Matter 29, 273002 (2017).
  * Cheng  _et al._ (2020) B. Cheng, R.-R. Griffiths, S. Wengert, C. Kunkel, T. Stenczel, B. Zhu, V. L. Deringer, N. Bernstein, J. T. Margraf, K. Reuter, and G. Csanyi, “Mapping Materials and Molecules,” Accounts of Chemical Research 53, 1981–1991 (2020), publisher: American Chemical Society.
  * Fraux  _et al._ (2023) G. Fraux, P. Loche, S. Kliavinek, F. Bigi, K. Hugeunin-Dumittan, A. Grisafi, R. K. Cersonsky, M. Kellner, D. Tisi, S. Shah, A. Goscinski, and M. Ceriotti, “rascaline,”  (2023).
  * Pozdnyakov  _et al._ (2020) S. N. Pozdnyakov, M. J. Willatt, A. P. Bartók, C. Ortner, G. Csányi, and M. Ceriotti, “On the Completeness of Atomic Structure Representations,” Physical Review Letters 125, 166001 (2020), arXiv:2001.11696 [cond-mat, physics:physics].
  * Note (2) Principal covariates regressionHelfrecht  _et al._ (2020), while ideal for regression-related dimensionality reductions, would again limit our ability to compare the _relevant information content_ of these representations.
  * Note (3) Using the same major semiaxis length as the width of the isotropic Gaussian field, and keeping all other hyperparameters consistent.
  * Helfrecht  _et al._ (2020) B. A. Helfrecht, R. K. Cersonsky, G. Fraux, and M. Ceriotti, “Structure-property maps with Kernel principal covariates regression,” Machine Learning: Science and Technology 1, 045021 (2020), publisher: IOP Publishing.
  * Fraux, Cersonsky, and Ceriotti (2020) G. Fraux, R. K. Cersonsky, and M. Ceriotti, “Chemiscope: interactive structure-property explorer for materials and molecules,” Journal of Open Source Software 5, 2117 (2020).
  * Talirz  _et al._ (2020) L. Talirz, S. Kumbhar, E. Passaro, A. V. Yakutovich, V. Granata, F. Gargiulo, M. Borelli, M. Uhrin, S. P. Huber, S. Zoupanos, C. S. Adorf, C. W. Andersen, O. Schütt, C. A. Pignedoli, D. Passerone, J. VandeVondele, T. C. Schulthess, B. Smit, G. Pizzi, and N. Marzari, “Materials Cloud, a platform for open computational science,” Scientific Data 7, 299 (2020), number: 1 Publisher: Nature Publishing Group.
  * (61) “DLMF: 14.10 Recurrence Relations and Derivatives: Real Arguments: Chapter 14 Legendre and Related Functions,” .
  * Goscinski  _et al._ (2021b) A. Goscinski, F. Musil, S. Pozdnyakov, J. Nigam, and M. Ceriotti, “Optimal radial basis for density-based atomic representations,” The Journal of Chemical Physics 155, 104106 (2021b).
  * Corbato (1969) F. Corbato, “A Paging Experiment with the Multics System,” in  _MIT Press_ (1969).
  * Matsakis and Klock (2014) N. D. Matsakis and F. S. Klock, “The rust language,” ACM SIGAda Ada Letters 34, 103–104 (2014).

##  Appendix A Evaluating the Coefficients

We have shown that, by choosing a MVG density and suitable basis functions, we
can analytically compute the coefficients

$\displaystyle\langle nlm|\rho_{i}\rangle,$ (27)

by evaluating of integrals of the form

$\displaystyle
I_{nlm}=\int_{\mathbb{R}^{3}}\mathrm{d}^{3}re^{-\frac{1}{2}(\bm{r}-\bm{a})^{T}A(\bm{r}-\bm{a})}r^{2n}R_{l}^{m}(\bm{r}).$
(28)

In this section, we will present an explicit algorithm to obtain all
$I_{nlm}$. The algorithm will require two key ingredients, namely converting
the spherical formulation of $r^{2n}R_{l}^{m}(\bm{r})$ into Cartesian
coordinates (A.2), then computing the corresponding moments across the subject
to our MVG (A.3).

###  A.1 Computing the effective quadratic form in the exponential

If we use the monomial basis, the integrand in

$\displaystyle I_{nlm}$
$\displaystyle=\int_{\mathbb{R}^{3}}\mathrm{d}^{3}re^{-\frac{1}{2}(\bm{r}-\bm{r}_{ij})^{T}A(\bm{r}-\bm{r}_{ij})}r^{2n}R_{l}^{m}(\bm{r})$
(29)

is already in a convenient form, since the Gaussian part is completely
specified by the center $\bm{r}_{ij}$ and the precision matrix $A$. If, on the
other hand, we use the GTO basis, we get an extra exponential factor

$\displaystyle I_{nlm}$
$\displaystyle=\int_{\mathbb{R}^{3}}\mathrm{d}^{3}re^{-\frac{1}{2}(\bm{r}-\bm{r}_{ij})^{T}A(\bm{r}-\bm{r}_{ij})}r^{2n}R_{l}^{m}(\bm{r})e^{-\frac{\bm{r}^{2}}{2\sigma^{2}}}$
(30)
$\displaystyle=\int_{\mathbb{R}^{3}}\mathrm{d}^{3}re^{-\frac{1}{2}(\bm{r}-\bm{r}_{ij})^{T}A(\bm{r}-\bm{r}_{ij})-\frac{\bm{r}^{2}}{2\sigma^{2}}}r^{2n}R_{l}^{m}(\bm{r})$
(31)

meaning that the Gaussian part

$\displaystyle
e^{-\frac{1}{2}(\bm{r}-\bm{r}_{ij})^{T}A(\bm{r}-\bm{r}_{ij})-\frac{\bm{r}^{2}}{2\sigma^{2}}}$
(32)

as a whole no longer has the convenient form. By completing the square,
ignoring the global factor of $-\frac{1}{2}$ in the exponent, we can rewrite

$\displaystyle(\bm{r}-\bm{r}_{ij})^{T}A(\bm{r}-\bm{r}_{ij})+\frac{1}{\sigma^{2}}\bm{r}^{2}$
$\displaystyle=(\bm{r}-\bm{r}_{0})^{T}\tilde{A}(\bm{r}-\bm{r}_{0})+c,$ (33)

with

$\displaystyle\tilde{A}$ $\displaystyle=A+\frac{1}{\sigma^{2}}$ (34)
$\displaystyle\bm{r}_{0}$
$\displaystyle=\tilde{A}^{-1}A\bm{r}_{ij}=\bm{r}_{ij}-\frac{1}{\sigma^{2}}\tilde{A}^{-1}\bm{r}_{ij}$
(35) $\displaystyle c$
$\displaystyle=\frac{1}{\sigma^{2}}\bm{r}_{ij}^{T}\tilde{A}^{-1}A\bm{r}_{ij}$
(36)

The second form of $\bm{r}_{0}$ is more convenient to obtain a qualitative
picture: the second term represents the deviation of the center due to the
addition of the second Gaussian. Depending on how sharp this Gaussian is, the
relative importance of this term will change.

Firstly, note that using the first representation of $\bm{r}_{0}$ and the fact
that both $A$ and $\tilde{A}$ are symmetric, we get

$\displaystyle\bm{r}_{0}^{T}\tilde{A}\bm{r}_{0}$
$\displaystyle=\bm{r}_{ij}^{T}A\tilde{A}^{-1}\tilde{A}\tilde{A}^{-1}A\bm{r}_{ij}$
(37) $\displaystyle=\bm{r}_{ij}^{T}A\tilde{A}^{-1}A\bm{r}_{ij}$ (38)
$\displaystyle=\bm{r}_{ij}^{T}A\tilde{A}^{-1}\left(\tilde{A}-\frac{1}{\sigma^{2}}\right)\bm{r}_{ij}$
(39)
$\displaystyle=\bm{r}_{ij}^{T}A\bm{r}_{ij}-\frac{1}{\sigma^{2}}\bm{r}_{ij}^{T}\tilde{A}^{-1}A\bm{r}_{ij}$
(40)

Expanding the quadratic terms, we get

$\displaystyle(\bm{r}-\bm{r}_{0})^{T}A(\bm{r}-\bm{r}_{0})$
$\displaystyle=\bm{r}^{T}\tilde{A}\bm{r}-2\bm{r}^{T}\tilde{A}\bm{r}_{0}+\bm{r}_{0}^{T}\tilde{A}\bm{r}_{0}$
(41)
$\displaystyle=\bm{r}^{T}A\bm{r}+\frac{1}{\sigma^{2}}\bm{r}^{2}-2\bm{r}^{T}\tilde{A}\left(\bm{r}_{ij}-\frac{1}{\sigma^{2}}\tilde{A}^{-1}\bm{r}_{ij}\right)+\bm{r}_{0}^{T}\tilde{A}\bm{r}_{0}$
(42)
$\displaystyle=\bm{r}^{T}A\bm{r}+\frac{1}{\sigma^{2}}\bm{r}^{2}-2\bm{r}^{T}\tilde{A}\bm{r}_{ij}+2\frac{1}{\sigma^{2}}\bm{r}^{T}\bm{r}_{ij}+\bm{r}_{0}^{T}\tilde{A}\bm{r}_{0}$
(43)
$\displaystyle=\bm{r}^{T}A\bm{r}+\frac{1}{\sigma^{2}}\bm{r}^{2}-2\bm{r}^{T}A\bm{r}_{ij}+\bm{r}_{0}^{T}\tilde{A}\bm{r}_{0}$
(44)
$\displaystyle=\bm{r}^{T}A\bm{r}+\frac{1}{\sigma^{2}}\bm{r}^{2}-2\bm{r}^{T}A\bm{r}_{ij}+\bm{r}_{ij}^{T}A\bm{r}_{ij}-\frac{1}{\sigma^{2}}\bm{r}_{ij}^{T}\tilde{A}^{-1}A\bm{r}_{ij}$
(45)
$\displaystyle=(\bm{r}-\bm{r}_{ij})^{T}A(\bm{r}-\bm{r}_{ij})+\frac{1}{\sigma^{2}}\bm{r}^{2}-\frac{1}{\sigma^{2}}\bm{r}_{ij}^{T}\tilde{A}^{-1}A\bm{r}_{ij}$
(46)

Thus, we do see that indeed,

$\displaystyle(\bm{r}-\bm{r}_{ij})^{T}A(\bm{r}-\bm{r}_{ij})+\frac{1}{\sigma^{2}}\bm{r}^{2}$
$\displaystyle=(\bm{r}-\bm{r}_{0})^{T}\tilde{A}(\bm{r}-\bm{r}_{0})+\frac{1}{\sigma^{2}}\bm{r}_{ij}^{T}\tilde{A}^{-1}A\bm{r}_{ij}.$
(47)

###  A.2 Spherical to Cartesian Transformation

Recall that for both the monomial and GTO bases, our integrand is a product of
a Gaussian and an expression of the form

$\displaystyle r^{2n}R_{l}^{m}(\bm{r})=\mathrm{poly}(x,y,z),$ (48)

where we emphasize that $r^{2n}R_{l}^{m}(\bm{r})$ is simply a polynomial in
the variables $(x,y,z)$.
For practical evaluations, however, it does not suffice to know that it is
equal to “some” polynomial. We will need to explicitly express the solid
harmonics $R_{l}^{m}(\bm{r})$ in monomial terms. In general, there will exist
a decomposition

$\displaystyle
R_{l}^{m}(\bm{r})=\sum_{n_{0}+n_{1}+n_{2}=l}T^{lm}_{n_{0},n_{1}n_{2}}x^{n_{0}}y^{n_{1}}z^{n_{2}},$
(49)

where $T^{lm}_{n_{0},n_{1}n_{2}}$ are some coefficients. These coefficients
depend only on $l,m$ and can thus be cached. For the full basis function,
there will then exist coefficients $T^{nlm}_{n_{0},n_{1}n_{2}}$ such that

$\displaystyle
r^{2n}R_{l}^{m}(\bm{r})=\sum_{n_{0},n_{1},n_{2}}T^{nlm}_{n_{0},n_{1}n_{2}}x^{n_{0}}y^{n_{1}}z^{n_{2}}.$
(50)

As soon as we are provided with a complete list of the $(l,n)$ pair we need,
we can precompute these coefficients.

###  A.3 Evaluation of Moments

Once we have decomposed $r^{2n}R_{l}^{m}(\bm{r})$ into a sum of monomial
terms, we have reduced the evaluation problem to integrals of the form

$\displaystyle\langle
x^{n_{0}}y^{n_{1}}z^{n_{2}}\rangle=\int_{\mathbb{R}^{3}}\mathrm{d}^{3}re^{-\frac{1}{2}(\bm{r}-\bm{a})^{T}A(\bm{r}-\bm{a})}x^{n_{0}}y^{n_{1}}z^{n_{2}}.$
(51)

We shall call this expression the $(n_{0},n_{1},n_{2})$-th moment of the
Gaussian, or just moment for short.

To provide some context for the terminology, “moment” is a term used in
probability theory. If we are given some probability density (in three
variables) $p(\bm{r})=p(x,y,z)$, its $(n_{0},n_{1},n_{2})$-th moment is
defined as

$\displaystyle\int_{\mathbb{R}^{3}}\mathrm{d}^{3}rp(\bm{r})x^{n_{0}}y^{n_{1}}z^{n_{2}}.$
(52)

The most common notation for this in the mathematical literature is to write
it as $\mathbb{E}[x^{n_{0}}y^{n_{1}}z^{n_{2}}]$. In the quantum field theory /
statistical mechanics literature, on the other hand, the notation $\langle
x^{n_{0}}y^{n_{1}}z^{n_{2}}\rangle$ is used more commonly instead.

Connecting this back to our integrals, we can see that the Gaussian part is
almost a probability density, apart from the normalization factor, motivating
the use of an analogous notation.

The normalization factor is a global factor, meaning that techniques that have
been developed to evaluate moments in probability theory can be applied to our
problem as well. This will allow us to compute all the moments for a given
Gaussian function specified by the precision matrix $A$ and the center
$\bm{a}$.

###  A.4 Putting Everything Together

Combining the two ingredients, we can now formulate the complete algorithm for
the evaluation of the above mentioned integrals, namely:

$\displaystyle I_{nlm}$
$\displaystyle=\int_{\mathbb{R}^{3}}\mathrm{d}^{3}re^{-\frac{1}{2}(\bm{r}-\bm{a})^{T}A(\bm{r}-\bm{a})}r^{2n}R_{l}^{m}(\bm{r})$
(53)
$\displaystyle=\int_{\mathbb{R}^{3}}\mathrm{d}^{3}re^{-\frac{1}{2}(\bm{r}-\bm{a})^{T}A(\bm{r}-\bm{a})}\sum_{n_{0},n_{1},n_{2}}T^{n,lm}_{n_{0},n_{1}n_{2}}x^{n_{0}}y^{n_{1}}z^{n_{2}}$
(54)
$\displaystyle=\sum_{n_{0},n_{1},n_{2}}T^{n,lm}_{n_{0},n_{1}n_{2}}\int_{\mathbb{R}^{3}}\mathrm{d}^{3}re^{-\frac{1}{2}(\bm{r}-\bm{a})^{T}A(\bm{r}-\bm{a})}x^{n_{0}}y^{n_{1}}z^{n_{2}}$
(55) $\displaystyle=\sum_{n_{0},n_{1},n_{2}}T^{n,lm}_{n_{0},n_{1}n_{2}}\langle
x^{n_{0}}y^{n_{1}}z^{n_{2}}\rangle$ (56)

We can therefore see that the computation of our features is completely
determined by the two ingredients mentioned before, namely:

  1. 1.

the transformation coefficients $T^{n,lm}_{n_{0},n_{1}n_{2}}$ that express
$r^{2n}R_{l}^{m}$ using monomials

  2. 2.

the moments $\langle x^{n_{0}}y^{n_{1}}z^{n_{2}}\rangle$ of the Gaussian

In the following two subsections, we will explain how to perform these two
steps, respectively.

###  A.5 Evaluation of Solid Harmonics, Spherical Harmonics Part

In this subsection, we present an algorithm that will allow us to decompose
the basis functions into its constituent monomial terms. For this, we start
from the definition

$\displaystyle
R_{l}^{m}(\bm{r})=\sqrt{\frac{4\pi}{2l+1}}r^{l}Y_{l}^{m}(\hat{r})=\sqrt{\frac{(l-m)!}{(l+m)!}}r^{l}P_{l}^{m}(\cos\theta)e^{im\phi}.$
(57)

The $\phi$ dependent part will be a polynomial of degree $m$, while the
associated Legendre polynomial will provide us with the remaining $l-m$
degrees. Starting with the exponential, we get

$\displaystyle e^{im\phi}$ $\displaystyle=\left(e^{i\phi}\right)^{m}$ (58)
$\displaystyle=\left(\cos\phi+i\sin\phi\right)^{m}$ (59)
$\displaystyle=\sum_{k=0}^{m}\begin{pmatrix}m\\\
k\end{pmatrix}i^{m-k}\cos^{k}\phi\sin^{m-k}\phi$ (60)

After multiplying by $r^{m}\sin^{m}\theta$, this will result in contributions
of the form

$\displaystyle\sum_{k=0}^{m}\begin{pmatrix}m\\\
k\end{pmatrix}i^{m-k}x^{k}y^{m-k}.$ (61)

For the associated Legendre polynomial part, it is possible to use the
recurrence relationnoa

$\displaystyle
P_{l+1}^{m}(x)=\frac{2l+1}{l-m+1}xP_{l}^{m}-\frac{l+m}{l-m+1}P_{l-1}^{m}$ (62)

and

$\displaystyle P_{l}^{l}(x)=(-1)^{l}(2l-1)!!(1-x^{2})^{\frac{l}{2}}$ (63)

The special case for $l=m$ can be used to initialize the recurrence relations
by computing all the associated Legendre polynomials at $l=m$ and $l=-m$ with
symmetry, while the recurrence relation can be used to generate all higher
order ones.

Combining Eqs. 57, 62 and 63, we get

$\displaystyle R_{l+1}^{m}$
$\displaystyle=\sqrt{\frac{(l+1-m)!}{(l+1+m)!}}\sqrt{\frac{(l+m)!}{(l-m)!}}\frac{2l+1}{l+1-m}zR_{l}^{m}-$
(64)
$\displaystyle\sqrt{\frac{(l+1-m)!}{(l+1+m)!}}\sqrt{\frac{(l-1+m)!}{(l-1-m)!}}\frac{l+m}{l+1-m}r^{2}R_{l-1}^{m}$
(65)

from where we can identify the prefactors used in the algorithm below The
equation be further simplified to obtain the recurrence as,

$\displaystyle R_{l+1}^{m}$
$\displaystyle=\sqrt{\frac{(l+1-m)}{(l+1+m)}}\frac{2l+1}{l+1-m}zR_{l}^{m}-$
(66)
$\displaystyle\sqrt{\frac{(l+1-m)!}{(l+1+m)!}}\sqrt{\frac{(l-1+m)!}{(l-1-m)!}}\frac{l+m}{l+1-m}r^{2}R_{l-1}^{m}$
(67)

###  A.6 Evaluation of Radial Dependence

If all transformation coefficients for $n=0$ are known, it is straight forward
to compute those for higher $n\geq 2$ using

$\displaystyle
B_{n+1,lm}=r^{2}B_{nlm}=(x^{2}+y^{2}+z^{2})B_{nlm}\equiv(x^{2}+y^{2}+z^{2})r^{2n}R_{l}^{m}$
(68)

Thus, if all transformation coefficients $T^{nlm}_{n_{0},n_{1},n_{2}}$ at some
$n$ are known, we obtain those for higher $n$ by running the iteration:

$\displaystyle\forall(l,m,n_{0},n_{1},n_{2}):$
$\displaystyle\mathrm{with}\quad n_{0}+n_{1}+n_{2}=l+2n=d$ (69) $\displaystyle
T^{n+1,lm}_{n_{0}+2,n_{1},n_{2}}+=T^{n,lm}_{n_{0},n_{1}n_{2}}$ (70)
$\displaystyle T^{n+1,lm}_{n_{0},n_{1}+2,n_{2}}+=T^{n,lm}_{n_{0},n_{1}n_{2}}$
(71) $\displaystyle
T^{n+1,lm}_{n_{0},n_{1},n_{2}+2}+=T^{n,lm}_{n_{0},n_{1}n_{2}}$ (72)

After the iteration, all coefficients at the radial channel $n+1$ will have
the correct values.

###  A.7 Computing the Moments

In this subsection, we explain how we can compute the moments

$\displaystyle\langle
x^{n_{0}}y^{n_{1}}z^{n_{2}}\rangle=\int_{\mathbb{R}^{3}}\mathrm{d}^{3}re^{-\frac{1}{2}(\bm{r}-\bm{a})^{T}\mathbf{A}(\bm{r}-\bm{a})}x^{n_{0}}y^{n_{1}}z^{n_{2}}.$
(73)

for a given precision matrix $A$ and center $\bm{a}$.
As we will show, this step will be significantly easier if the matrix
$\mathbf{A}$ is diagonal, i.e. if we work in the basis of the principal axes.
In our current implementation, we still use the original coordinate frame in
which the matrix $\mathbf{A}$ keeps its general form. Despite the slightly
more complicated computation of the moments, this will reduce some
computational cost associated with the rotation of the coordinate frame
(summing over the Wigner matrix) later on. Nevertheless, there might still be
some benefits to have an implementation that uses the diagonalization step. We
will, therefore, start by discussing the simpler diagonal case, as well as
extensions that would become possible for this special case. We then move on
to the algorithm for general $A$, which is what we are using in the current
implementation.

##### Diagonal Case

In the diagonal case, the integral separates into an x-, y- and z-dependent
part. The total integral can thus easily be evaluated by using the explicit
formula for the one-dimensional Gaussian integral, namely:

$\displaystyle\langle x^{n_{0}}y^{n_{1}}z^{n_{2}}\rangle$
$\displaystyle=\int_{\mathbb{R}^{3}}\mathrm{d}^{3}rx^{n_{1}}y^{n_{2}}z^{n_{3}}e^{-\frac{1}{2}\left(\frac{x^{2}}{\sigma_{1}^{2}}+\frac{y^{2}}{\sigma_{2}^{2}}+\frac{z^{2}}{\sigma_{3}^{2}}\right)}$
(74)
$\displaystyle=\left(2\sigma^{2}_{1}\right)^{n_{1}^{\prime}}\Gamma(n_{1}^{\prime})\left(2\sigma^{2}_{2}\right)^{n_{2}^{\prime}}\Gamma(n_{2}^{\prime})\left(2\sigma^{2}_{3}\right)^{n_{3}^{\prime}}\Gamma(n_{3}^{\prime}).$
(75)

An advantage of the diagonal formalism is that it can easily be adapted to
non-Gaussian densities. We could, for instance, consider densities of the form

$\displaystyle
g_{0}(\bm{r})=\exp\left[-\left(\frac{|x|}{\sigma_{x}}+\frac{|y|}{\sigma_{y}}+\frac{|z|}{\sigma_{z}}\right)\right]$
(76)

This could be used to have densities that have the symmetry of a cube (if all
$\sigma_{j}$ are equal) or other rectangular shapes. While integrating these
in spherical coordinates would be a nightmare due to the absolute values, in
Cartesian coordinates, the computation is relatively simple once the problem
is brought into this diagonal form. There could, therefore, be some value in
having an implementation that works for the diagonalized case. It should be
noted, however, that we are not directly working with the density. Instead, we
are approximating the density with our basis functions. If the resolution of
our basis functions is low, it might not be possible to properly distinguish
Gaussians from rectangular densities. This should be kept in mind before
spending too much time developing complicated schemes with minimal additional
improvements.

##### General Case

Here, we use an iterative updating scheme. For later convenience, we define
the covariance matrix $C=A^{-1}$ and work with it instead. To keep the
notation closer to the implementation, we will index its components as $0,1,2$
rather than $1,2,3$, corresponding to the coordinate axes $x,y,z$.
We can initialize the first few moments, up to a global factor, by

$\displaystyle\langle 1\rangle$ $\displaystyle=\langle
x^{0}y^{0}z^{0}\rangle=1$ (77) $\displaystyle\langle x\rangle$
$\displaystyle=a_{x}$ (78) $\displaystyle\langle y\rangle$
$\displaystyle=a_{y}$ (79) $\displaystyle\langle z\rangle$
$\displaystyle=a_{z}$ (80)

We can compute all higher-order moments by using three recurrence
relationsNigam, Pozdnyakov, and Ceriotti (2020), one in the x-, y-, and
z-directions, respectively. These are given by:
x-iteration:

$\displaystyle\langle x^{n_{0}+1}y^{n_{1}}z^{n_{2}}\rangle$
$\displaystyle=a_{0}\langle
x^{n_{0}}y^{n_{1}}z^{n_{2}}\rangle+C_{00}n_{0}\langle
x^{n_{0}-1}y^{n_{1}}z^{n_{2}}\rangle$ (81) $\displaystyle+C_{01}n_{1}\langle
x^{n_{0}}y^{n_{1}-1}z^{n_{2}}\rangle+C_{02}n_{2}\langle
x^{n_{0}}y^{n_{1}}z^{n_{2}-1}\rangle$ (82)

y-iteration:

$\displaystyle\langle x^{n_{0}}y^{n_{1}+1}z^{n_{2}}\rangle$
$\displaystyle=a_{1}\langle
x^{n_{0}}y^{n_{1}}z^{n_{2}}\rangle+C_{10}n_{0}\langle
x^{n_{0}-1}y^{n_{1}}z^{n_{2}}\rangle$ (83) $\displaystyle+C_{11}n_{1}\langle
x^{n_{0}}y^{n_{1}-1}z^{n_{2}}\rangle+C_{12}n_{2}\langle
x^{n_{0}}y^{n_{1}}z^{n_{2}-1}\rangle$ (84)

z-iteration:

$\displaystyle\langle x^{n_{0}}y^{n_{1}}z^{n_{2}+1}\rangle$
$\displaystyle=a_{2}\langle
x^{n_{0}}y^{n_{1}}z^{n_{2}}\rangle+C_{20}n_{0}\langle
x^{n_{0}-1}y^{n_{1}}z^{n_{2}}\rangle$ (85) $\displaystyle+C_{21}n_{1}\langle
x^{n_{0}}y^{n_{1}-1}z^{n_{2}}\rangle+C_{22}n_{2}\langle
x^{n_{0}}y^{n_{1}}z^{n_{2}-1}\rangle.$ (86)

Warning: Please note that for any of the above iterations, it is possible that
some exponents will become negative. A concrete example, in which we naively
use the above formulae, would be the evaluation of $\langle x^{2}\rangle$ from
the $x$-iteration, which would give

$\displaystyle\langle x^{2}\rangle$ $\displaystyle=a_{0}\langle
x^{1}\rangle+C_{00}\cdot 1\cdot\langle 1\rangle+C_{01}\cdot 0\cdot\langle
xy^{-1}\rangle+C_{02}\cdot 0\cdot\langle xz^{-1}\rangle$ (87)

In such cases, the terms containing negative exponents, strictly speaking, do
not exist. As the multiplication by zero suggests, however, these terms do not
contribute to the final result. The correct iteration in this case therefore
is

$\displaystyle\langle x^{2}\rangle$ $\displaystyle=a_{0}\langle
x^{1}\rangle+C_{00}n_{0}\langle 1\rangle.$ (88)

In practice, depending on the specific implementation, it might not even be
necessary to make special cases for this if all array elements are initialized
to zero. Please also note that the initializations only apply to normalized
Gaussians that have a proper interpretation as a probability density. For
other normalizations, the coefficients will all be multiplied by the same
global factor. Thus, in practice, it makes sense to include this global factor
in the initialization, since the iterative scheme will then automatically
ensure that all higher order moments have the correct prefactor as well.

###  A.8 Orthonormalization

The coefficients $\langle nlm|\rho_{i}\rangle$ from the above procedures
result from using non-orthonormal bases. While not strictly necessary, working
in orthonormal bases ensures minimal overlap in mutual information between
featuresGoscinski  _et al._ (2021b). This discussion only applies to basis
functions that are square integrable (e.g. GTO); non-square-integrable basis
functions (e.g. monomial basis) cannot be orthonormalized. Below, we discuss
how to scale our coefficients, taking the GTO basis as a representative
example for any square-integrable basis.

We first choose the coefficients of our density expansion $\langle
nlm|\rho_{i}\rangle$ to use the unnormalized GTO basis
$R_{nl}(r)=r^{l+2n}e^{\frac{-r^{2}}{2\sigma^{2}}}$ for our radial expansion.
The general unnormalized GTO
$\phi_{d}=r^{d}e^{\frac{-r^{2}}{2\sigma_{d}^{2}}}$ has a finite square-
integral over $\mathbb{R}^{3}$:

$I_{d}=\int_{0}^{\infty}{|\phi_{d}(r)|^{2}*r^{2}dr}=2^{-1}\sigma_{d}^{2d+3}*\Gamma(\frac{2d+3}{2})$

and can hence be normalized with the constant $N_{d}=1/\sqrt{I_{d}}$. In other
words, our normalized GTO is $\Phi_{d}=N_{d}*\phi_{d}$, and
$\int_{0}^{\infty}{|\Phi_{d}(r)|^{2}*r^{2}dr}=1$. After normalizing all the
bases, we can take orthogonalize them using Löwdin Symmetric
OrthonormalizationLöwdin (1970). We first find the overlap matrix between two
normalized GTOs:

$G_{ij}=\int_{0}^{\infty}\Phi_{i}\Phi_{j}r^{2}dr.$

The orthonormalization matrix is the inverse square root of the overlap matrix
$G^{-1/2}$, which is guaranteed to exist because the overlap matrix is a gram
matrix and is hence symmetric positive definite. Specifically, $G^{-1/2}$ is
calculated by diagonalizing $G$, then taking the recipricol of the square root
of the diagonal matrix: $G=MDM^{T}$, $G^{-1/2}=MD^{-1/2}M^{T}$. Then, we can
apply this matrix to obtain a set of orthonormal GTO basis vectors
$\hat{\Phi}_{i}$: $\hat{\Phi}_{i}=G^{-1/2}_{ij}\Phi_{j}$. We note that in our
previous procedure, we could in theory construct an overlap matrix from
unnormalized GTOs, but practically, calculating the overlap between
unnormalized GTOs of high order, then inverting them to find the
orthonormalization matrix, is not numerically stable.

###  A.9 Tricks for Efficient Implementation

The current implementation of AniSOAP is performance limited, both by
inefficient recalculations of the Clebsch-Gordan Matrices, and by the large
number of computations and nested iterations performed in Python when
calculating high-order moments. Below we outline two strategies under active
development to address these inefficiencies. The computational costs of the
current implementation and the proposed future implementations are shown in
Fig. S1.

Figure S1: Time measurement of major bottlenecks with and without caching. The
same computation was performed for all N by N by M crystal structures with N
and M both ranging from 1 to 16, and with fixed
$l_{max}=6,\sigma=5.0,r_{cut}=1.0$. With caching, after the first computation,
all subsequent computations of Clebsch-Gordan matrices are near instantaneous
($\sim 10\mu s$, the time required to access a list). The gray dashed line
show the expected $\sim$80x speedup when we port the moments calculation to a
lower-level language like Rust. The red and gray lines approach a slope of 1
on the log-log plot, indicating linear scaling with system size.

####  A.9.1 Caching the Clebsch-Gordan (CG) Matrix

The construction of Clebsch-Gordan matrices depends on only the hyperparameter
$l_{max}$. Originally, the matrices were stored within a Python dictionary
indexed by $(l_{1},l_{2},L)$. Since $0\leq l_{1},l_{2}\leq l_{max}$ and
$|l_{1}-l_{2}|\leq L\leq\min(l_{max},(l_{1}+l_{2}))$, the number of matrices
required to compute grows as $\mathcal{O}(l_{max}^{3})$, which is expensive
even for modest values of $l_{max}$.

However, since the construction of CG matrices only depends on $l_{max}$, we
can cache the matrices and re-use them internally. For caching, we use a
simple cyclic list that stores (key, value) pair with key storing $l_{max}$
and value storing corresponding matrices, as shown below, where $M$
corresponds to the appropriate CG matrices.

$key_{1}$$M_{1}$$key_{2}$$M_{2}$$\cdots$$key_{n-1}$$M_{n-1}$$key_{n}$$M_{n}$

The cache has a finite and fixed number of (key, value) pairs it can store to
prevent memory overflow, and whenever the cache is full, it decides the entry
for replacement using an algorithm that mimics the CLOCK algorithm for page
replacementCorbato (1969). In accordance with the algorithm, the key for the
CG list also contains a replacement bit to store the time-indexed usage of the
entry. The format of the CLOCK algorithm does place implicit limits on
$l_{max}$, such that $l_{max}\leq 2^{31}-1$, although this value is beyond
practical usage. The utility of caching these CG matrices is apparent when
performing repeated higher-body order expansions with the same $l_{max}$, such
as in a molecular dynamics simulation.

####  A.9.2 Moments Array

Given the simplicity of the moments calculations and the frequency of invoking
this section of code, we are currently porting this functionality to
RustMatsakis and Klock (2014), which gives the benefits of a compiled language
while retaining easy interfacing with Python.

While the Rust code is almost a one-to-one translation, there were two minor
changes. Firstly, the computation of the inverse of dilation matrix $\bm{D}$
(9) is changed. By definition, $\bm{D}$ is a 3 by 3 symmetric matrix.
Therefore, in Rust we compute the analytical formula for inverting this matrix
to avoid unnecessary operations while minimizing package dependency. Secondly,
the code was re-organized in Rust to maintain clarity in the new syntax.

####  A.9.3 Results of Optimization

With the CG-matrix caching and proposed porting of the moments above, we have
performed benchmark timings (Fig. S1) to compare the speed of execution. These
benchmarks were computed on the High-Performance Computers (HPC) within the
Wisconsin Center for High-Throughput Computation. Each HPC had a 3GHz AMD EPYC
7763 64-Core Processor and 514 GB Total RAM. Note that none of this code is
currently parallelized and therefore only a single core was used, yielding
timings that are similar to what is obtained on a local machine.

##  Appendix B Gay-Berne Case-Study

###  B.1 Gay-Berne Interactions

Continuing our discussion of Eq. (22), the second term, $\eta_{12}$, is a
function of each ellipsoids’ geometry and orientation but not a function of
distance, and is given as follows:

$\eta_{12}(\mathcal{R}_{1},\bm{S}_{1},\mathcal{R}_{2},\bm{S}_{2})=\left[\frac{2s_{1}s_{2}}{\det\bm{\Xi}_{12}(\mathcal{R}_{1},\mathcal{R}_{2})]}\right]^{\nu/2}$
(89) $s_{i}=[a_{i}b_{i}+c_{i}c_{i}][a_{i}b_{i}]^{1/2}$ (90)

This correction term describes the interaction strength between two ellipsoids
at 0 separation, whose influence is tuned by parameter $\nu$.

The last term, $\chi_{12}$, is only a function of each ellipsoid’s individual
orientations ($\mathcal{R}_{1},\mathcal{R}_{2}$) and relative orientation
($\hat{r_{12}}$), but not of their geometries ($\bm{S}_{1},\bm{S}_{2}$). It is
given as follows:

$\chi_{12}(\mathcal{R}_{1},\mathcal{R}_{2},\hat{r}_{12})=\left[\hat{r}_{12}^{T}\bm{B}_{12}^{-1}(\mathcal{R}_{1},\mathcal{R}_{2})\hat{r}_{12}\right]^{\mu}$
(91)
$\bm{B}_{12}(\mathcal{R}_{1},\mathcal{R}_{2})=\mathcal{R}_{1}^{T}\bm{E}_{1}\mathcal{R}_{1}+\mathcal{R}_{2}^{T}\bm{E}_{2}\mathcal{R}_{2}$
(92) $\bm{E}_{i}=\begin{pmatrix}&e_{ai}^{-1/\mu}&0&0\\\
&0&e_{bi}^{-1/\mu}&0\\\ &0&0&e_{ci}^{-1/\mu}\\\ \end{pmatrix}$ (93)

Note that for Gay-Berne, the rotations $\mathcal{R}_{i}$ transform from lab
frame to the body frame, while AniSOAP defines rotations to transform from the
body frame to the lab frame. Hence,
$\mathcal{R}_{i,{AniSOAP}~{}}\equiv\mathcal{R}_{i,GB}^{-1}$. The above
definitions use $\mathcal{R}_{i}\equiv\mathcal{R}_{i,GB}$.

$\chi_{12}$ corrects for the well-depth by interpolating the relative well
depth between the side-to-side, face-to-face, and end-to-end interactions,
given by $e_{ai}$, $e_{bi}$, $e_{ci}$. Generally, these pole-pole relative
well depths are arbitrarily fitted, but can be assumed to be equal to the
gaussian curvature of the ellipsoids at each pole, provided that $\mu=1$:

$\bm{E}_{i}=\sigma\begin{pmatrix}&\frac{a_{i}}{b_{i}c_{i}}&0&0\\\
&0&\frac{b_{i}}{a_{i}c_{i}}&0\\\ &0&0&\frac{c_{i}}{a_{i}b_{i}}\\\
\end{pmatrix}$ (94)

In our case-study with $a_{i}=1,b_{i}=1.5,c_{i}=2$ ellipsoids, we set $\mu=1$,
enabling the use of 94 to calculate $\bm{E}_{i}$. We furthermore set $\nu=1$,
completely specifying the Gay-Berne hyperparameters.

##  Appendix C Supplementary Information

###  C.1 Liquid Crystals

AniSOAP is well-suited for datasets containing multiple ellipsoid types;
however, it is not necessarily a fair comparison with the established order
parameters to demonstrate their inability to reconstruct such complex
datasets. Here we include additional analyses for identical frames as those
used in Sec. III.1, augmenting by another 1000 frames that represent $L/D=2$
spheroids.

Figure S2: Analogous PCA to Fig. 1, with both L/D=3 (left) and L/D=2 (right)
spheroids. Figure S3: Analogous to Fig. 2, the reconstruction of different
representations of our augmented dataset.

Our GFRE results emphasize that AniSOAP is able to contain analogous
information to both the Steinhardt and Orientational Order Parameters;
however, the reconstruction of AniSOAP decreases by either representation when
our dataset contains multiple particle types.

Figure S4: Learning curve from Gay-Berne dataset, showing RMSE values for
different training set sizes. The solid line is the average RMSE over 5 folds,
while the minimum and maximum RMSE denoted by the boundaries of the shaded
region. Note that $N_{train}=2000$ corresponds with a training size consistent
with 10% of our full dataset.

###  C.2 Gay-Berne Learning Curve

We computed the learning curves for the Gay-Berne dataset contained in the
main text using the learning_curve functionality in scikit-learn, assuming
cross-validated ridge regression as our underlying estimator, five-fold cross-
validation, and RMSE as the scoring mechanism. The learning curves for the
Gay-Berne dataset demonstrate that even small train-test splits can be used to
train a highly-accurate ridge-regression predictor on the AniSOAP
representation, showing its powerful interpolation abilities on both $h_{12}$
and energies.

###  C.3 Benzene Crystals Case-Study

####  C.3.1 Learning Curve

We computed learning curves for the dataset contained in the main text using
the scikit-learn learning curve functionality, assuming cross-validated ridge
regression as our underlying estimator, five-fold cross-validation, and RMSE
as the scoring mechanism. Our plot shows that we are within the saturation
regime with the learning of this dataset using the AniSOAP representation.
This learning curve also demonstrates that AniSOAP can effectively learn
benzene energetics within minimal training data.

Figure S5: Learning curve for benzene dataset, showing root-mean-squared-error
values for different training set sizes. Circles denote outliers in the
learning curve, as determined by those falling outside the second quartile of
the distribution by matplotlib.pyplot.boxplot.

####  C.3.2 Hyperparameter Tuning: Additional Results

Figure S6: Additional results from hyperparameter running for
$\sigma_{\text{GTO}}=0.5,1.0,1.5,2.0,2.5,3.0$Å. Color of the scatter and
contour lines denote $GFRE(X_{\text{AniSOAP}},X_{\text{SOAP}})$, where contour
levels are consistent across all plots and in the main text. For
$\sigma_{GTO}=0.5$, many points at $\sigma_{1,2}<3.0$ had results greater than
0.7; however, we truncated the color map to be consistent with the main text
and show the gradation of the oblate region more clearly.