Completeness of Atomic Structure Representations

In this paper, we address the challenge of obtaining a comprehensive and symmetric representation of point particle groups, such as atoms in a molecule, which is crucial in physics and theoretical chemistry. The problem has become even more important with the widespread adoption of machine-learning techniques in science, as it underpins the capacity of models to accurately reproduce physical relationships while being consistent with fundamental symmetries and conservation laws. However, some of the descriptors that are commonly used to represent point clouds -- most notably those based on discretized correlations of the neighbor density, that underpin most of the existing ML models of matter at the atomic scale -- are unable to distinguish between special arrangements of particles in three dimensions. This makes it impossible to machine learn their properties. Atom-density correlations are provably complete in the limit in which they simultaneously describe the mutual relationship between all atoms, which is impractical. We present a novel approach to construct descriptors of \emph{finite} correlations based on the relative arrangement of particle triplets, which can be employed to create symmetry-adapted models with universal approximation capabilities, which have the resolution of the neighbor discretization as the sole convergence parameter. Our strategy is demonstrated on a class of atomic arrangements that are specifically built to defy a broad class of conventional symmetric descriptors, showcasing its potential for addressing their limitations.

In this paper, we address the challenge of obtaining a comprehensive and symmetric representation of point particle groups, such as atoms in a molecule, which is crucial in physics and theoretical chemistry.The problem has become even more important with the widespread adoption of machinelearning techniques in science, as it underpins the capacity of models to accurately reproduce physical relationships while being consistent with fundamental symmetries and conservation laws.However, some of the descriptors that are commonly used to represent point clouds -most notably those based on discretized correlations of the neighbor density, that underpin most of the existing ML models of matter at the atomic scale -are unable to distinguish between special arrangements of particles in three dimensions.This makes it impossible to machine learn their properties.Atom-density correlations are provably complete in the limit in which they simultaneously describe the mutual relationship between all atoms, which is impractical.We present a novel approach to construct descriptors of finite correlations based on the relative arrangement of particle triplets, which can be employed to create symmetry-adapted models with universal approximation capabilities, which have the resolution of the neighbor discretization as the sole convergence parameter.Our strategy is demonstrated on a class of atomic arrangements that are specifically built to defy a broad class of conventional symmetric descriptors, showcasing its potential for addressing their limitations.
The construction of data-driven models of physics starts by defining a suitable mathematical representation of the problem, which is then used as the input to different types of model architectures [1].In many cases, this involves describing the spatial arrangement of point particles [2,3] -e.g.atoms in the construction of microscopic models of matter -in a way that is consistent with the basic physical symmetries governing the relationship between a structure and its properties [4,5].It would be desirable to determine a complete set of invariants, i.e. a set of geometric descriptors that are equal for two configurations if and only if these configurations are the same, modulo translations, rotations, inversion, and permutations of equivalent atoms (see also Sec. 5 of Ref. 6 for a discussion of different definitions of completeness).Unfortunately, it was recently shown that many popular structural representations are incomplete, [7,8] leading to the outright impossibility of predicting different values for the properties of structures even though they are unrelated by symmetry, and more broadly to numerical instabilities that extend to nearby configurations [9,10].
Even though the significance of this problem for datadriven atomistic modeling has been recognized only recently, the problem of describing an arrangement of particles based on limited and/or symmetrized information is a long-standing one [11].This issue is the subject of a substantial body of literature, touching upon open problems in invariant theory [12] and the Euclidean distance geometry problem [13][14][15], with applications to fields as diverse as signal processing and the determination of protein structure by NMR.The broad connection to geomet- * michele.ceriotti@epfl.chric deep learning and graph convolution [2], as well as to the construction of equivariant physical models [5] is also well-understood.In the application domain that we focus on here, symmetry has been widely used since the early efforts to construct machine-learning interatomic potentials [16][17][18], and it also underlies efforts geared towards recognizing the (dis)similarity between structures in simulations and databases [19,20].
A systematic study of the representations used to characterize atomic structures [6,21] has revealed that many of these frameworks rely on the same physicsinspired descriptors to characterize groups/clusters of atoms relative to a central one: discretized versions of structural correlation functions of the neighbor density with increasing body order [22] -which are equivalent to the unordered list of interatomic distances, angles, dihedrals, relative to the central atom.While it was shown that keeping all orders of correlations provides a complete, symmetry-adapted linear basis to expand structure-property relations [23,24], the low-order versions -that are most commonly adopted in practicecannot distinguish groups of degenerate configurations [7].Similar problems affect other machine-learning architectures, such as distance-based graph convolutional networks [8], and even equivariant neural networks [4,25], that are, in theory, universal approximators, but rely on constructing a hierarchy of symmetric polynomials that are closely related to the density-correlation features discussed above [26,27].Completeness, in this context, is only guaranteed in the limit of infinitely deep networks, which makes it harder to allocate computational effort between different architecture parameters, such as network depth and basis resolution [24].In this Letter, we discuss a construction that relies only on correlations between a finite number of neighbors to define invariant descriptors of a local atomic environment, that are complete in the sense that they can be used to build simple universal approximators for atomic properties and to differentiate between arbitrary configurations while scaling favorably with respect to system size and allowing for efficient implementations.
To appreciate the significance of this result without delving into the theory of density-based representations (interested readers may refer to the discussion in Ref. 6, or to the summary in the SM) it is useful to consider the concrete case of two-neighbor representations (Fig. 1a).To evaluate a property y(A) of an atomic structure A, it is convenient to express it as a sum of atom-centered contributions y(A i ), that in turn can be expanded into body-ordered terms y (ν) (A i ) depending simultaneously on the relative position of ν of the neighbors of the central atom i.In line with previous literature, we indicate terms associated with correlations between ν neighbors as having body order (ν +1).In the case of a scalar target (e.g. the interatomic potential) the two-neighbor term can be expressed as a sum of a function of two distances and one angle, y (2) (r i1i , r i2i , ri1i • ri2i ) over all neighbors pairs.Alternatively, and equivalently [21,23,24,28], we can expand the histogram of (r, r ′ , θ, ) on a discrete basis ⟨q|, and approximate y (2) (A i ) as a linear combination of these coefficients, q y (2) q ⟨q|ρ ⊗2 i ⟩.Many widespread descriptors, from atom-centered symmetry functions [29] to the SOAP power spectrum [30] can be written in one of these forms.The resolution of the basis is a convergence parameter that in practical applications [24,31,32], and in particular when using the descriptors as inputs of a non-linear model [33], can be usually converged with ease.Higher-order terms can be obtained by computing (ν + 1)-body symmetric density correlations -that we indicate with the shorthand |ρ ⊗ν i ⟩.This body-ordered expansion of y(A) converges systematically with ν but is rather cumbersome.Instead, many frameworks rely on low-ν descriptors and use them as input of artificial neural networks [16].The non-linear functional form generates some, but not all, the terms with higher body order [26,34], leading to more flexible models and making more effective use of the features restricted to a finite body order ν.
Finite body-order descriptors, however, have a fundamental limitation: low-ν invariants are not complete, in a way that cannot be remedied by non-linear transformations.For example, if one uses only interatomic distances ( |ρ ⊗1 i ⟩), it is easy to realize multiple structures with neighbors placed in arbitrary directions but the same list of distances {r ji } (Fig. 1b).All these structures will have identical |ρ ⊗1 i ⟩ features, and therefore also the same output for any non-linear model built on them.Less trivial constructions have been found that lead to degenerate pairs of configurations for triangle-histogram features (the power spectrum |ρ ⊗2 i ⟩) as well as for tetrahedrahistogram features (the bispectrum |ρ ⊗3 i ⟩) [7].Counterexamples for higher ν are not known, but there is no proof of completeness either: for a general function no finite order ν of rotationally symmetric correlations suffices for a complete reconstruction [35].The existence of degeneracies means that, for finite ν, there are pairs of distinct structures for which it is impossible to build an ML model that predicts the differences in their properties, regardless of the resolution of the descriptors.It also determines numerical instabilities in their vicinity, and therefore models built on these incomplete features can suffer from degraded accuracy and transferability [10,36,37].
To avoid these fundamental and practical flaws we want to construct a framework based on neighbor densities that is provably complete while using a finite order of correlations between neighbors.To do so, we use an extension of the density-correlation framework to (N +1)center quantities, i.e. descriptors of N -tuples of atoms in the environment A i , that are invariant to the ordering of their neighbors, but equivariant to the labeling of the N tagged particles [38].The central idea in our construction is to build invariant descriptors representing the histogram of one neighbor (ν = 1) in the coordinate system defined by the central atom and two tagged neighbors (N = 2), translating complete frameworks based on local coordinate systems [39][40][41][42][43] in the language of density correlation descriptors.Using Dirac's bra-ket notation to enumerate the components of the descriptor, we compute features of the following form: This expression involves an expansion 2. (a) 3-center-1-neighbor features describe the position of atoms relative to two tagged atoms around the center i.This is equivalent to a list of tetrahedra sharing the ii1i2 triangle.(b) Tagging two atoms defines a local orthogonal coordinate system -except when the tagged atoms are collinear, and the azimuthal directions are ill-defined.(c) An encoder-decoder architecture can be used to determine nonlinear environment features ξb that contain enough information to reproduce high-order representations.
R nl (r ji )Y m l (r ji ), in radial functions R nl and spherical harmonics Y m l , of the position of neighbors i 1 and i 2 .⟨nlm|ρ i ⟩ = j ⟨nlm|r ji ⟩ is the expansion of the neighbor density on the same basis.The contraction with Clebsch-Gordan coefficients ⟨l 1 m 1 ; l 2 m 2 |l 3 m 3 ⟩ selects the rotationally invariant parts.The first ingredient in our approach is the realization that Eq. (1) provides a complete description of the environment A i .A detailed proof is reported in the SM, but the reasoning can be explained in relatively simple terms.As long as atoms do not overlap, the peaks in the neighbor density precisely encode information on the position of all the atoms.Averaging |ρ i ⟩ directly over rotations discards angular information.However, the presence of descriptors of the position of two tagged neighbors (i 1 , i 2 ) means that rotational averaging of the entire tensor product preserves the full information on the position of all neighbors relative to these two atoms (Fig. 2a).In other terms, the descriptors in Eq. ( 1) are equivalent to knowledge of the position of the neighbors in the coordinate system defined by (r i1i , r i2i ).Provided that r i1i and r i2i are not collinear, it is possible to reconstruct the position of all the atoms within the environment relative to the center i, modulo rotations.Thus, by feeding a sufficiently fine-grained discretization of |ρ ⊗1 ii1i2 ⟩ to a universal approximator, any function of the coordinates of the atoms in A i can be approximated in a way that is invariant to rotations, and to neighbor permutations, except for the tagged atoms i 1 and i 2 .As we shall see, the convergence with the number of three-center descriptors and non-linear features can be achieved easily.In order to restore invariance and recover a description of the single-center environment, we sum over these two indices.Performing this sum for Eq.(1) recovers the bispectrum |ρ ⊗3 i ⟩ (recall that j |r ji ⟩ = |ρ i ⟩) that we know to be incomplete [7].Instead, we apply a non-linear layer on top of the three-center (N = 2) features before summing over i 1 and i 2 .Information on individual pair-dependent features is discarded, but any property of the environment A i can be learned, provided that one uses a sufficiently expressive architecture.Defining the shorthand ξ q (A ii1i2 ) = ⟨q|ρ ⊗1 ii1i2 ⟩, where q is a compound index that encompasses the labels (n 1 , l 1 , n 2 , l 2 , n 3 , l 3 ), we can write the environment features ( If one chooses ψ p to be a scalar universal approximator, Eq. ( 2) implements a deep set architecture [44], which is known to provide permutation invariant universal approximation.Combining information on all (i 1 , i 2 ) pairs within A i also ensures completeness in the presence of collinearity of certain pairs.For pairs corresponding to collinear vectors, |ρ ⊗1 ii1i2 ⟩ cannot resolve the position of the neighbors along the azimuthal direction (Fig. 2b).However if some of the neighbors in A i are not collinear, there will be at least one pair that forms a basis for three dimensions and provides a complete description, and the compound features ξ incorporate this information.If all the atoms are aligned along the same axis, then there is no azimuthal information to be encoded, and the description remains complete.
It is important to note the relationships between this procedure and several ideas that have been used in the literature, built around the definition of a rigid coordinate system relative to which one can build local (or global) geometry descriptors [39,45], or to that of summing over all (or some) neighbor index permutations, a representation that is invariant to rotations but not to neighbor ordering [46,47].These approaches typically suffer from a lack of smoothness, because the procedure that determines the reference environment, or the index permutations that should be included in the summation, leads to discontinuous variations as a structure is deformed continuously.The recently-introduced Weighted Matrices Invariant (WMI) framework [42,48] is the closest to the spirit of our construction.It uses the unordered set of local atomic coordinates as an invariant and provides a practical algorithm to compare structures using a suitably defined smoothly-varying metric.Our approach, instead, effectively defines body-order-complete, smooth, invariant descriptors that can be readily understood and used in relation with several popular schemes based on discretized correlations of the neighbor density.A more in-depth comparison of the current approach with alternatives can be found in the SM. and 'NL' linear and non-linear layers) for learning the high-ν features, and the energy of a validation set of 500 bispectrumdegenerate B8 pairs.We separately report the error on the mean of the features and energies of degenerate pairs (A + , A − ), and on their difference, normalized by the standard deviation in the validation set -i.e. for the error on the mean ϵ(⟨E±⟩) 2 = 1

