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 coarsegrained 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 (variational 1,2 and hierarchical 3 auto-encoders, different forms of neural networks [4][5][6][7][8] , ensemble learning 9 ) 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 architectures [10][11][12] . 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 SMILES 13 or SELFIES 14 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 features 12,15 , including Behler-Parinello symmetry functions 16 , smooth overlap of atomic positions (SOAP) 17 , and molecular graphs 18 .
So, for coarse-grained systems, how do we apply a featureforward 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., SOAP 17 , NICE 19 ) 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 ex-tending 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 performancedefining 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.
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 formalism 20,21 .
Using the braket notation prescribed in Ref. 22, the contribution of an atom j to center atom i's expansion coefficients are given by
General Variables r r r a vector in Cartesian space r magnitude of vector r r r r direction of vector r r r r r r i j the vector between point i and point j R a 3x3 rotation matrix x x x a feature vector for one configuration X X X a matrix containing, as rows, feature vectors for multiple configurations Density Expansions g(r r r) a potentially anisotropic density function centered at r r r ρ TABLE I. Notation Guide. Throughout the text, we adopt the notation typical of the atom-centered symmetrized representation community as detailed in the table above.
where g(r r r) is a localized function (usually a Gaussian), R nl (r) is a radial basis function, and Y m l (r) are spherical harmonics. For a general g(r r 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 f (x)dx, which is why, after all, there are so many books on integral tables.
If we use any arbitrary density g(r r 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 Lebedev 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) basis 23 . The main downside of the numerical approach, at least for sufficiently general densities g(r r 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 r r r i j onto the zaxis of our coordinate system, eliminating any dependence of g(r r r -r r r i j ) on the angular integrands. This simplifies (1), even for complicated bases, to
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0210910
where F(r i j ) 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(r r r i j ) on the full relevant interval for efficient spline evaluations 22 .
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 σ 1 , σ 2 , σ 3 and corresponding axes. In order to maintain consistency with (1), we replace all instances of r r r with RR T r r r, where R is our rotation matrix into this component axis. We then simplify by changing our axes of integration from r r r to r r r ′ ≡ R T r r r and noting that r = RR T r to obtain
To perform the rotation of our spherical harmonics, one commonly sums over Wigner-D matrices 24 to obtain
With this step, we can reduce our calculation in the special case in which the matrix D l mm ′ (R) is diagonal, which proves useful for analytical evaluation techniques. We thus focus on the factor
which will be the focus of our discussion from hereon.
The simplest anisotropic function that can be put in place of g(r) in ( 5) is the multivariate Gaussian density (hereon MVG). The goal of this subsection is to show that by choosing g(r r 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
2 (r r r-r r r i j ) T D D D(r r r-r r r i j ) p(r r r),
where p(r r r) is some polynomial expression in r r r ≡ (x, y, z) and D D D is some 3 × 3 diagonal matrix. Consider the three-dimensional Gaussian defined by
where A is a symmetric and (strictly) positive definite 3 × 3 matrix. Any matrix satisfying these two conditions can be orthogonally diagonalized, also called the principal axis decomposition, allowing us to write
where R ∈ SO(3) is a rotation matrix that specifies the three principal axes and
is a diagonal matrix containing the widths σ j of the Gaussian along the three principal directions. With this decomposition, we can write
where r r r ′ = R T r r r are the coordinates with respect to the principal axes.
Putting the equation for our MVG into (5), we get
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, . . . , we have
where Γ is the gamma function. If we could reduce the evaluation of the integral (6) to such (multivariate) polynomials, this might provide us with analytical expressions for the expansion coefficients. Thus, following (13), we require that R nl (r)Y m l (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
Given that r l Y m l (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.
With one extra step, we can also show that a suitable modification of the GTO basis,
2σ 2 , provides an equally well-suited basis. . This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
By separating out all the exponential and polynomial factors, we can write the integrand as
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
with R m l (r r r) = r l Y m l (r) and where the diagonal matrix D D 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 Slatertype 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
where ⟨nlm|ρ i ⟩ = ∑ j ⟨nlm|ρ i j ⟩ . We go through the details on analytically computing eq. 17 in Appendix A 2-A 7. We first detail how to convert r 2n R m l into a cartesian polynomial (A 2, A 5), then detail how to evaluate the integrals (i.e. moments of the MVG) in A 3 and A 6. A 4 gives an overview of how these steps are combined together to result in evaluating a weighted sum of high-order moments of the MVG. Finally, we discuss imposing orthonormality within the basis set in Appendix A 7 via Löwdin symmetric orthogonalization 25 .
Even for general point clouds, determining the optimal anisotropic analog is non-trivial, and the focus of ongoing work 26 . 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 particles 27 , as many complex behaviors within molecular systems can be explained by analyzing the steric interactions of their molecular volumes [28][29][30][31][32] . Future studies will focus on optimizing AniSOAP for less trivial cases and further varieties of molecular anisotropy.
In Sections III A-III B, 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 C, 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. 18, but note the simplicity of continuing onto higher body-order terms for greater accuracy and transferability 19 .
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. 33 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 34 to classify the orientational alignment of particles, Steinhardt order parameters (SOPs 35 ) 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:
where j is a neighbor of i, and weights w i j can be introduced based on neighbor distances or Voronoi tesselations 36 . The l th SOP is computed by combining these terms such that
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Thus, for comparison, for all configurations, we also compute the orientational order parameter X X X OOP (traditionally called the nematic order parameter) and distance-weighted Steinhardt order parameters X X X SOP with freud 34 , using similar l 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 X X X SOP and X X 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 σ GT O = 3.0, n max = 6, l 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. 1 of the supplementary material, 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 interactions 37 , as these aspects (shape, mutual orientation, and translational order) are often interrelated and correlated in the ways in which they influence molecular behavior.
Quantitatively, we can compare the information density of AniSOAP with other order parameters for our configurations using the global feature reconstruction error (GFRE 38 , computed using scikit-matter 39 ). For two representations X X X 1 and X X X 2 of the same dataset, GFRE(X 1 , X 2 ), determines how much information X X X 1 contains relative to X X X 2 , where 0.0 infers X X X 1 can perfectly reconstruct X X 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 construct the analogous X X 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. 2 of the supplementary material, and its ability to reconstruct the AniSOAP representation will decrease when a dataset contains multiple particle types.
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 40 and first introduced by Berardi, Fava, and Zannoni 41 . This formulation calculates the pairwise potential between (potentially dissimilar) ellipsoids, i.e., ellipsoids with three unequal semiaxes: a i , b i , c i . These three semi-axes define a diagonal struc-ture matrix for particle i:
where S S S is analogous to D D D -2 of Eq. ( 12), although (a, b, c) are discrete semiaxes, rather than Gaussian widths, like σ . 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 r r r i and a 3×3 rotation matrix R i , respectively. The relative position is r r r 12 = r r r 2r r 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:
LJ-like term ( 22)
anisotropy corrections (23) The first term U r is given as follows:
where the distance h 12 between two ellipsoids is estimated as the Perram distance of closest approach 42 h 12 (R 1 , R 2 , S S S 1 , S S S 2 , r r r 12 ) = r 12s 12 (R 1 , R 2 , S S S 1 , S S S 2 , r12 ) (25) where
(26) and
The terms ε GB and σ GB are the potential well depth and location analogous to LJ. γ 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. ( 23) account for the ellipsoidal geometry (η 12 ) and misorientation (χ 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 42 is computationally efficient, but open to failure in some cases, like in the case of two bodies with very unequal radii 40 . It is worth noting that it is the nature of density expansions to implicitly contain the distance of closest approach, provided that the This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 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. 3 43 . And thus, in AniSOAP , we alleviate the requirement to explicitly compute h 12 or rely on its approximations.
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. 3 of the supplementary material.
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 PyXtal 45 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 46 using Perdew-Burke-Ernzerhof (PBE) pseudopotentials and cutoff parameters reported by Prandini et al. 47 , a Grimme D3dispersion correction 48 , and a 3x3 Monkhorst-Pack k-point grid 49 . Computations were managed with the signac and signac-flow packages 50,51 , and computed using the ASE computational front-end 52 .
With any new featurization, there is always valid concern about hyperparameter optimization and the sensitivity of the featurization to changes in these values 53 . 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 σ 1 = σ 2 , and prove that, in line with conventional knowledge on benzene geometry, best results are obtained with σ 1 > σ 3 (an oblate ellipsoid).
To tune these parameters, we again use the GFRE of a given AniSOAP representation and an analogous SOAP rep-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. resentation of our benzene crystal configurations. By taking GFRE(X AniSOAP , X 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 structureaveraged representation using rascaline 54 , with a cutoff radius of 7.0Å, Gaussian density width of 0.5Å, l max = 10, and n 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 σ 1 = σ 2 , σ 3 , and σ GTO . Each AniSOAP vector has length 146.
From the results in Fig. 5, we see how the oblate MVG with σ 1 = σ 2 =4.0 and σ 3 =0.5 minimizes the information loss compared to the atomistic SOAP representation with a σ GT O = 1.5. Prolate (σ 1 < σ 3 ) and spherical (σ 1 = σ 3 ) Gaussian densities perform notably worse. It is interesting, however, that there is a wide range of σ 1 , σ 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. This indicates that fitting on AniSOAP vectors is robust to small changes in hyperparameters; without this robustness, models are prone to overfit. We include in the SI similar analyses for σ GTO = {0.5, 1.0, 1.5, 2.0, 2.5, 3.0}, demonstrating that σ GTO ∈ [1.0Å, 3.0Å] obtains similar results to σ GTO = 1.5Å.
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, incomplete 55 with respect to many atom-atom correlations.
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 σ 1 = σ 2 = σ 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 analysis 56 .
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) 57 , has heavily limited regression performance, with a typical RMSE on the order of 45 -50meV/atom, and an R 2 ≤ 0.4. By introducing intermolecular anisotropy, our performance jumps R 2 ≈ 0.9, RMSE ≈ 15 -20meV/atom, much closer to the all-atom regression in the left panel of Fig. 6.
While AniSOAP can achieve reasonable performance, the performance of AniSOAP cannot match that of all-atom SOAP due to the nature of the underlying dataset. In this dataset, we specifically include configurations that will always be challenging to represent via molecular 58 coarse-graining, such as those with unstable hydrogen-hydrogen repulsion, or stabilizing hydrogen-hole interactions. As shown in the dynamic chemiscope 59 visualization in the corresponding Materials Cloud entry 60,61 , these configurations are those with the largest error when regressed using the AniSOAP representation. This shows that the limits of AniSOAP correspond to the limits of coarse-graining itself.
The tools developed under the umbrella of "machinelearning" are powerful, and revolutionizing society and the scientific community at an aggressive pace. One important advancement for chemical sciences is the concept of treating This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 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.
By feeding AniSOAP descriptors into various functional forms to learn molecular energy, we envision these anisotropic descriptors to be foundational for building general coarsegrained potentials of anisotropic molecules, similar to how many machine-learned interatomic potentials are built off atomistic descriptors 16,17,21 . Future work will focus on the efficient implementation of these anisotropic potentials to explore complex, multiscale molecular systems and compare the performance of AniSOAP-based potentials with traditional methods focused on reproducing the potential of mean force. For full integration into molecular dynamics, further derivations will also be necessary to calculate spatial derivatives of these descriptors for use in predicting the forces and torques on each particle.
The supplementary material contains additional analysis on our case studies. For the liquid crystals, we show how AniSOAP represents differently-shaped spheroids, and that it can delineate differently-shaped particles (unlike existing order parameters). For the Gay-Berne and Benzene case studies, we show learning curves and further results of our GFREbased hyperparameter optimization.
Data for this paper will be made available via MaterialsCloud 60 upon acceptance, including all raw input files and analysis scripts. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0210910 VII. 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.
There are no conflicts to declare.
The coefficients ⟨nlm|ρ i ⟩ 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 features 64 . 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 ⟨nlm|ρ i ⟩ to use the unnormalized GTO basis R nl (r) = r l+2n e
and can hence be normalized with the constant
After normalizing all the bases, we can take orthogonalize them using Löwdin Symmetric Orthonormalization 25 . We first find the overlap matrix between two normalized GTOs:
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 Φi : Φi = G -1/2 i j Φ 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.
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. 7.
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 ≤ l 1 , l 2 ≤ l max and |l 1l 2 | ≤ L ≤ min(l max , (l 1 + l 2 )), the number of matrices required to compute grows as O(l 3 max ), 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.
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 replacement 65 . 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 ≤ 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. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 7. 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, σ = 5.0, r cut = 1.0. With caching, after the first computation, all subsequent computations of Clebsch-Gordan matrices are near instantaneous (∼ 10µs, the time required to access a list). The gray dashed line show the expected ∼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.
Given the simplicity of the moments calculations and the frequency of invoking this section of code, we are currently porting this functionality to Rust 66 , 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 D D D (10) is changed. By definition, D D 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.
With the CG-matrix caching and proposed porting of the moments above, we have performed benchmark timings (Fig. 7) 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. Continuing our discussion of Eq. ( 23), the second term, η 12 , is a function of each ellipsoids' geometry and orientation but not a function of distance, and is given as follows:
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
This correction term describes the interaction strength between two ellipsoids at 0 separation, whose influence is tuned by parameter ν.
The last term, χ 12 , is only a function of each ellipsoid's individual orientations (R 1 , R 2 ) and relative orientation ( r 12 ), but not of their geometries (S S S 1 , S S S 2 ). It is given as follows:
B4) E E E i = e -1/µ ai 0 0 0 e -1/µ bi 0 0 0 e -1/µ ci (B5)
Note that for Gay-Berne, the rotations 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, R i,AniSOAP ≡ R -1 i,GB . The above definitions use R i ≡ R i,GB . χ 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 µ = 1:
In our case-study with a i = 1, b i = 1.5, c i = 2 ellipsoids, we set µ = 1, enabling the use of B6 to calculate E E E i . We furthermore set ν = 1, completely specifying the Gay-Berne hyperparameters. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE ASDOI: 10.1063/5.0210910