Model Relative Validation Error Name
We report the error relative to the standard deviations σ( 12 (E(A + ) + E(A − ))) and σ(E(A + ) − E(A − )), respectively.Absolute energy errors can be obtained by multiplying the values by 4.26 eV (Mean A ± ) and 0.245 eV (Diff.A ± ).
The expression in Eq. ( 2) does in principle have sufficient flexibility to fit arbitrary symmetric functions of the atoms within A i .For example, a single universal approximator ψ would suffice to make a complete predictor of a (short-ranged) interatomic potential, that would effectively be obtained as a sum of triplet energies.From this point of view, it is interesting to note the excellent performance that has been recently shown [49] for ML potentials that are built as a sum of pair terms, which goes one step towards the provably-complete framework we propose here.It is easy to see, however, that a single feature cannot provide a complete description of an environment: a simple counting argument implies that at least (3n i − 6) features are needed to completely describe the degrees of freedom of an environment that contains n i atoms.We propose to use an encoder-decoder architecture to generate a set of non-linear features ξ that can predict a large set of higher-order invariants (ideally with ν approaching the number of neighbors), as a way to generate a complete description (Figure 2b).
To demonstrate this idea, we take the family of degenerate pairs of environments from Ref. 7 that are indistinguishable with respect to all atom-centered invariants up to the ν = 3 bispectrum.This construction involves 7 neighbors around a central atom and generates a manifold with six continuous geometric parameters, and can be decorated with up to three atom types.For simplicity, we realize it with a single element, choosing boron clusters B 8 that are known to exhibit remarkable structural diversity [50,51].By varying the parameters in an appropriate range, we generate a dataset of 4000 pairs, computing their energy with PBE-DFT in the cc-pVDZ basis, using PySCF [52] (details and resulting structures are given in the SM).We then compute both 3-center-1-neighbor descriptors of A ii1i2 triplets, as well as atomcentered invariant correlation features of order ν = 7, based on the NICE framework [53].We emphasize that -given that the degeneracy affects the environment and is resolved by including multiple centers -all descriptors and models are built for a single atom at the origin.We build two simple encoder/decoder architectures that combine 3-center descriptors, contract them into environment descriptors, and then apply a multi-layer perceptron (MLP) with SiLU activation to approximate the high-order features.One architecture (A L model) uses only a linear transformation before summing over (i 1 , i 2 ), while the other (A NL model) also includes a non-linear layer (another MLP) before the contraction.Details of the architectures are given in the SM, but are not essential to our findings.We use the intermediate ξ(A i ) features in the two models to learn the energies of the B 8 clusters.We also build "direct" models for the energy as a sum of linear and non-linear triplet energies (B L and B NL ), full-body-ordered atom-centered features (C L and C NL ), and bispectrum features (D L and D NL ).Table I summarizes the validation set accuracy of all these architectures.The accuracy in reproducing the full-bodyorder |ρ ⊗7 i ⟩ correlation features can be interpreted as a measure of the resolving power of features, while the ac-curacy in predicting energies measures the impact of lack of resolution on a common regression task.
As expected, all models that rely on a linear contraction of the 3-center descriptors |ρ ⊗1 ii1i2 ⟩, being equivalent to models based on bispectrum, are incapable of resolving the degenerate pairs, leading to a fractional error of 1 on the difference of features and energies.The lack of resolving power, and the consequent loss of regression accuracy, are also evident in the sensitivity analysis in (Fig. 3a), and in the energy prediction parity plot (Fig. 3b).Even though the performance of individual models depends in a non-trivial way on the model details and the train-set composition and size, one can make a few observations based on the fact that we kept a consistent architecture for the non-linear layers among all models.For starters, the relative errors on feature and energy differences among degenerate pairs are large also for the models that have -in principle -universal approximating power.This is due in part to the fact that these differences are small compared to the overall variation in energy, and that the magnitude of correlation features is dominated by low-order terms.Furthermore, the practical model accuracy depends on the fine-tuning of the architecture, which we deliberately avoided to facilitate a qualitative comparison between models.For example, one can clearly see that the NICE discretization of |ρ ⊗7 i ⟩ is not sufficiently converged to provide a complete linear basis: the C L model has a validation set error of almost 90% for the energy difference between pairs.As for many frameworks in the literature, this lack of basis set convergence can be mitigated by applying a non-linear layer to the high-order features: both the C NL and the A NL models (the latter being computed on top of the latent features that allow the reconstruction of |ρ ⊗7 i ⟩ to a very good accuracy) achieve an error on the energy differences around 45%. From this point of view, the fact that a latent layer with only 128 features performs similarly to the full set of |ρ ⊗7 i ⟩ (that is iteratively contracted to hold ≈8000 features) underscores the effectivenes of non-linearities in tackling issues of basis set convergence.The error in approximating |ρ ⊗7 i ⟩ converges quickly with latent-space size (see SM), which is consistent with the notion that a minimal set of symmetric descriptors does not have to be much larger than the number of degrees of freedom in the atomic environment being described, even though we cannot provide a rigorous upper bound to the number of non-linear features needed to be complete.Finally, we observe that the B NL model based directly on non-linear functions of the 3-center features performs much better than all other options -not only in terms of accuracy in describing pair differences, but also the mean energy of each pair.The impact of degeneracies on realistic data sets, that do not explicitly include them, depends on the details of the model [37], and previous studies of power spectrum degeneracies showed that a clear effect is only observed in the limit of very large datasets [7].
Summarizing our findings, we propose a practical solution to the open problem of determining a complete local representation of an atomic structure (more generally, a 3D point-cloud) based on a finite-dimensional and low body-order O(3)-invariant description of an environment.We demonstrate the idea on an artificial, but very challenging, set of B 8 cluster configurations, designed to defy linear atom-centered descriptors of the same body order.The accuracy and computational cost of an ML model in real-life scenarios depend on many different considerations.This approach indicates a way forward to designing ML frameworks that have universal approximation power while being smooth, symmetric, and free from the restrictions of architectures based on O(3) equivariant layers (see e.g.Ref. 43 for a very recent example of a model that achieves excellent performance, enforcing equivariance by averaging over an ensemble of coordinate systems).The close connection with established, well-understood body-ordered descriptors means that the ideas we expose here can be easily incorporated into existing frameworks, such as those based on atom-centered symmetry functions or those that build non-linear models based on invariant density-correlation features, which can now be safely truncated to a fixed and low order by evaluating features and models at the level of neighbor triplets.

I. SUMMARY OF DENSITY-BASED DESCRIPTORS
A structure A containing n at independent atoms with positions {r i } and chemical types {a i } can be described within the framework of atom-centered density correlations (ACDC) as follows.As the name suggests, this family of descriptors or representations decomposes the structure into environments A i = {a i , {(a j , r ji ) : r ji < r cut }} centered on each atom, consisting of neighboring atoms (expressed by their inter-atomic distance vectors r ji = r j − r i and species a j ) within a certain distance r cut away from the center.Usually, the neighbors or equivalently their separation from the center r ji is captured through a localized function g(x − r ji ) that peaks at the position of the neighbor, for example, g could be a Gaussian function or a Dirac-delta distribution (when this is the case we use the notation g → δ).
This description of environments in terms of interatomic displacements is invariant to the choice of an arbitrary origin (equivalently, invariant to rigid translations).By invariance to a transformation, we mean that the description can be written as a function that remains unchanged when the transformation is applied to the arguments/variables.Historically, these descriptions are then made independent of the order in which the neighbors are enumerated.This is easily achieved by an "invariant pooling", e.g. a summation over all neighbors in the environments A i , resulting in a translation and permutation invariant atomic neighbor density ρ i (x).In the following discussion, we will use Dirac's bra-ket notation to capture a broad class of descriptions (such as SOAP, ACE, FCHL, etc.) following a similar procedure, and attribute the specificities of each representation to the choice of basis and the localized function used to enumerate the neighboring atoms.So g(x − r ji ) ≡ ⟨x|r ji ⟩, and ρ i (x) ≡ ⟨x|ρ i ⟩.It is easy to see that in the g → δ limit this representation provides a linear basis to expand properties that depend separately on the position of individual neighbors, y (1) Notice that the description of each atomic environment still depends on the orientation of the reference system chosen to describe the coordinates r ji of the neighbor j as the translational symmetrization only fixed the origin to the central atom.However, one might be interested in modeling scalar properties such as energies that are invariant to global rotations of the structure, or alternatively target an equivariant property such as the dipole moments that have a specific rotational behavior.To capture the behavior of the target with respect to rotations (or equivalently any other symmetry operation) of the structure, the structural descriptions can be endowed with the very same behavior by integrating over the corresponding symmetry operations.
For instance, in the case of rotational equivariance, to integrate the neighbor density over rotations, it is easier to first write it on a basis of spherical harmonics ( |lm⟩) that have a precisely determined rotational dependence.This change of basis from ⟨x| to ⟨nlm|, where ⟨n| is a discretized description of the absolute distance |x| and ⟨lm| captures the direction x, can be written as ⟨nlm|r ji ⟩ ≡ dx ⟨x|r ji ⟩ ⟨nl|x⟩ ⟨lm|x⟩ (S4) so that ⟨nlm|ρ i ⟩ = j∈Ai ⟨nlm|r ji ⟩ A representation with the desired symmetry (of rotation like the spherical harmonic |λµ⟩) is then achieved by in which we emphasize the SO(3) character of the representation by indicating λµ in the ket.Following this procedure, we have obtained a symmetric (equivariant) description of the target property in terms of contributions of individual neighbors.Eq.S3 still holds in the symmetrized form with ⟨n|ρ ⊗1 i ; λµ⟩, providing a linear basis to expand a function with the symmetry of Y µ λ that depends on the coordinates of individual neighbors.This procedure can be easily extended to body-ordered expansions of atom-centered properties [21,23,28,34], that use the tensor products of |ρ i ⟩ as a linear basis to express quantities that simultaneously depend on the position of ν neighbors.We use the shorthand ⟨q|ρ ⊗ν i ; λµ⟩ to indicate such a description that is also equivariant.For more details on these atom-centered descriptions, we advise the reader to consult Refs.6, 21.Notice that in constructing |ρ i ⟩, we summed over all the neighbors in the environment, thereby losing the identity of an individual atom in the neighborhood, but retaining the identity of the single central atom i.An equivalent way of saying this is that |ρ i ⟩ is a histogram of neighbor positions in A i as it describes the distribution of the neighbor displacements without their explicit identities.Similarly, for |ρ ⊗ν i ; λµ⟩, we started with a description of individual atoms that was then symmetrized over translations, permutations, and finally over rotations.This representation describes the histogram of correlations of the central atom with ν neighbors.One might be interested in descriptions of atomic environments retaining the identity of more than one atom, such as pairs, triplets or N-tuples of atoms.We proposed an extension [38] of the ACDC formalism reviewed above to tag N other particles to obtain a representation which explicitly addresses the identity of each of these atoms |ρ ⊗ν ii1i2...i N ; λµ⟩.
These N-centered representations are translationally invariant and are also symmetrized in the same order-they are invariant to the permutation of the ν neighbors but equivariant to the permutation of the labelling of the N -centers; and are rotationally symmetrized the same way as their one-center counterpart.We direct the interested reader to see the accompanying SI of Ref. 38 for a detailed derivation.Note that with an N -centered description, we effectively describe the correlations of 1 central atom with N -1 other centers and ν neighbors.One can restore the N-centered representation to an atom-centered one by summing over the N-1 additional centers This provides a way one can promote neighbors to special centers or demote special centers to neighbors.Following the discussion above, a three-centered 1 neighbor representation is comparable to an atom-centered bispectrum, both of which describe correlations of one central atom with 3 other atoms, thus having similar descriptive power in terms of the correlation order captured.However, the former singles out indices of two of the neighbors while the latter is permutation invariant to all of the indices but the center i.

II. PROOF OF COMPLETENESS
In this section, we provide a proof of the central result of the main text, namely that any continuous function f (r i1i , r i2i , . . ., r i N i ) defined within the local neighborhood of a center atom i that is invariant under rotations and permutations can be approximated to arbitrary accuracy using the features shown in Eq. (2), The 3-vectors (r i1i , r i2i , . . ., r i N i ) indicate the positions of N particles in B3 rcut = {r ∈ R 3 : ∥r∥ ≤ r cut } relative to some center atom i (r ji = r j − r i ) with coordinates given in an arbitrary coordinate system, and ⟨nlm|□⟩ indicates discretization on a basis of radial functions and spherical harmonics.The proof rests on two key steps: (1) using relative coordinate frames for a complete rotationally invariant description and (2) constructing a permutationally invariant function using these inputs.We devote subsections II A and II B to the two steps, and discuss in Section III the connections to other approaches.

A. Complete Rotationally Invariant Descriptors using Relative Coordinate Frames
For the first step, we will follow the presentation in Ref. 42 and use, for simplicity, the notation (r i1i , . . ., r i N i ) to indicate the position of the neighbors relative to the central atom.Let us first assume that the first two particles at positions r i1i and r i2i are not collinear.These then define a new orthonormal basis (v 1 , v 2 , v 3 ) of R 3 by taking the first two basis vectors from the Gram-Schmidt orthonormalization procedure, i.e. by setting v 1 = ri 1 i ∥ri 1 i∥ and the second basis vector v 2 to be the normalized version of r i2i − (r i2i • v 1 )v 1 , so that it lies in the plane spanned by r i1i and r i2i .Finally, we set v 3 = v 1 × v 2 leading to a unique right-handed coordinate system.The coordinates of an arbitrary point r ji with respect to this relative basis are simply given by the inner products r ji •v i .In particular, these coordinates are invariant under rotations R ∈ SO(3) as this reference frame rotates with the underlying structure.(A similar procedure is adopted in Ref. 55 for describing atomic environments.)Apart from the non-collinearity assumption, which will be addressed shortly, we can repeat this construction by choosing any pair (r i1i , r i2i ) to form the basis vectors from, leading to N (N − 1) rotationally invariant ways to describe the position of all the particles, which (after considering permutation invariance) maps to the weighted matrix invariants (WMIs) that have been proven to provide a metric that is able to differentiate between any pair of atomic structures [42].
We now show that this is precisely the information contained in the 3-center descriptors in the main text.Consider the function (strictly speaking, a tempered distribution) which peaks whenever (x, x ′ , x ′′ ) are arranged in the same way as (r i1i , r i2i , r i3i ).By averaging over all rotations, we obtain precisely the 3-center feature for the special case of infinitely sharp δ-like densities (i.e.g → δ limit of atomic densities) There is a single peak whenever (x, x ′ , x ′′ ) are arranged as (r i1i , r i2i , r i3i ), up to a rotation.In particular, |ρ ii1i2i3 ⟩ determines the position of r i3i relative to r i2i and r i1i , just as the basis we have constructed before.Summing over the third index i 3 , we obtain This is essentially a histogram of the relative position of all the neighbors, that has (provided the atomic density peaks are sufficiently sharp, and the resolution fine enough) complete information on the position of all atoms within the environment A i relative to i 2 and i 1 .A projection onto a complete discrete basis preserves this information, which proves the completeness of ⟨n 1 l 1 ; n 2 l 2 ; n 3 l 3 |ρ ⊗1 ii1i2 ⟩.The size of this basis is a further convergence parameter, that is common to all atom-density correlation descriptors.
To fill in a final gap, we briefly address issues related to collinearity of the particles.Since we consider all (r i1i , r i2i ) pairs in the environment around i (A i ), if only some of the particles are collinear but A i still contains at least one non-collinear pair (r i1i , r i2i ), the coordinate frame relative to these two atoms provides a complete description of the structure.A more formal proof that at least one such pair suffices can be found in Ref. 42.On the other hand, if all atoms in an environment are collinear, the environment can be completely reconstructed from the set of pairwise distances.This information is still contained in the descriptor |ρ ⊗1 ii1i2 ⟩ due to the location of the peaks so that the 3-center descriptor allows reconstruction of the structure up to rotations.

B. Permutation Invariance
Recalling that our goal is to have a general expression for a function f (r i1i , r i2i , . . ., r i N i ) that is invariant under rotations and permutations, we still need to provide the missing link from the rotationally but not yet permutationally invariant representation |ρ ⊗1 ii1i2 ⟩.The crucial step is to perform a transformation that symmetrizes over the N (N − 1) possible choices of the two vectors i 1 and i 2 that form the reference basis.
For the weighted matrix invariants (WMIs) in Ref. 42, the permutation invariance is taken into account by simply considering the N (N − 1) coordinate systems to be unordered when comparing two structures.In practical terms, checking whether the ordered tuple (x 1 , . . ., x N ) is equal to (y 1 , . . ., y N ) requires N comparisons at most while checking whether two unordered sets {x 1 , . . ., x N } and {y 1 , . . ., y N } are equal naively suffers from an additional cost of checking through N !ways to rearrange the elements.While Ref. 42 proposed more efficient methods to test the equality of two unordered sets, the methods are not directly translatable to common problems in condensed matter physics and related fields, in which we want to specify the descriptors for one structure rather than to compare two configurations.
At the same time, properties of functions invariant under permutations of the input variables have been studied in a more general context [44].In particular, rephrasing the results for the case of three center representations, it has been shown (theorem 9 in Ref. 44) that any permutationally invariant continuous function f perm that takes a discretization of the three-center features ξ q (A ii1i2 ) = ⟨q|ρ ⊗1 ii1i2 ⟩ as inputs can be approximated arbitrarily well by functions of the form where F now is an arbitrary continuous function, and the ψ p are (non-linear) maps that transform the initial features, which are then summed over the indices i 1 and i 2 to generate the permutationally invariant intermediate features ξ = ( ξ1 , . . ., ξp , . . .).These precisely correspond to the components of the encoded intermediate features in Eq. ( 2) in the main text.This completes the proof that any function f (A i ) that depends on an atomic environment can be approximated arbitrarily well by the presented approach.

III. CONNECTION TO RELATED RESULTS
In this section, we summarize the connection between the three center descriptors proposed in this work and alternative viewpoints on the problem of achieving a complete description of an atomic environment.

A. Connection to the Gram matrix and Weighted Matrix Invariants
Let us denote the N displacement vectors by (r i1i , . . ., r i N i ).If we only require descriptors that are invariant under the action of some improper rotation R ∈ O(d), the N × N Gram-matrix defined by is well known to be an overcomplete set of invariants (see e.g.Ref. 5 for how this classical result can be used in an ML context).
In fact, for d = 3, if the first three vectors r i1i , r i2i , r i3i are not coplanar, the first three rows (or equivalently, columns) of the Gram-matrix suffice, hence the overcompleteness.While we need three (non-coplanar) rather than two (non-collinear) vectors here, this is consistent with the previous observation that we can reconstruct the entire structure from the relative coordinates of only r i1i and r i2i .This is because the third of the three basis vectors (v 1 , v 2 , v 3 ) mentioned in section II A is constructed using the first two.Thus, being able to take the cross-product to upgrade the vectors to a basis is what makes two vectors sufficient for a complete description of an environment.For permutational invariance, the purely theoretical approach would be to simply consider the Gram-matrix up to permutation of the rows and columns, which is closely related to the WMIs in [42] and other similar constructions [48], that however avoid enumerating all permutations.Regardless of whether one uses inner products (Gram matrix) or the relative coordinates (WMIs) between the atomic positions, these sets of invariants are quite different from many of the more commonly used descriptors, in the sense that one is working with unordered sets.While defining a function on an unordered set is simple mathematically, a practical algorithm will in many cases require the storage of the elements in an ordered manner.Explicitly parametrizing a function on an unordered set would then require an additional specification of how to enforce the unordered nature.One such example is the summation over all permutations of indices, as is done with the invariants in this work.Thus, one cannot just use these "invariants" as naively as many of the other "descriptors" and, for instance, feed them directly into a neural network.These invariants are nevertheless useful to perform comparisons between structures due to the existence of a metric [42,48].
A further issue with using invariants such as the unordered set of Gram-matrices or the WMIs is the fact that the number of features changes with the number of neighbor atoms.In principle, the two aforementioned methods could either be used to describe the structure as a whole, leading to a computational cost scaling as O(N 2 ) in the number of atoms N , or be restricted to a finite environment using a cutoff radius.The second approach would be more favorable in terms of scaling, which then however introduces an extra complication due to discontinuities associated with atoms leaving or entering the local environment.In this case, not only the invariant matrices per se, but also their shape and number would change discontinuously.
As mentioned in the main text, similar approaches have been used to introduce approximate permutation invariance for models based on the distance matrix or its non-linear transformations.[46,47] This idea is also related to approaches, such as Social Permutation Invariant coordinates [56] or the Overlap Matrix fingerprints [19] that compute non-linear functions of the distance matrix and then build descriptors from the sorted eigenvalues of the transformed matrix.While there is no proof that these descriptors are complete, they are empirically capable of differentiating between the known degenerate structures [9].

B. Connection to other ML schemes in atomistic modeling
Within the subfield of condensed matter physics and materials science, there are two other approaches providing proofs of complete approximation capability starting from invariant descriptors.The first such approach was the moment tensor potential (MTP) [23] proposed in 2016, which uses the Cartesian coordinates of the particles, followed by the atomic cluster expansion (ACE) in 2019 [24] which can be seen as a related approach relying more on spherical coordinates.We focus in particular on ACE and the mathematically equivalent methods including the NICE framework [53] which are closer in spirit to our approach.All such models start by constructing an atom-centered neighbor density |ρ i ⟩ around some atom i, either with Gaussians or infinitely sharp delta-functions, which by itself contains complete information about the atomic positions (for sufficiently sharp densities, the atomic positions are precisely the location of the peaks of the density function).As explained in section I, this representation is already permutationally invariant, but not yet invariant under rotations.In practice, the densities are projected onto a basis, the coefficients of which ⟨nlm|ρ i ⟩ are used as inputs to further models.If a sufficiently large basis and sufficiently many coefficients of the expansion are retained, these still allow us to completely reconstruct the position of all atoms in the neighborhood.Methods like ACE and NICE now use a linear model on all possible rotationally invariant combinations of these coefficients, which is mathematically equivalent to taking an arbitrary rotationally invariant (nonlinear) function of the inputs f rot (⟨nlm|ρ i ⟩), which can be written in a form similar to eq. (S12) using the shorthand ζ q (A i ) = ⟨q ≡ (nlm)|ρ i ⟩ as where F now is an arbitrary continuous function, and the φ p are (non-linear) maps that transform the initial features, which are then averaged over all possible rotations to generate the rotationally invariant intermediate features ζ = ( ζ1 , . . ., ζp , . . ., ).In practice, these intermediate features would correspond to the ACE/NICE features that are used in the expansion.In our approach, we start from a density |ρ ⊗1 ii1i2 ⟩ that is invariant under rotations only, and then apply a generic permutationally invariant nonlinear function f perm |ρ ⊗1 ii1i2 ⟩ .At the purely conceptual level, apart from switching the orders between permutation and rotation symmetries, we can see that the two approaches are quite similar in spirit.
One key difference is that in many of the practical implementations including the original formulations of ACE and NICE, the body-order expansion is truncated at a finite order ν.Thus, while complete in the limit of arbitrarily high body orders, the expansion is incomplete in practice.Applying a nonlinear model would technically generate higher order invariants, but in order to preserve a complete description, the model would still require that the atom-centered invariants are complete for some ν < N , which is at present unproven.From a mathematical point of view, this difference regarding implementations can be seen as a consequence of the fact that generating a generic permutationally invariant function, having the form shown in Eq. (S12), is significantly simpler than a generic function having rotational invariance.While we can apply arbitrary nonlinear functions and perform a sum over indices, rotationally invariant combinations are restricted to those obtained from a careful study of irreducible representations of the rotation group.The closest conceptual analogue to our approach would be to only take algebraically independent coefficients in the ACE/NICE expansion and use those as inputs into a nonlinear model.The MTP approach essentially performs both steps at once, by constructing features that are invariant under both rotations and permutations from the get-go.More precisely, they are suitable rotationally invariant contractions of the atomic coordinates summed over all atoms in the structure, which can also be related to ideas of deep sets applied, with some modifications, to components of the Gram-matrix.
For all of the methods presented above, it should be noted that completeness for an arbitrary number of particles N (excluding the center atom) is only guaranteed in a certain limit.For MTP and ACE/NICE, this is more directly seen as the maximal order of the used expansion, as well as (for ACE/NICE) the size of the basis used in the expansion.Even though a large number of features is usually necessary when these descriptors are used in a linear model, far fewer terms are necessary to achieve the notion of completeness that is used in this work, as is underscored by the aggressive truncation that can be used when introducing nonlinearities or message-passing operations in the architecture of related models e.g. in MACE [57], NeQUIP [33] or PAINN [58].For our invariants starting from the three center features, this truncation is hidden in the number of required encoded intermediate features ξ p as well as in the resolution of the discretization of the 3-centers-1-neighbor features.Considering that N atoms in an environment have 3N − 3 degrees of freedom modulo rotations, any complete model would need at least as many coefficients.In practice, apart from simple enough group actions, no simple bound for the minimal number of such features required to keep the representation complete is known.Some numerical experiments on the convergence of model performance as a function of the number of intermediate encoded features are shown in Fig. S4.
In parallel, some of the authors in this work are also pursuing an alternative way to use the local coordinate systems defined by the triplets.Similar to the WMIs [42], one makes use of the coordinates of the atoms with respect to this local frame.Instead of using unordered sets however, the coordinates are used more directly, albeit with some modifications, as inputs to a neural network for the final target property, which is discussed in a separate article [43].

IV. DATASET GENERATION AND MODEL DETAILS
To generate a pair of bispectrum-degenerate structures we follow the construction in Ref. 7. We begin by placing an atom at the origin.A ring of radius r with n atoms (one oriented parallel to the x axis, the others having azimuthal angles ϕ i i ∈ {1, . . .n−1}) can then be placed with center (0,0,+z 1 ) and subsequently, the same ring with an azimuthal phase shift of ψ is placed at a height of −z 1 .An additional atom is then placed at coordinates (0,0, ±z 2 ) to generate the A ± structures.From here we can see that the manifold of bispectrum-degenerate structures composed of 2n + 2 atoms is parameterized by {r, ϕ 1 , . . .ϕ n−1 , ψ, z 1 , z 2 }.Additionally, the species of the central atom, the ring of atoms at z = ±z 1 , and the atom at (0,0,±z 2 ) can be different.
The dataset of 4000 pairs of bispectrum-degenerate B 8 structures was generated by fixing n = 3 (yet still being able to randomly position the atoms on the ring) and sampling these parameters randomly within a reasonable range (essentially selecting structures with minimal B -B distances in a range comparable to that found in metastable cluster geometries).The energies for these structures were obtained by PBE-DFT calculations performed using PySCF in the cc-pVDZ basis.As an atom-centered descriptor of correlation order ν = 7 can in principle simultaneously describe all atoms in the structure, this representation is complete (in the limit of a complete basis).In the following, we compare the proposed three-center one-neighbor representation, which is universally complete, with the ACDC feature that is complete for this dataset.As a reminder, by completeness, we mean that the features take on unique values for structures unrelated by symmetric transformations that are already encoded in their construction (here, permutation of identical atoms, rigid translations and (im)proper rotations.
A. Feature computation ACDC features up to ν = 7 and the three-center, one neighbor correlation (N = 2, ν = 1) features were computed for 4000 pairs of degenerate B 8 structures with a cutoff radius of 2.0 Å (that encompasses all atoms in the cluster), Gaussian smearing of 0.2 Å n max = 6, l max = 3.Finally, ACDCs from ν = 1 to up to ν = 7 were stacked together and compressed using PCA so that the final atom-centered descriptor had as many features as the number of structures.

B. Details of the feature reconstruction models
The encoder-decoder architectures for reconstructing the atom-centered feature (as described above, thus fixing the output size to the total number of structures ) from the three-centered description were implemented in PyTorch [59] and differed in the transformations employed in the encoder.Here, we describe two types of models, namely linearly encoded and nonlinearly encoded triplet features for a given latent space size (which we fix to be 128 for all models discussed in the main text).Linearly encoded models employed just a single linear layer to contract the input, whereas a nonlinear encoder sequentially contracted the input to the final latent size over three layers.The first linear layer compressed the dimension to a hidden size of 512 neurons, following which the output was normalized (using LayerNorm) and transformed with a SiLU nonlinearity.The second layer preserved the hidden dimension but also involved normalization and nonlinearity.The final layer condensed the features to the latent space size.
Both the encoders described above were subsequently used with a decoder comprising three linear layers interspersed with a normalization layer (LayerNorm) and SiLU nonlinearities as in the case of the nonlinear encoder, except that the dimensions sequentially increased from 128 to 512 to the output size.
The validation error in reconstructing |ρ ⊗7 i ⟩ features converges quickly with respect to the latent-space size, suggesting that a relatively small number of descriptors is necessary to fully characterize the environment structure in a rotationally-invariant fashion (Fig. S4).

C. Details of energy models
We consider as targets, the difference between the energy of the structure and the one obtained for the minimum energy structure (a centered heptagon).Just as in the case of feature models, we considered linear and nonlinear models.The linear models consist of two linear layers sequentially reducing the input features to the a latent-space size of 128 features and finally to a single scalar output per structure.We use this structure even though obviously a single weight vector would suffice to mimic the behavior of the contracted features model during training.Nonlinear models on the other hand first compress the features to a hidden size of 128 following which a SiLU nonlinearity is applied and the procedure is repeated for another intermediate layer before producing a single energy output.

FIG. S4.
a) Convergence of overall RMSE in the prediction of the |ρ ⊗7 i ⟩ features using encoder/decoder architectures with changing hidden layer dimension (latent space sizes) considering both structures in the degenerate pairs (A + , A − ).b) Highlights, in particular, the convergence of relative RMSE in the difference of the features of the degenerate pairs.

D. Training Errors
We report the errors incurred while training the various models described on a dataset of 3500 pairs of the B 8 clusters in Table S2.For simplicity, we denote by the subscript L if the models employed were linear, or NL if the model was nonlinear.A ρ L , A ρ NL with superscript ρ indicates that the model built on three-center features was targeting the high ν atom-centered feature.The encoded features from A ρ L and A ρ NL models were subsequently used as inputs with the nonlinear energy models, yielding A L , A NL models respectively.If instead the energies were targeted directly by using the N = 2, ν = 1 features, we obtain B L and B NL respectively.We also consider both linear and nonlinear energy models when using |ρ ⊗7 i ⟩ as the inputs, C L and C NL respectively.Additionally, we construct D L and D NL that are linear and nonlinear energy models using |ρ ⊗3 i ⟩ as the input.

FIG. 1 .
FIG. 1. (a)Representations of an atomic environment Ai in terms of neighbor-atom pairs, or symmetrized three-point correlations.The two schemes are equivalent in the complete basis set limit and generalize to higher-order representations.(b) A pair of environments that are degenerate to ν = 1 (listof-distances) representations.

FIG. 3 .
FIG. 3. (a) Correlation plot of the distances between pairs of structures from the B8 dataset, computed based on the bispectrum (d3, red), the intermediate encoded features from Model A ρ NL (dξ, blue), and the full-body features (d7), taken as fully-discriminating descriptors.Configurations are randomly distorted ( Ã+ in the image), to reveal the behavior in the vicinity of the singular points (a bispectrum degenerate pair A + and A − is also shown in the inset), for which the bispectrum distance would be exactly zero.(b) Predictions for a representative 250 pairs of the validation set from model BNL(blue) and DNL (red).Degenerate pairs are joined by a line.The bispectrum-based predictions are identical for the pair, whereas the triplet features ξ can resolve the difference in energy of the pair and achieve better accuracy overall.

TABLE I
. Accuracy of different models (we indicate with 'L'

TABLE S2 .
Accuracy of different models for learning the high-ν features, and the energy of the train set of 3500 bispectrumdegenerate pairs of B8 structures.Similar to the main text, we report separately the error on the mean of the features and energies of degenerate pairs (A + , A − ), and on their difference, normalized by the standard deviation of the train set.Absolute energy errors can be obtained by multiplying the values by 3.84 eV (Mean A ± ) and 0.250 eV (Diff.A ± ).