Representations of molecules and materials for interpolation of quantum-mechanical simulations via machine learning
Marcel F. Langer Alex Goeßmann Matthias Rupp
physics.comp-ph cond-mat.mtrl-sci cs.LG physics.chem-ph







[\begin@twocolumnfalse




    
	Marcel F. Langer[1][2] Alex Goeßmann[1][3] Matthias Rupp[1][4][5]

    
	
	  NOMAD Laboratory, Fritz Haber Institute of the Max Planck Society, Berlin, Germany;    Machine Learning Group, Technische Universität Berlin, Germany;    Institute of Mathematics, Technische Universität Berlin, Germany;    Citrine Informatics, Redwood City, CA, USA    Present address: Department of Computer and Information Science, University of Konstanz, Germany.
 

===========================================================================================================================================================================================================================================================================================================================================================================================





  §.§.§ Abstract

Computational study of molecules and materials from first principles is a cornerstone of physics, chemistry, and materials science, but limited by the cost of accurate and precise simulations.
In settings involving many simulations, machine learning can reduce these costs, often by orders of magnitude, by interpolating between reference simulations.
This requires representations that describe any molecule or material and support interpolation.

We comprehensively review and discuss current representations and relations between them,
using a unified mathematical framework based on many-body functions, group averaging, and tensor products.
For selected state-of-the-art representations, we compare energy predictions for organic molecules, binary alloys, and Al-Ga-In sesquioxides in numerical experiments controlled for data distribution, regression method, and hyper-parameter optimization.



\end@twocolumnfalse]



-4ex
	[main]
	

§ 

	[main]1[2]


+4ex


§ GLOSSARY



	Acronym     Meaning 
 
	2@l@Representations 

	ACE        atomic cluster expansion 

	BoB        bag of bonds 

	BS         bispectrum 

	CM         Coulomb matrix 

	DECAF      density-encoded canonically-aligned fingerprint 

	FCHL       Faber-Christensen-Huang-von Lilienfeld 

	GM         Gaussian moments 

	HDAD       histograms of distances, angles, and dihedral angles 

	IDMBR      inverse-distance many-body representation 

	MBTR       many-body tensor representation 

	MOB        molecular orbital basis 

	MTP        moment tensor potential 

	NICE       N-body iterative contraction of equivariants 

	OMF        overlap matrix fingerprint 

	SF         symmetry function 

	SOAP       smooth overlap of atomic positions 
	
	WST        wavelet scattering transform 

	

	2@l@Methodology 

	GPR        Gaussian process regression 

	HP         hyperparameter (free parameter) 

	KRR        kernel ridge regression 

	MAE        mean absolute error 

	ML         machine learning 

	QM         quantum mechanics 

	QM/ML      ML model for accurate prediction of QM data 

	RMSE       root mean squared error 

	system     poly-atomic system, e.g., a molecule or a crystal 

	

	2@l@Other 

	SI         supplementary information

	








-2



§ INTRODUCTION


Quantitative modeling of atomic-scale phenomena is central for scientific insights and technological innovations in many areas of physics, chemistry, and materials science.
Solving the equations that govern quantum mechanics (QM), such as Schrödinger's or Dirac's equation, allows accurate calculation of the properties of molecules, clusters, bulk crystals, surfaces, and other polyatomic systems.
For this, numerical simulations of the electronic structure of matter are used, with tremendous success in explaining observations and quantitative predictions. 

However, the high computational cost of these ab initio simulations (<ref>) often only allows investigating
from tens of thousands of small systems with a few dozen atoms to a few large systems with thousands of atoms, particularly for periodic structures.
In contrast, the number of possible molecules and materials grows combinatorially with the number of atoms: 
13 or fewer C, N, O, S, Cl atoms can form a billion possible molecules, <cit.>
and for 5-component alloys, there are more than a billion possible compositions when choosing from 30 elements. (<ref>)

This limits systematic computational study and exploration of molecular and materials spaces.
Similar considerations hold for ab initio dynamics simulations, which are typically restricted to systems with a few hundred atoms and sub-nanosecond timescales.

Such situations require many simulations of systems correlated in structure, implying a high degree of redundancy.
Machine learning <cit.> (ML) exploits this redundancy to interpolate between reference simulations <cit.> (<ref>).
This ansatz replaces most ab initio simulations by ML predictions, based on a small set of reference simulations.
Effectively, it maps the problem of repeatedly solving a QM equation for many related systems onto a regression problem.
This approach has been demonstrated in benchmark settings <cit.> and applications <cit.>,
with reported speed-ups between zero to six orders of magnitude. <cit.>
It is currently regarded as a highly promising avenue towards extending the scope of ab initio methods.




The most relevant aspect of ML models for interpolation of QM simulations (QM/ML models) after data quality (<ref>) is the definition of suitable input features, that is, representations of atomistic systems.
Representations define how systems relate to each other for regression and are the subject of this review.



 §.§ Scope and structure



QM/ML models require a space in which interpolation takes place.
Such spaces can be defined explicitly, often as vector spaces, or implicitly, for example, via a kernel function in kernel-based machine learning. <cit.>
This work reviews and compares explicit Hilbert-space representations of finite and periodic polyatomic systems for accurate interpolation of QM observables via ML,
focusing on representations that satisfy the requirements in <ref> and energy predictions.

This excludes features that do not encode all input information, such as atomic numbers and coordinates,
for example, descriptors or fingerprints used in cheminformatics and materials informatics to interpolate between experimental outcomes, <cit.> 
and implicit representations learned by end-to-end deep neural networks <cit.>
or defined via direct kernels between systems. <cit.> (<ref>)

Characteristics and requirements of representations are discussed in <ref>.
<ref> describes a unified mathematical framework for representations.
Specific representations are delineated (benchmarked ones in <ref>, others in <ref>), qualitatively compared (<ref>), and empirically benchmarked (<ref>).
We conclude with an outlook on open problems and possible directions for future research in <ref>.



 §.§ Related work


Studies of QM/ML models often compare their performance estimates with those reported in the literature.
While such comparisons have value, they entertain considerable uncertainty due to different datasets, learning algorithms, including choice of hyperparameters (HPs, free parameters), sampling, validation procedures and reported quantities.
Accurate, reliable performance estimates require a systematic comparison that controls for the above factors, which we perform in this work.

Several recent studies systematically measured and compared prediction errors of representations (<ref>).
We distinguish between studies that automatically (as opposed to manually) optimize
numerical HPs of representations, for example, the width of a normal distribution;
structural HPs of representations, for example, choice of basis functions;
and HPs of the regression method, for example, regularization strength.
<ref> discusses the individual studies in <ref>.













§ ROLE AND TYPES OF REPRESENTATIONS


An N-atom system formally has 3N-6 degrees of freedom.[4]
Covering those with M samples per dimension requires M^3N-6 reference calculations, which is infeasible except for the smallest systems.
How then is it possible to learn high-dimensional energy surfaces?

Part of the answer is that learning the whole energy surface is unnecessary, as configurations high in energy become exponentially unlikely—it is sufficient to learn low-energy regions.
Another reason is that the regression space's formal dimensionality is less important than the data distribution in this space. (<ref>)
Representations can have thousands of dimensions, but their effective dimensionality <cit.> can be much lower if they are highly correlated.

The role of representations is, therefore, to map atomistic systems to spaces amenable to regression. 
These spaces, together with the data's distribution, determine the efficiency of learning.

We distinguish between local representations that describe parts of an atomistic system, such as atoms in their environment, <cit.> and global ones that describe the whole system.
For global representations, represented systems are either finite, such as molecules and clusters, or periodic, such as bulk crystals and surfaces.
(<ref>)



Local representations are directly suitable for local properties, such as forces, nuclear magnetic resonance shifts, or core-level excitations, <cit.> which depend only on a finite-size environment of an atom.
Extensive global properties (<ref>) such as energies can be modeled with local representations via additive approximations, summing over atomic contributions (<ref>).
Since local representations require only finite support, it does not matter whether the surrounding system is finite or periodic.
Global representations are suited for properties of the whole system, such as energy, band gap, or the polarizability tensor.
Since periodic systems are infinitely large, global representations usually need to be designed for or adapted to these.
Trade-offs between local and global representations are discussed in <ref>.

Historically, interpolation has been used to reduce the effort of numerical solutions to quantum problems from the beginning.
Early works employing ML techniques such as Tikhonov regularization and reproducing kernel Hilbert spaces in the late 1980s and throughout the 1990s were limited to small systems. <cit.>
Representations for high-dimensional systems appeared a decade later, <cit.>
underwent rapid development, and constitute an active area of research today.
<Ref> presents an overview.










§ REQUIREMENTS


The figure of merit of ML models for fast, accurate interpolation of ab initio properties is sample efficiency:
The number of reference simulations required to reach a given target accuracy.
Imposing physical constraints on representations improves their sample efficiency by removing the need to learn these constraints from the training data.
The demands of speed, accuracy, and sample efficiency give rise to specific requirements,
some of which depend on the predicted property:


	
  * Invariance to transformations that preserve the predicted property, including 
		
			
  *  changes in atom indexing (input order, permutations of like atoms),
		
		and often 
		[resume]
			
  *  translations
			
  *  rotations 
			
  *  reflections. 
		 
		Predicting tensorial properties requires
		[resume]
			
  *  covariance (equivariance) with rotations.
		 <cit.>
		
		Dependence of the property on a global frame of reference, for example, due to the presence of a non-isotropic external field, can affect variance requirements.
	
	
  * Uniqueness, 
		that is, variance against all transformations that change the predicted property:
		Two systems that differ in property should be mapped to different representations.
		
		Systems with equal representation that differ in property introduce errors: <cit.>
		Because the ML model cannot distinguish them, it predicts the same value for both, resulting in at least one erroneous prediction.
		Uniqueness is necessary and sufficient for reconstruction, up to invariant transformations, of an atomistic system from its representation. <cit.>
		
	
  *  
  *  Continuity,  
		and ideally [resume] 
  *  differentiability, 
		with respect to atomic coordinates.
		
		Discontinuities work against the regularity assumptions in ML models, which try to find the least complex function compatible with the training data.
		Intuitively, continuous functions require less training data than functions with jumps.
		Differentiable representations enable differentiable ML models.
		If available, reference gradients can further constrain the interpolation function (“force matching”), improving sample efficiency. <cit.>
		
	
  * Computational efficiency relative to the reference simulations. 

        For an advantage over simulations alone (without ML), overall computational costs should be reduced by one or more orders of magnitude to justify the effort.
        The difference between running reference simulations and computing representations usually dominates costs. (<ref>)
		
		Therefore, the results of computationally sufficiently cheaper simulations, for example, from a lower level of theory, can be used to construct representations <cit.> or to predict properties at a higher level of theory (“Δ-learning”). <cit.>
		
	
  * Structure
		of representations and the resulting data distribution should be suitable for regression. (<ref>)
		It is useful if feature vectors always have the same length.<cit.>.

		Representations often have a Hilbert space structure, featuring an inner product, completeness, projections, and other advantages.
		Besides the formal space defined by the representation, the structure of the subspace spanned by the data is critical. <cit.>
		This requirement is currently less well understood than <ref>–<ref> and evaluated mostly empirically (see <ref>).
		
	
  * Generality,
		in the sense of being able to encode any atomistic system.

		While current representations handle finite and periodic systems,
		less work was done on charged systems, excited states, continuous spin systems, isotopes, and systems subjected to external fields.


Simplicity, both conceptually and in terms of implementation, is, in our opinion, a desirable quality of representations, albeit hard to quantify.



The above requirements preclude direct use of Cartesian coordinates, which violate requirement <ref>,
and internal coordinates, which satisfy <ref>–<ref> but are still system-specific, violating <ref> and possibly <ref> if not defined uniquely.
Descriptors and fingerprints from cheminformatics <cit.> and materials informatics violate <ref> and <ref>.

Simple representations such as the Coulomb matrix (<ref>) either suffer from coarse-graining, violating <ref>, or from discontinuities, violating <ref>.

In practice, representations do not satisfy all requirements exactly (<ref>) but can achieve high predictive accuracy regardless;
for example, for some datasets, modeling a fraction of higher-order terms can be sufficiently unique already. <cit.>
The optimal interaction orders to utilize in a representation also depend on the type and amount of data available. <cit.>








§ A UNIFIED FRAMEWORK


Based on recent work <cit.> we describe concepts and notation towards a unified treatment of representations.
For this, we successively build up Hilbert spaces of atoms, k-atom tuples, local environments, and global structures,
using group averaging to ensure physical invariants and tensor products to retain desired information and construct invariant features.



 §.§ Representing atoms, environments, systems



Information about a single atom, such as position and proton number, is represented as an abstract ket |α⟩ in a Hilbert space ℋ_α.
Relations between k atoms, where their order can matter, are encoded as k-body functions g_k : ℋ_α^× k→ℋ_g. (<ref>)
These functions can be purely geometric, such as distances or angles, but could also be of (al)chemical or mixed nature.
Tuples of atoms and associated many-body properties are thus elementary tensors of a space ℋ≡ℋ_α^⊗ k⊗ℋ_g,


    |𝒜_α_1...α_k⟩≡|α_1⟩⊗ ... ⊗|α_k⟩⊗ g_k(|α_1⟩,...,|α_k⟩).

A local environment of an atom |α⟩ is represented via the relations to its k-1 neighbors by keeping |α⟩ fixed:


    |𝒜_α⟩≡∑_α_1,…,α_k-1|𝒜_α,α_1,…,α_k-1⟩ .


Weighting functions can reduce the influence of atoms far from |α⟩;
we include these in g_k.

An atomistic systems as a whole is represented by summing over the local environments of all its atoms:


    |𝒜⟩ = ∑_α_i|𝒜_α_i⟩ = ∑_α_1,…,α_k|𝒜_α_1,…,α_k⟩ .


For periodic systems, this sum diverges, which requires either 
exploiting periodicity, for example, by working in reciprocal space, or 
employing strong weighting functions and keeping one index constrained to the unit cell. <cit.>



 §.§ Symmetries, tensor products, and projections


Representations incorporate symmetry constraints (<ref>) by using invariant many-body functions g_k, such as distances or angles, 
or through explicit symmetrization via group averaging <cit.>.
Explicit symmetrization transforms a tensor |T⟩ by integrating over a symmetry group 𝒮 with right-invariant Haar measure dS,


    |T⟩_𝒮≡∫_𝒮 S|T⟩dS ,


where symmetry transformations S∈𝒮 act separately on each subspace of ℋ or parts thereof.
For example, for rotational invariance, only the atomic positions in ℋ_α change.

Rotationally invariant features can be derived from tensor contractions, as any full contraction of contravariant with covariant tensors yields rotationally invariant scalars.<cit.>

Sometimes group averaging can integrate out desired information encoded in |T⟩. 
To counter this, one can perform tensor products of |T⟩ with itself, effectively replacing ℋ by ℋ^⊗ν.
Together, this results in a generalized transform


    |T^ν⟩_𝒮≡∫_𝒮 (S|T⟩)^⊗νdS .


To retain only part of the information in 𝒜, 
one can project onto orthogonal elements {|h_l⟩}_l=1^m of ℋ 
via an associated projection operator 𝒫 = ∑_l |h_l⟩⟨h_l|. 
Inner products and induced distances between representations are then given by


    ⟨𝒜 | 𝒫 | 𝒜'| ⟩and 
    	d_𝒫(|𝒜⟩,|𝒜'⟩) = ||𝒫|𝒜⟩-𝒫|𝒜'⟩||_ℋ .









§ SELECTED REPRESENTATIONS


We discuss three representations that fulfill the requirements in <ref> and for which an implementation not tied to a specific regression algorithm and supporting finite and periodic systems was openly available.
These representations are empirically compared in <ref>.






 §.§ Symmetry functions


Symmetry functions<cit.> (SFs) describe k-body relations between a central atom and the atoms in a local environment around it. (<ref>)
They are typically based on distances (radial SFs, k=2) and angles (angular SFs, k=3).
Each SF encodes a local feature of an atomic environment, for example the number of H atoms at a given distance from a central C atom.

For each SF and k-tuple of chemical elements, contributions are summed.
Sufficient resolution is achieved by varying the HPs of an SF.
For continuity (and differentiability), a cut-off function ensures that SFs decay to zero at the cut-off radius.

Two examples of SFs from b2011eq (see <ref> for further references and SFs) are


    G^2_i = ∑_jexp(-η (d_ij-μ)^2 ) f_c(d_ij)         
    
    	
    	G^4_i =  2^1-ζ∑_j,k ≠ i (1 + λcosθ_ijk)^ζ ·          
    exp(-η (d_ij^2 + d_ik^2 + d_jk^2)) 
    	        f_c(d_ij)   f_c(d_ik)   f_c(d_jk)


where η, μ, ζ, λ are numerical HPs controlling radial broadening, shift, angular resolution, and angular direction, respectively, d_ij is a distance, θ_ijk is the angle between atoms i, j, k, and f_c is a cut-off function.
<Ref> illustrates the radial SFs in <ref>.
The choice of which SFs to use is a structural HP.

Variants of SFs include partial radial distribution functions, <cit.> SFs with improved angular resolution <cit.> and reparametrizations for improved scaling with the number of chemical species <cit.>.

In terms of the unified notation, SFs use invariant functions g_k based on distances and angles, multiplied by a cut-off function, to describe local environments |𝒜_α⟩.
Projections 𝒫 onto tuples of atomic numbers then separate contributions from different combinations of chemical elements.









 §.§ Many-body tensor representation


The global many-body tensor representation <cit.> (MBTR) consists of broadened distributions of k-body terms, arranged by element combination. 
For each k-body function and k-tuple of elements, all corresponding terms (for example, all distances between C and H atoms) are broadened and summed up (<ref>). 
The resulting distributions describe the geometric features of an atomistic system:


    f_k(x, z_1, …, z_k) = ∑_i_1,…,i_k w_k  𝒩(x|g_k) ∏_j=1^k δ_z_j,Z_i_j ,


where w_k is a weighting function that reduces the influence of tuples with atoms far from each other, and g_k is a k-body function; both w_k and g_k depend on atoms i_1, …, i_k.
𝒩(x|μ) denotes a normal distribution with mean μ, evaluated at x.
The product of Kronecker δ-functions restricts to the given element combination z_1,…,z_k.



Periodic systems can be treated by using strong weighting functions and constraining one index to the unit cell.

In practice, <ref> can be discretized.

Structural HPs include the choice of w_k and g_k; numerical HPs include variance of normal distributions.

Requiring one atom in each tuple to be the central atom results in a local variant. <cit.>

In terms of the unified notation, MBTR uses distribution-valued functions g_k, including weighting, with distributions centered on k-body terms such as (inverse) distances or angles.
The outer-product structure of |𝒜⟩ corresponds to the product of δ-functions in <ref>, which selects for specific k-tuples of chemical elements.







 §.§ Smooth overlap of atomic positions


Smooth overlap of atomic positions <cit.> (SOAP) representations expand a central atoms' local neighborhood density, approximated by Gaussian functions located at atom positions, in orthogonal radial and spherical harmonics basis functions (<ref>):


    ρ(r⃗) = ∑_n,l,m c_n,l,m  g_n(r⃗)   Y_l,m(r⃗) ,


where c_n,l,m are expansion coefficients, g_n are radial, and Y_l,m are (angular) spherical harmonics basis functions.

From the coefficients, rotationally invariant quantities can be constructed, such as the power spectrum 


    p_n,n',ℓ = ∑_m c_n,ℓ,m c^*_n',ℓ,m  ,


which is equivalent to a radial and angular distribution function, <cit.> and therefore captures up to three-body interactions.

Numerical HPs are the maximal number of radial and angular basis functions, the broadening width, and the cut-off radius.



An alternative to the power spectrum is the bispectrum <cit.> (BS), a set of invariants that couples multiple angular momentum and radial channels. 
The Spectral Neighbor Analysis Potential (SNAP) includes quadratic terms in the BS components. <cit.>
Extensions of the SOAP framework include recursion relations for faster evaluation <cit.> and alternative radial basis functions g_n, such as third- and higher-order polynomials, <cit.> Gaussian functions, <cit.> and spherical Bessel functions of the first kind. <cit.>

In terms of the unified notation, SOAP uses vector-valued g_k to compute the basis set coefficients in <ref>.
Analytic group-averaging (symmetry integration) then results in invariant features such as the power spectrum (ν=2, <ref>) or bispectrum (ν=3).








§ OTHER REPRESENTATIONS


Many other representations were proposed.


The Coulomb matrix <cit.> (CM) globally describes a system via inverse distances between atoms but does not contain higher-order terms.
It is fast to compute, easy to implement, and in the commonly used sorted version (see footnote reference 25 in rtml2012q), allows reconstruction of an atomistic system via a least-squares problem.
However, its direct use of atomic numbers to encode elements is problematic, and it suffers either from discontinuities in the sorted version or from information loss in the diagonalized version as its eigenspectrum is not unique <cit.>.
A local variant exists. <cit.>


The bag-of-bonds <cit.> (BoB) representation uses the same inverse distance terms as the CM but arranges them by element pair instead of by atom pair.
The “BA-rep­res­ent­a­tion” <cit.> extends this to higher-order interactions by using bags of dressed atoms, distances, angles, and torsions.
The inverse-distance many-body representation <cit.> (IDMBR) employs higher powers of inverse distances and separation by element combinations.


Histograms of distances, angles, and dihedral angles <cit.> (HDAD) are histograms of geometric features organized by element combination.
This global representation is similar to MBTR but typically uses fewer bins, without broadening or explicit weighting. 

The Faber-Christensen-Huang-von Lilienfeld representation <cit.> (FCHL)
describes atomic environments with normal distributions over 
row and column in the periodic table (k=1), interatomic distances (k=2),
and angles (k=3), scaled by power laws. 
In the FCHL18 variant, <cit.> the full continuous distributions are used, requiring an integral kernel for regression.
Among other optimizations, FCHL19 <cit.> discretizes these distributions, 
similar to the approach taken by SFs, and can be used with standard vector kernels.


Wavelet scattering transforms <cit.> (WST) use a convolutional wavelet frame representation to describe variations of (local) atomic density at different scales and orientations.
Integrating non-linear functions of the wavelet coefficients yields invariant features, where second- and higher-order features couple two or more length scales.
Variations use different wavelets (Morlet, <cit.> solid harmonic, or atomic orbital <cit.>) and radial basis functions (exponential <cit.>, Laguerre polynomials <cit.>).


Moment-tensor potentials <cit.> (MTP) describe local atomic environments using a spanning set of efficiently computable, rotationally and permutationally invariant polynomials derived from tensor contractions.
Related representations include
Gaussian moments <cit.> (GM), based on contractions of tensors from (linear combinations of) Gaussian-type atomic orbitals;
the N-body iterative contraction of equivariants (NICE) framework, <cit.> which uses recursion relations to compute higher-order terms efficiently;
and atomic cluster expansion <cit.> (ACE), which employs a basis of isometry- and permutation-invariant polynomials from trigonometric functions and spherical harmonics.


Overlap-matrix fingerprints <cit.> (OMF) and related approaches <cit.> employ the sorted eigenvalues (and derived quantities) of overlap matrices based on Gaussian-type orbitals as representation.
Eigenvalue crossings can cause derivative discontinuities, requiring post-processing <cit.> to ensure continuity.
Using a molecular orbital basis (MOB <cit.> and related approaches <cit.>) adds the cost of computing the basis, for example, localized molecular orbitals via a Hartree-Fock self-consistent field calculation.
Other matrices can be used, such as Fock, Coulomb, and exchange matrices, or even the Hessian, for example, from a computationally cheaper reference method.


Density-encoded canonically-aligned fingerprints <cit.> (DECAF) represent the local density in a canonical, invariant coordinate frame found by solving an optimization problem related to kernel principal component analysis.


Tensor properties require covariance (equivariance).
Proposed solutions include 
local coordinates from eigendecompositions, <cit.> which exhibit discontinuities when eigenvalues cross, 
related local coordinate systems, <cit.>
and internal vectors <cit.> (IV), based on inner products of summed neighbor vectors at different scales,
as well as covariant extensions of SOAP <cit.> and ACE <cit.>.





*2




§ ANALYSIS


We discuss
relationships between specific representations,
to which degree they satisfy the requirements in <ref>,
trade-offs between local and global representations, 
and relationships to other models and modeling techniques, including systematic selection and generation of features.



 §.§ Relationships between representations


All representations in <ref> and most representations in <ref> are related through the concepts in <ref>.
We distinguish two primary strategies to deal with invariances,
the use of invariant k-body functions (BoB, CM, FCHL, HDAD, IDMBR, MBTR, SF)
and explicit symmetrization (ACE, BS, GM, MOB, MTP, NICE, OMF, SOAP, WST).

A similar distinction can be made for kernels. <cit.>
Some representations share specific connections:

An evenly-spaced grid of SFs can be seen as a histogram of distances, angles, or higher-order terms, similar to MBTR and HDAD. 
This suggests a local MBTR or HDAD variant by restricting summation to atomic environments, <cit.> and a global variant of SFs by summing over the whole system.

A difference is that MBTR explicitly broadens k-body terms, whereas SFs implicitly broaden them via the exponential functions in <ref>.

The original formulation of SFs represents each chemically distinct central atom separately,
whereas MBTR represents each (unique) tuple of k elements in separate tensor components.
Both approaches correspond to using Kronecker δ functions on element types.

ACE, BS, GM, MTP, NICE, and SOAP share the idea of generating tensors that are then systematically contracted to obtain rotationally invariant features.
These tensors should form an orthonormal basis, or at least a spanning set, for atomic environments.
Formally, expressing a local neighborhood density in a suitable basis before generating derived features avoids asymptotic scaling with the number of neighboring atoms, <cit.>
although HPs, and thus runtime, still depend on it.
Within a representation, recursive relationships can exist between many-body terms of different orders. <cit.>
d2019bq,d2020q,dbcdeoo2020 discuss technical details of the relationships between ACE and SFs, BS, SNAP, SOAP, MTP.



 §.§ Requirements


Some representations, in particular early ones such as the CM, do not fulfill all requirements in <ref>.
Most representations fulfill some requirements only in the limit,
that is, absent practical constraints such as truncation of infinite sums, short cut-off radii, and restriction to low-order interaction terms.

The degree of fulfillment often depends on HPs,
such as truncation order, the length of a cut-off radius, or the highest interaction order k used.
Effects can be antagonistic;
for example, in <ref>, both <ref> uniqueness and <ref> computational effort increase with n,l,m. <cit.>

In addition, not all invariances of a property might be known or require additional effort to model, for example, symmetries <cit.>.

Mathematical proof or systematic empirical verification that a representation satisfies a requirement or related property are sometimes provided:

The symmetrized invariant moment polynomials of MTPs form a spanning set for all permutationally and rotationally invariant polynomials; <cit.> basis sets can also be constructed. <cit.>
For SOAP, systematic reconstruction experiments demonstrate the dependence of uniqueness on parametrization. <cit.>

While <ref> uniqueness guarantees that reconstruction of a system up to invariances is possible in principle, accuracy and complexity of this task vary with representation and parametrization.
For example, reconstruction is a simple least-squares problem for the global CM as it comprises the whole distance matrix D_ij = ||r⃗_⃗i⃗-r⃗_⃗j⃗||_2,
whereas for local representations, (global) reconstruction is more involved.

If a local representation comprises only up to 4-body terms then there are degenerate environments that it cannot distinguish, <cit.> but that can differ in property.
Combining representations of different environments in a system can break the degeneracy.
However, by distorting feature space <ref> structure, these degeneracies degrade learning efficiency and limit achievable prediction errors, even if the training set contains no degenerate systems. <cit.>
It is currently unknown whether degenerate environments exist for representations with terms of order k>4.
The degree to which a representation is unique can be numerically investigated through the eigendecomposition of a sensitivity matrix based on a representation's derivatives with respect to atom coordinates. <cit.>



 §.§ Global versus local representations


Local representations can be used to model global properties by assuming that these decompose into atomic contributions.
In terms of prediction errors, this tends to work well for energies. (<ref>)

Learning with atomic contributions adds technical complexity to the regression model and is equivalent to pairwise-sum kernels on whole systems, (<ref>)
with favorable computational scaling for large systems (see <ref>).
Other approaches to creating global kernels from local ones exist. <cit.>

Conversely, using global representations for local properties can require modifying the representation to incorporate locality and directionality of the property. <cit.> 

A general recipe for constructing local representations from global ones is to require interactions to include the central atom, starting from k=2. <cit.>



 §.§ Relationships to other models and techniques


Two modeling aspects directly related to representations are which subset of the features to use and the construction of derived features.
Both modulate feature space dimensionality and <ref> structure.

Adding products of 2-body and 3-body terms as features, for example,
can improve performance, <cit.> as these features relate to higher-order terms, (<ref>)
but can also degrade performance if the features are unrelated to the predicted property,
or if there is insufficient data to infer the relationship.
Feature selection tailors a representation to a dataset by selecting a small subset of features that still predict the target property accurately enough.
Optimal choices of features depend on the data's size and distribution.

In this work, we focus exclusively on representations.
In kernel regression, however, kernels can be defined directly between two systems, without an explicit intermediate representation.
For example, n-body kernels between atomic environments can be systematically constructed from a non-invariant Gaussian kernel using Haar integration, or using invariant k-body functions (<ref>), yielding kernels of varying body-order and degrees of freedom. <cit.>

Similarly, while neural networks can use representations as inputs, their architecture can also be designed to learn implicit representations from the raw data (end-to-end learning).

In all cases, the Requirements in <ref> apply.








§ EMPIRICAL COMPARISON


We benchmark prediction errors for the representations from <ref> on three benchmark datasets.
Since our focus is exclusively on the representations, we control for other factors, in particular for 
data distribution, regression method, and HP optimization.



 §.§ Data


The  consensus benchmarking dataset<cit.> comprises 133 885 organic molecules composed of H, C, N, O, F with up to 9 non-H atoms. (<ref>)
Ground state geometries and properties are given at the DFT/B3LYP/6-31G(2df,p) level of theory.
We predict U_0, the atomization energy at 0 K.

The  dataset <cit.> (<ref>)
contains the ten binary alloys AgCu, AlFe, AlMg, AlNi, AlTi, CoNi, CuFe, CuNi, FeV, and NbNi.
For each alloy system, it comprises all structures with up to 8 atoms for face-centered cubic (FCC), body-centered cubic (BCC), and hexagonal close-packed (HCP) crystal types, 15 950 structures in total.
Formation energies of unrelaxed structures are given at the DFT/PBE level of theory.

The  challenge <cit.> dataset <cit.> (<ref>) contains 3 000 ternary (Al_x-Ga_y-In_z)_2O_3 oxides, x+y+z=1, of potential interest as transparent conducting oxides.
Formation and band-gap energies of relaxed structures are provided at the DFT/PBE level of theory.
The dataset contains both relaxed (, used here) and approximate () structures as input.
In the challenge, energies of relaxed structures were predicted from approximate structures.

Together, these datasets cover finite and periodic systems, organic and inorganic chemistry, and ground state as well as off-equilibrium structures.
See <ref> for details.




 §.§ Method


We estimate prediction errors as a function of training set size (learning curves, <ref>).
To ensure that subsets are representative, we control for distribution of elemental composition, size, and energy. (<ref>)
This reduces the variance of performance estimates and ensures the validity of the independent-and-identically-distributed data assumption inherent in ML.

All predictions are on data never seen during training.

We use kernel ridge regression <cit.> (KRR; predictions are equivalent to those of Gaussian process regression, <cit.> GPR) with a Gaussian kernel as an ML model. (<ref>)
KRR is a widely-used non-parametric non-linear regression method.
There are two regression HPs, the length scale of the Gaussian kernel and the amount of regularization. (<ref>)

In this work, training is exclusively on energies; in particular, derivatives are not used.
All HPs, that is, regression HPs, numerical HPs (e.g., a weight in a weighting function), and structural HPs (e.g., which weighting function to use), are optimized with a consistent and fully automatic scheme based on sequential model-based optimization and tree-structured Parzen estimators. <cit.> (<ref>)

This setup treats all representations on equal footing.
See <ref> for details on the optimized HPs.




 §.§ Results


<Ref> presents learning curves for SF, MBTR, SOAP on datasets , ,  (see <ref> for tabulated values).
For each dataset, representation, and training set size, we trained a KRR model and evaluated its predictions on a separate hold-out validation set of size 10 k (), 1 k (), and 0.6 k ().
This procedure was repeated 10 times to estimate the variance of these experiments.



Boxes, whiskers, horizontal bars, and crosses show interquartile ranges, minimum / maximum value, median, and mean, respectively, of the root mean squared error (RMSE) of hold-out-set predictions over repetitions.

We show RMSE as it is the loss minimized by least-squares regression such as KRR, and thus a natural choice.
For other loss functions, see <ref>.

From statistical learning theory, RMSE decays as a negative power of training set size (a reason why learning curves are preferably shown on log-log plots). <cit.>
Lines show corresponding fits of mean RMSE, weighted by the standard deviation for each training set size.

<Ref> reveals dependencies between the time to compute representations for a training set (horizontal axis) and RMSE (vertical axis).
When comparing observations in two dimensions, here time t and error e, there is no unique ordering <,
and we resort to the usual notion of dominance:
Let x⃗, x⃗'⃗∈ℝ^d;
then x⃗ dominates x⃗'⃗ if x_i ≤ x'_i for all dimensions i and x_i < x'_i for some i.
The set of all non-dominated points is called the Pareto frontier,
shown by a line, with numbers indicating training set sizes.
<Ref> presents compute times for representations (<ref> for kernel matrices).


-



 §.§ Discussion


Asymptotically, observed prediction errors for all representations on all datasets relate as


    SF-2,3   ≺SF-2 ,        MBTR-2,3   ≼MBTR-2 ,   
    SOAP   ≺SF-2,3 ,    SOAP   ≺MBTR-2,3 , 
    SF-2,3   ≼MBTR-2,3 ,    SF-2   ≺MBTR-2 ,


where A ≺ B (A ≼ B) indicates that A has lower (or equal) estimated error than B asymptotically.
Except for MBTR-2,3⋠SF-2 on dataset ,


    SOAP≺SF-2,3≼MBTR-2,3≺SF-2≺MBTR-2 .


We conclude that, for energy predictions, accuracy improves with modeled interaction order and for local representations over global ones.
The magnitude of, and between, these effects varies across datasets.

Dependence of predictive accuracy on interaction order has been observed by others <cit.>

and might be partially due to a higher resolution of structural features. <cit.>
The latter would only show for sufficient training data, such as for dataset  in <ref>.
We do not observe this for dataset , possibly because angular terms might be immediately relevant for characterizing organic molecules' carbon scaffolds. <cit.>

Better performance of local representations might be due to higher resolution and better generalization 
(both from representing only a small part of the whole structure), and has also been observed by others <cit.>.
The impact of assuming additivity is unclear but likely depends on the structure of the modeled property. (<ref>)
Our comparison includes only a single global representation (MBTR), warranting further study of the locality aspect.

For additional analysis details, see <ref>.

Computational costs tend to increase with predictive accuracy.
Representations should therefore be selected based on a target accuracy, constrained by available computing resources.

Converged prediction errors are in reasonable agreement with the literature (<ref>)
considering the lack of standardized conditions such as sampling, regression method, HP optimization, and reported performance statistics.

In absolute terms, prediction errors of models trained on 10 k samples are closer to the differences between DFT codes than the (systematic) differences between the underlying DFT reference and experimental measurements. (<ref>)










§ CONCLUSIONS AND OUTLOOK



We review representations of atomistic systems, such as molecules and crystalline materials, for machine-learning of ab initio quantum-mechanical simulations.
For this, we distinguish between local and global representations and between using invariant k-body functions and explicit symmetrization to deal with invariances.

Despite their apparent diversity, many representations can be formulated in a single mathematical framework based on k-atom terms, symmetrization, and tensor products.

Empirically, we observe that when controlling for other factors, including distribution of training and validation data, regression method, and HP optimization, both prediction errors and compute time of SFs, MBTR and SOAP improve with interaction order, and for local representations over global ones.



Our findings suggest the following guidance:


	
  * If their prediction errors are sufficient for an application, we recommend two-body versions of simple representations such as SF and MBTR as they are fastest to compute. 
	
  * For large systems, local representations should be used.
	
  * For strong noise or bias on input structures, as in dataset , performance differences between representations vanish, (<ref>)
		and computationally cheaper features that do not satisfy the requirements in <ref> (descriptors) suffice.


We conclude by providing related current research directions, grouped by topic.



Directly related to representations:


	
  * Systematic development of representations via extending the mathematical framework (<ref>) to include more state-of-the-art representations.
		  This would enable deriving “missing” variants of representations (see <ref>), such as a global SOAP <cit.> and local MBTR <cit.>, on a principled basis,
		  as well as understanding and reformulating existing representations in a joint framework, perhaps to the extent of an efficient general implementation. <cit.>

	
  * Representing more systems.
		  Develop or extend representations for atomistic systems currently not representable, or only to a limited extent, such as
		  charged atoms and systems, <cit.> excited states <cit.>, spin systems, isotopes, and systems in an applied external field <cit.>.

	
  * Alchemical learning.
		  Further understand and develop alchemical representations <cit.> that incorporate similarity between chemical species to improve sample efficiency. 
		  What are the salient features of chemical elements that need to be considered, also with respect to charges, excitations, spins, and isotopes?

	
  * Analysis of representations to better understand structure and data distribution in feature spaces and how they relate to physics and chemistry concepts.
		  Possible approaches include
		  quantitative measures of structure and distribution of datasets in these spaces,
		  dimensionality reduction methods, 
		  analysis of data-driven representations from deep neural networks,
		  and construction, or proof of non-existence, of non-distinguishable environments for representations employing terms of order higher than four.

	
  * Explicit complexity control.
		  Different applications require different trade-offs between computational cost and predictive accuracy.
		  This requires determination, and automatic adaptation as an HP, of the capacity (complexity, dimensionality) and computational cost of a representation to a dataset, for example, through selection, combination, <cit.> or systematic construction of features <cit.>.


Related to benchmarking of representations:


	
  * Extended scope.
		  We empirically compare one global and two local representations on three datasets to predict energies using KRR with a Gaussian kernel.
		  For a more systematic coverage, other representations (<ref>) and datasets, training with forces, <cit.> and more properties should be included
		  while maintaining control over regression method, data distribution, and HP optimization.
		  Deep neural networks <cit.> could be included via representation learning.
		  Comparison with simple baseline models such as k-nearest neighbors <cit.> would be desirable.
		  
	
  * Improved optimization of HPs:
		  The stochastic optimizer used in this work required multiple restarts in practice to avoid sub-optimal results, and reached its limits for large HP search spaces.
		  It would be desirable to reduce the influence and computational cost of HP optimization.
		  Possible means include reducing the number of HPs in representations, employing more systematic and thus more robust optimization methods, and providing reliable heuristics for HP default values.
		  
	
  * Multi-objective optimization.
		  We optimize HPs for predictive accuracy on a single property.
		  In practice, though, parametrizations of similar accuracy but lower computational cost would be preferable, and more than one property can be of interest.
		  HPs should, therefore, be optimized for multiple properties and criteria, including computational cost and predictive uncertainties (see below).
		  How to balance these is part of the problem. <cit.>
		  
	
  * Predictive uncertainties.
		  While prediction errors are frequently analyzed, and reasonable guidelines exist, this is not the case for predictive uncertainties.
	      These are becoming increasingly important as applications of ML mature, for example, for human assessment and decisions, learning on the fly, <cit.> and active learning.
	      
	      Beyond global analysis of uncertainty estimates, local characterization (in input or feature space) of prediction errors is relevant. <cit.>


Related through context:


	
  * Long-range interactions.
		  ML models appear to be well-suited for short- and medium-ranged interactions, but problematic for long-ranged interactions due to the increasing degrees of freedom of larger systems and larger necessary cut-off radii of atomic environments.
		  Two approaches are to integrate ML models with physical models for long-range interactions, <cit.> and to adapt ML models to learn long-range interactions directly. <cit.>

	
  * Relationships between QM and ML.
		  A deeper understanding of the relationships between QM and kernel-based ML could lead to insights and technical progress in both fields.
		  As both share concepts from linear algebra, such relationships could be formal mathematical ones.
		  For example, QM concepts such as matrix product states can parameterize non-linear kernel models. <cit.>









 §.§ Acknowledgments


This work received funding from the European Union's Horizon 2020 Research and Innovation Programme, Grant Agreements No. 676580, the NOMAD Laboratory CoE, and No. 740233, ERC: TEC1P.
Part of the research was performed while the authors visited the Institute for Pure and Applied Mathematics (IPAM), which is supported by the National Science Foundation (Grant No. DMS-1440415).

The authors thank Profs. Matthias Scheffler, 
Klaus-Robert Müller,
Jörg Behler,
Gábor Csányi,
O. Anatole von Lilienfeld,
Carsten Baldauf,
Matthew Hirn,

as well as
Emre Ahmetcik,
Lauri Himanen,
Yair Litman,
Dmitrii Maksimov,
Felix Mocanu,
Wiktor Pronobis,
and Christopher Sutton
for constructive discussions.



 §.§ Author contributions


M.F.L. and M.R. designed numerical experiments and analyzed results.
M.F.L. developed software and conducted numerical experiments.
All authors contributed to writing, with emphasis by A.G. on the mathematical framework and M.F.L. on representations and benchmarking.
M.R. supervised the study.



 §.§ Competing interests


The authors declare no competing financial or non-financial interests.



 §.§ Data and code availability


The data that support the findings of this study are publicly available,
datasets at <https://qmml.org>, 
hyperparameter search spaces, machine-learning models, program code, and results at <https://marcel.science/repbench>.

A tutorial introduction to the  Python framework developed for this work is part of the NOMAD Analytics Toolkit. <cit.>






999 
tocsectionReferences



 









§ INVARIANT K-BODY FUNCTIONS

Symmetry invariances can already be implemented by certain choices of k-body functions, which shall be discussed in this chapter oriented on two representation schemes. 
After projections of the representing tensors onto the invariant atom spaces, the group averaging (<ref>) leaves the tensor invariant.




 §.§ Geometric basis functions and moment tensors

Geometric k-body functions g_k are a widely used class of functions depending only on the position vector r_i of the atoms |α_i⟩. For the following we thus set ℋ_α=ℝ^3 and ℋ_g=ℝ. Position vectors depend however on the choice of the origin, which does not influence the properties of a crystal and should therefore not influence the k-body function. Our first assumption on g_k is thus the invariance under translations t∈ℝ^3

    g_k(r_1,...,r_k)    =g_k(r_1+t,...,r_k+t)  
        = g_k(r_1-r_k,...,r_k-1-r_k,0)
       ≡g̃_k-1(u_1,...,u_k-1)

With this assumption the k-body function is described by an effective (k-1)-body function g̃_k dependent on varianbles u_i≡ r_i-r_k. A further group, which should leave the functions invariant, is the orthogonal group O(3) transforming the chosen coordinate frame. We thus assume:

    g̃_k-1(u_1,...,u_k-1)=g̃_k-1(Ru_1,...,Ru_k-1)   ∀ R ∈ O(3)
 
If g̃_k-1 is polynomial dependent on the arguments u_i, the First Fundamental Theorem of the orthogonal group can be applied:

g̃_k can be represented as a polynomial h of n(n+1)/2 scalar variables v_ij≡⟨u_i|u_j|$⟩, where1≤i ≤j ≤k-1.

The monomial basis for polynomials of the variablesv_ijare indexed by symmetric matricesα∈ℕ_0^(k-1)×(k-1):

    ℳ_α(u_1,...,u_k-1)≡∏_1≤ i ≤ j ≤ k-1⟨u_i|u_j|^⟩α_ij

The contraction pattern of the variablesu_ileads to an graphical interpretation of the monomialsℳ_α. To understand this structure further, we define moment tensorsg^μνas2-body functions valued inℝ^ν:

    g^μν(r,r̃)=||r-r̃||^2μ· (r-r̃)^⊗ν

Withr=r_iandr̃=r_k, the moment tensors can be understood as a tensors withνlegs representing the argumentu_i. Scalar productsv_ijbetween argumentsu_iandu_jcan be visualized by a contraction of different tensorial legs. To reproduce eachℳ_α, a collection{g^μ_i ν_i(r_i,r_k)}of moment tensors withμ_i=α_iiandν_i=∑_l<iα_liis fully contracted. The integerα_ijrepresents the number of contracted legs of theith with thejth moment tensor.



Invariant geometrick-body functions have been widely used to represent chemical structures. In the following we will briefly review two representation schemes and their formulation in our framework. 



 §.§ Moment tensor potentials

Moment tensor potentials use the monomial basis decomposition (<ref>) to optimizek-body functions for a particular learning problem, where permutation invariance is demanded in addition to the above. Drawing on their flexibility,k-body functions are used for model building by sums of argument tuples. With the above discussion, geometrick-body functions are parametrized by:

    g_k^{λ_α}(r_1,...,r_k)≡∑_αλ_αℳ_α(r_1-r_k,...,r_k-1-r_k)

For each atom|α_k⟩one takes(k-1)nearest neighbors and sums thek-body contributions for all numberings(α_1,...,α_k-1)of the atoms for the representation of a local environment:

    |𝒜_α_k⟩=∑_α_1,...,α_k-1^dis. neighbors of α_k g_k^{λ_α}(r_1,...,r_k)

Distinguished from the representation scheme (<ref>), we understand (<ref>) as the projection of the underlying tensor representation onto its first dimension. A representation|𝒜⟩of a global structure is achieved by sums of the local terms|𝒜_α_k⟩:

    |𝒜⟩=∑_α_k∑_α_1,...,α_k-1^dis. neighbors of α_k g_k^{λ_α}(r_1,...,r_k)

Because of the omission of the further tensor dimensions in the scheme (<ref>) the representation is a real function on the space of materials and optimization of the parameterλ_αin a supervised setup leads to learning a model. We note that learning instead a linear model on an analog representation (<ref>) would establish alchemical resolution of the model.




 §.§ Many body tensor representation

Another representation scheme based on geometrick-body function is the family of Many body tensor representation, which use all dimensions of a global representation (<ref>) to define distances between structures. A key idea is the restructuring of the spaceℋ_gas a functional spaceL^2(ℝ)instead ofℝ.k-body functionsg̃_kare therefore modelled as Gaussian distributions with their mean given by an associated geometric functiong_k:

    g̃_̃k̃(r_1,...,r_k)[r]=1/C(σ)exp[-||r-g_k(r_1,...,r_k)||/2σ^2]

With the Euclidean scalar product onL^2(ℝ)dissimilarities are measured by mismatch integrals of differently centered distributions. With this trick, the dissimilarity measure provides a sensitivity to the distance between the distribution centers.  
[main]








Supplementary material for

    
	Marcel F. Langer[1][2] Alex Goeßmann[1][3] Matthias Rupp[1][4][5]

    
	
	  NOMAD Laboratory, Fritz Haber Institute of the Max Planck Society, Berlin, Germany;    Machine Learning Group, Technische Universität Berlin, Germany;    Institute of Mathematics, Technische Universität Berlin, Germany;    Citrine Informatics, Redwood City, CA, USA    Present address: Department of Computer and Information Science, University of Konstanz, Germany.
 

===========================================================================================================================================================================================================================================================================================================================================================================================

subsectionsectionseccntformat#1ifundefined#1@cntformat the#1   #1@cntformatcntformatsubsection|[si]


§ CONTENTS
section[0em]
	
[]0emsubsection[2em]1.5em1pc[si]1[2]

§ INTRODUCTION


 §.§ Cost of electronic structure calculations

Although the computational cost of ab initio methods scales only polynomially in system sizeN(measured, for example, in number of electrons or orbitals), it remains a strongly limiting factor. For example, the currently most-widely used approach, Kohn-Sham density functional theory, scales asO(N^3)for (semi)local andO(N^4)for hybrid functionals: 
DoublingNthus increases compute time by roughly an order of magnitude, and a few such doublings will exhaust any computational resource. 
Advances in large-scale computing facilities, such as current exascale computing inititatives, will move this “computational wall” to larger systems, but cannot remove it. 
In practice, the large prefactor hidden in the asymptotic runtimes is relevant as well.



 §.§ Size of molecular and materials spaces

Various estimates of the size of chemical compound spaces exist, popular ones pmv2013x,bmg1996x including10^33and10^60molecules. 
Reymond et al.fr2007x,br2009x,rdbr2012x systematically enumerate all small organic molecules with up to 11 C, N, O, F atoms, 13 C, N, O, S, Cl atoms, and 17 C, N, O, S, F, Cl, Br, I atoms, yielding 26 million, 970 million and 166 billion molecules, respectively. 

Following Cantor, c2014cx we estimate the number of possible compositions (not considering unit cell size or symmetry) for an alloy system to be 
the multinomial coefficient(n-1,k)! = n-1+kk, wherenis number of components andk = 100/xis determined by the tolerancex %to which the amount of a species is specified.
Forn = 5and a very conservative choice ofx=5 %, removing combinations that contain only 4 or fewer components and multiplying by all ways to choose 5 out of 30 elements yields( 5-1+2020 - 4-1+2020 ) 305 ≈1.5·10^9. 



 §.§ The role of data quality for QM/ML models

Data are the basis for data-driven models, and errors in them can only be corrected to a limited extent (“garbage in, garbage out”).
Even dealing with simple errors like independent identically distributed noise requires additional data, and more severe errors lead to qualitative problems such as outliers.
Conversely, problems in fitting a ML model can be indicative of problems in the data.



 §.§ Explicit and implicit features

Features used for regression can be defined explicitly via representations, or implicitly, for example, via kernels or deep neural networks.

In this work, we focus on explicit Hilbert-space representations in conjunction with kernel-based regression with a Gaussian kernel.
Technically, the features used for regression are the components of the kernel feature space, that is, the non-linear transformations of the representations' components via the Gaussian kernel.
While used implicitly in this sense, the representations are still defined explicitly. 

This is in contrast to implicitly defined representations, for example, feature spaces of kernels defined directly on “raw inputs“ such as atomic coordinates and numbers, without an intermediate explicit Hilbert-space representation,
or, the layers of deep neural networks (end-to-end learning).
For the latter, the requirements in <ref> can be imposed via the network architecture,
which can be seen as the conceptual analog to explicitly conformant representations or kernels.



 §.§ Related work


Faber et al.fhhgsdvkrv2017x compare combinations of representations and regression methods for atomization energies of organic molecules ( dataset, see <ref>).
Only some of the tested features are representations that satisfy the requirements in <ref>; their HPs are not optimized.



Himanen et al.hjmfrgrf2020qx

investigate the representations in <ref>, also using kernel regression, 
to predict ionic charges of molecules from the  dataset, as well as formation energies in a custom dataset of inorganic crystals obtained from the Open Quantum Materials Database. 
They optimize numerical HPs of representations and regression method, but not structural ones.

Zuo et al.zcldcbcstwo2020x focus on dynamics simulations, and therefore include forces and stresses in training and evaluation.
They also evaluate predictions of derived physical quantities, such as elastic constants or equation-of-state curves.
Different combinations of representation, regression method, and HP tuning 
are evaluated on a dataset of elemental solids. 
Timings are discussed.

Schmitz et al.sgc2019x compare regression methods for potential energy surfaces of 15 small organic molecules, using non-redundant internal coordinates as features. 
HPs of the representation are not optimized.



Nyshadham et al.nrbsmrcwh2019x compare selected combinations of representations and regression methods on binary alloys ( dataset, see <ref>).
HPs of the representations are not optimized.


Stuke et al.strkghr2019qx evaluate prediction of molecular orbital energies with kernel regression on three datasets: 
organic molecules ( dataset, see <ref>), amino acids and dipeptides, as well as opto-electronically active molecules. 
Numerical HPs of representations and regression method are optimised via local grid search.

Onat et al.ook2020x empirically investigate effective dimensionality and sensitivity (to small perturbations of the underlying system) of representations using both materials and molecular datasets.
HPs are not optimized.
They do not investigate prediction errors.

Parsaeifard et al.pdcfkdblg2020x also study the sensitivity of representations to infinitesimal geometric perturbations (uniqueness), as well as correlations between the distances induced by different representations and with physical properties.

Käser et al.kkclm2020qx empirically evaluate representations with kernel-based regression and a deep neural network for prediction of energies, forces, vibrational modes and infrared spectra of formaldehyde at different levels of theory.

Jäger et al.jmfhf2018x compare SOAP, MBTR, CM and SFs with KRR for prediction of the adsorption free energy of hydrogen on the surface of nanoclusters.
Numerical HPs of KRR, and some numerical HPs of representations are optimized.

Goscinski et al.gfc2020x develop metrics to compare the feature spaces generated by different representations (SFs, SOAP, NICE), exploring the impact of HP choices in datasets of methane and solid carbon.





§ ROLE AND TYPES OF REPRESENTATIONS


 §.§ Structure and distribution of data
<Ref> illustrates the importance of representation space structure for regression with a toy example.
Low-dimensional (here, essentially one-dimensional) data is embedded into a high-dimensional (here, two-dimensional) space.
The spiral embedding is not suited for linear regression, whereas the linear embedding is.



 §.§ Extensive and intensive properties

A property whose magnitude is additive in the size (extent or mass) of an object is called extensive; 
a property whose magnitude is independent of the size of an object is called intensive.
For example, internal energy is an extensive property, band gap energy an intensive one.

Originating from thermodynamics, t1917x,hk1965x the application of these terms to microscopic quantities is limited by allowed changes in “size“ of a system:
For finite systems such as molecules, a propertypis extensive if for any two non-interacting systemsAandB,p(A+B) = p(A) + p(B), jskohrm2020qx
and intensive ifp(A) = p(A+A).
For periodic systems such as bulk crystals, we takeAandBto be supercells of the same unit cell.
In this minimal sense, total and atomization energy of atomistic systems are extensive.

However, energies are not additive for general changes in a system, such as changes in atomic position, and addition or removal of atoms.
With respect to the requirements in <ref>,
ML models for energies should be size-extensive in the (minimal) sense above.

For global representations, this can be achieved via normalization in conjunction with the linear kernel, jskohrm2020qx
whereas local representations as described in <ref> automatically satisfy this requirement.



 §.§ Learning with atomic contributions

One ansatz to scale prediction of global properties to large atomistic systems is to predict atomic contributions. 
This assumes additivity, as the predicted property is a sum of predicted atomic contributions, and locality, as efficient scaling requires representations of atoms in their environment to have local support, often achieved through a finite-radius cut-off function.

Predicting atomic contributions requires a modification of the basic kernel regression scheme, 
which we derive here building on References bc2015x and m2015x:

Let𝒜_1,…,𝒜_adenote atoms of systemsℳ_1,…,ℳ_mand letD ∈{0,1}^m ×abe their incidence matrix, that isD_i,j = 1if𝒜_jbelongs toℳ_iand 0 otherwise.
Letk̃denote a kernel function on atoms.
The prediction for the  system is the sum of its predicted atomic contributions,

    f(ℳ_i) = ∑_j=1^a f(𝒜_j) D_i,j = ∑_j,ℓ=1^aα̃_ℓk̃(𝒜_ℓ,𝒜_j) D_i,j .

Minimizing quadratic loss yields


    _α⃗̃⃗∈^a∑_i=1^m( y_i - ∑_j,ℓ=1^aα̃_ℓk̃(𝒜_ℓ, 𝒜_j) D_i,j)^2 + λ ||f||^2_ℋ
    
    	= _α⃗̃⃗∈^ay⃗ - DK̃α⃗̃⃗ | y⃗ - DK̃α⃗̃⃗ + λ α⃗̃⃗^T K̃α⃗̃⃗ .


Since this is a quadratic form, it suffices to set its gradient to zero and solve forα⃗̃⃗:

    ∇_α⃗̃⃗ = - 2 α⃗̃⃗^T K̃D^T y⃗ + 2 K̃D^T DK̃α⃗̃⃗ + 2 λK̃α⃗̃⃗ = 0⃗
    ⇔ α⃗̃⃗ = ( D^T DK̃ + λI )^-1D^T y⃗
    ⇔ α⃗̃⃗ = D^T ( DK̃D^T + λI )^-1y⃗ ,

where the last expression is preferable for numerical evaluation.
Predictions form'new systemsℳ'witha'atoms𝒜'can be expressed efficiently as


    y⃗'⃗= ( ∑_j=1^a' D'_i,j∑_ℓ=1^a α̃_ℓk̃(𝒜_ℓ, 𝒜'_j) )_i=1,…,m'= D' L̃^T α⃗̃⃗ ,


whereD'is the incidence matrix for the predicted systems andL̃is thea ×a'kernel matrix between atoms𝒜and𝒜'.



This approach is equivalent to kernel regression (<ref>) on whole systems with a kernelkgiven by the sum of the atom kernelk̃over all pairs of atoms in two systems,


    k(ℳ_i,ℳ_j) = ∑_p,q=1^a D_i,pK̃_p,qD^T_q,j .


This follows fromK = D K̃ D^T,L = D L̃ D'^Tandα⃗̃⃗ = D^T α⃗(<ref>):
Predictions for whole systems usingkare identical to predictions usingk̃:y⃗'⃗ = L^T α⃗ = D' L̃^T D^T α⃗ = D' L̃ α⃗̃⃗.

In particular, atomic weightsα⃗̃⃗are blocks of system weightsα⃗.

Computing full atom kernel matricesK̃and incidence matricesDcan require large amounts of memory.
In practice, we compute blocks ofK̃on the fly and directly sum over its entries.
Learning with atomic contributions is extensive (<ref>).





§ REQUIREMENTS


 §.§ Computational cost

Letf^QMandf^MLdenote the total computational cost when using only ab initio simulations and a ML-augmented model, respectively:


    f^QM(n,m)     = (n+m)  ref
    
    	f^ML(n,m)     = n  ref + (n+m)  repr + train(n) + m  pred,


wherenandmare number of training and predicted systems,ref,repr,predare the cost of one simulation, representation calculation and prediction,
andtrainis the cost of training the ML model.
Ifn ≈ϵ  mfor smallϵ, and costs of training and prediction are negligible, total savings in compute time are


    f^QM(n,m) - f^ML (n,m) ≈ m   (ref - repr) .


Bothrefandreprdepend on system size, often polynomially, with differences in asymptotic runtime as well as constant factors relevant in practice.
Local representations require computing more kernel matrix entries (<ref>) than global representations, which can noticeably influence compute time (<ref>), but enable scaling with system size:

Letcdenote the (average) number of atoms per system, anddthe (average) number of atoms in a local environment.
In the following, we assumedto be constant (bounded from above), and representations to have constant size.

Total computational effort to compute representations (first term) and kernel matrices (second term) is then given by


    𝒪( (n+m)   c^k + n m ) 
    	  and  𝒪( (n+m)   c   d^k + n m c^2 )


for global (left) and local (right) representations, whered^kis constant.
For small systemsc ≈d, and the additional overhead in computing kernel matrices will dominate runtime for smallk. 
In the limitc →∞of large systems, thec^kterm will dominate for global representations, while local representations enjoy quadratic scaling.

This can be observed to some extent in <ref>.

The above analysis applies to representations that explicitly computek-body terms.
Representations that computek-body terms as contractions of basis set coefficients enjoy better scaling in the number of atoms. d2019bqx
Recursion relations betweenk-body terms can improve scaling as well. c2019dqx,npc2020qx,dbcdeoo2020x

 §.§ Role of representations

The role of the representation is to map atomistic systems into a space amenable to regression (linear interpolation).
Strictly speaking, for kernel regression this is the kernel feature space, that is, representation space transformed by the kernel.
We limit our discussion to the representation itself—for the linear kernel this is exact as the transformation is the identity,
and many non-linear kernels like the Gaussian kernel act on the representation space, relying on its structure and implied metric.





§ A UNIFIED FRAMEWORK


 §.§ k-body functions

Informally, ak-body function maps information aboutkdistinct atoms|α_1⟩, …, |α_k⟩, where order can matter, to an output space, here the real numbers, or a distribution on them.
Atom information|α⟩typically includes coordinates and proton number, but is not limited to those; for example, it could include neutron number to model isotopes.

Typicalk-body functions include 
atomic number counts (k=1),
distances, sometimes inverted or squared (k=2),
angles or their cosine (k=3),
dihedral or torsional angles, volume-related terms (k=4). 
Less common, (al)chemical relationships can be included, for example, based on atoms' period and group in the periodic table. fchl2018qx

In this work, we do not usek=4or higher-order interactions due to the computational cost from combinatorial growth of number of terms when enumerated directly,
which becomes a limiting factor for larger systems, such as in the  dataset.

More formally, Glielmo et al.gzv2018x,gzfd2020x define the order of a kernel of two local atomic environments as the smallest integerkfor which differentiating bykdifferent atomic coordinates always yields zero.
Conceptually, in our notation, the body-order of a global representation is

    k = _k∂^k+1|𝒜⟩/∂|α_i_1⟩…∂|α_i_k+1⟩ = 0

for all distinct|α_i_1⟩, …, |α_i_k+1⟩.
For local representations, one of the atomsα_iis fixed, reducing the order of the derivative by one.

This definition is not concerned with the extent to whichk-body terms utilisek-body information. For instance, in a local representation, products of 2-body terms, say, distances from the central atom to one other atom, depend on two distances (from the central atom to two other atoms), and are therefore formallyk=3, but cannot resolve angular information because the distance between the two other atoms is not known. gzv2018x,gzfd2020x

Powers ofk-body terms gzv2018x,gzfd2020x and products ofk,k'-body terms jkvak2020x can improve performance, but are less expressive than fullk^ζ-body andk k'-body terms, respectively.






§ REPRESENTATIONS


 §.§ Local atomic neighbourhoods

Local representations are computed for a local neighbourhood of a central atom,
usually defined as{   |α_i⟩ |   d_i ≤r_c }, where|α_i⟩denotes atomi,d_iis distance of|α_i⟩to the central atom, andr_c ≥0is a cut-off radius. 
Usually, a cut-off function that smoothly approaches zero ford_i →r_cis employed to prevent discontinuities at the threshold.

Both the  and  implementations of SOAP include the central atom in the neighbourhood, and thus in the neighbourhood density, dbcc2016qx
in contrast to the original definition bkc2013qx.

SFs do not take the central atom into account explicitly. 

In periodic systems, the unit cell is replicated up to the cut-off radius to ensure that all interactions within the neighbourhood are included.
In practice, some implementations may internally use a modified effective cut-off radius.
For instance,  ensures that atoms up to the tail of the radial basis function are taken into account. 



sectionsection

§ EMPIRICAL COMPARISON


 §.§  dataset

The  dataset, rdrl2014x,rdrl2015x also known as , contains 133 885 small organic molecules composed of H, C, N, O, F with up to 9 non-H atoms.
It is a subset of the “generated database 17” (). rdbr2012x
Molecular ground state geometries and properties, including energetics, are computed at density functional level of theory using the Becke 3-parameter Lee-Yang-Parr (B3LYP) sdcf1994x hybrid functional with 6-31G(2df,p) basis set.

We use the version available at <qmml.org>, which offers a convenient format for parsing, and exclude all structures in the  file and those listed in the  file as “difficult to converge”, as those are potentially problematic.
Total energies were converted to energies of atomization by subtracting the atomic contributions given in file .



 §.§  dataset

The  dataset, nrbsmrcwh2019x also known as , contains unrelaxed geometries and their enthalpies of formation for the 10 binary alloys AgCu, AlFe, AlMg, AlNi, AlTi, CoNi, CuFe, CuNi, FeV, and NbNi.
For each alloy system, unrelaxed geometries with lattice parameters from Vegard's rule v1921x,da1991x and energies are computed for all possible unit cells hf2008x with 1–8 atoms for FCC and BCC lattices, and 2–8 atoms for HCP lattices, 
using the generalized gradient approximation (GGA) of Perdew, Burke and Ernzerhof (PBE) with projector-augmented wave (PAW) potentials and generalized regulark-point grids wmm2016x,mjhh2018x.
The dataset contains 631 FCC, 631 BCC, and 333 HCP structures per alloy system, yielding 15 950 structures in total. 
We use the version available at <qmml.org>.



 §.§  dataset

The  dataset sgylbhglzs2019x is a Kaggle challenge kagglex dataset containing 3 000 ternary (Al_x-Ga_y-In_z)_2O_3oxides,x+y+z=1, of potential interest as transparent conducting oxides.
We predict formation and band-gap energies of relaxed structures, using either relaxed () or approximate () structures from Vegard's rule as input.
Geometries and energies are computed at the density functional level of theory using the PBE functional as implemented in the all-electron code bghhhrrs2009x with  settings.

The challenge scenario is to predict formation and band-gap energies of relaxed structures from unrelaxed geometries obtained via Vegard's rule.
This is equivalent to strong noise or bias in the inputs.
Unlike pure benchmarking scenarios, where computationally expensive relaxed geometries are given, the challenge scenario is closer to a virtual screening application in that Vegard's rule geometries are computationally inexpensive to obtain.

The dataset contains all structures from the challenge training and leaderboard data.
Unless otherwise noted, we report RMSE, not the root mean square logarithmic error used in the challenge.

*1

 §.§ Learning curves

Plots of empirical prediction errorϵas a function of training set sizenare called “learning curves”. 
Asymptotically, we assume the error to decay as a negative power, afs1992xϵ= a' n^-b. 
On a log-log plot,ϵis therefore linear,logϵ= a - b log(n), and the offseta = loga'and slopebcan be used to characterize predictive performance of models. hl2016x
For QM/ML models the estimated quantities are noise-free (except for numerical noise, which is negligible for converged calculations)
and representations are unique.

We use base-10 logarithms for learning curves.

For asymptotic fits we weight training set sizes by the standard deviation over their respective splits to attenuate for small sample effects,
as the above equation is valid only in the limitn →∞.
All learning curves in <ref> show linear behaviour after at most a few hundred samples.
See <ref> for estimated offsetsaand slopesb.



 §.§ Subsets

For training and validation, data subsets were sampled as follows:
An outer validation set[In the literature, the terms “test set” and “validation set” are sometimes used with different meaning.
To avoid confusion, we use “outer” for the subset employed to measure performance, and “inner” for the subset employed to optimize HPs.]
was randomly drawn (10 k molecules for , 1 k structures for , 600 structures for ).
From the remaining entries, outer training sets of sizes 100, 250, 650, 1 600, 4 000 and 10 000 for datasets ,  and 100, 160, 250, 400, 650, 1 000 and 1 600 for dataset  were randomly drawn.
These sizes were chosen to be equidistant in log-space.
Each outer training set was then split into an inner training set and an inner validation set by randomly drawing the latter.
We used an 80  / 20 % split, yielding inner validation sets of size 20, 50, 130, 320, 800, 2 000 for datasets ,  and 20, 32, 50, 80, 130, 200, 320 for .
The whole procedure was repeated 10 times.
We excluded structures with few atoms (6 or fewer non-H atoms for , 5 or fewer atoms per unit cell for , 10 atoms per unit cell for ) as there are not enough of these for statistical learning.



 §.§ Sampling

To reduce variance, remove bias and ensure that subsets faithfully represent the distribution of the whole dataset, subsets were drawn using Monte-Carlo sampling such that differences to the parent dataset in selected statistics were below pre-defined fractional thresholds.

For dataset , these were number of N, O and F atoms, number of molecules with 7, 8 and 9 non-H atoms, binned number of atoms (with H), and binned energy.
For dataset , these were number of all constituting elements, unit cells with 6, 7, 8, and 9 atoms, binned sizes and energies. 
For dataset , these were number of Al, Ga, In, O atoms, unit cells with 20, 30, 40, 60, 80 atoms, and binned energies.



 §.§ Kernel regression

We use kernel ridge regression htf2009x or Gaussian process regression rw2006x (the two are equivalent in terms of predictions).
A detailed derivation can be found in r2015fx.
In summary, predictions are basis set expansions of the formf(x⃗) = ∑_i=1^n α_i k(x⃗_⃗i⃗,x⃗),
wherex⃗is the system to predict,x⃗_⃗1⃗,…,x⃗_⃗n⃗are the training systems, andkis a symmetric positive definite function (kernel).
The regression coefficientsα⃗are obtained by minimizing the regularized quadratic loss∑_i=1^n ( y_i - f(x⃗_⃗i⃗) )^2 + λ||f||^2_ℋ,
wherey⃗are property values of the training data and the regularization strengthλis a HP that controls the smoothness of the predictor.
In this work, we use the Gaussian kernelk(x⃗,z⃗) = exp(-||x⃗-z⃗||_2^2 / 2σ^2), where the length scaleσis a HP (<ref>).
We used the r2015fx,qmmlpackx implementation.



 §.§ HP optimization

For model selection we optimized the HPs of representations, kernel and regression method, 
including structural HPs (for example, whichk-body functions to use) and numerical HPs (for example, the Gaussian kernel length scale).
Specifically, the RMSE of an inner validation set was minimized using tree-structured Parzen estimators bbbk2011x,byc2013x in combination with local grid search.
The same optimization scheme was used for all representations, using consistent grid spacings and parameter ranges to reduce human bias.
Our corresponding  package cmlkitx provides interfaces to the  optimization package bbbk2011x and to each representation's implementation(s);
it is freely available under an open source license. 

The space of possible models (“HP search space”) is a tree-structured set of choices, for instance, between differentk-body functions, or different values of a numerical HP.
Tree-structured Parzen estimators treat this search space as a prior distribution over HPs, updated every time a loss is computed to increase prior weight around HP settings with better loss. 
We use uniform priors throughout, discretizing numerical HPs on logarithmic or linear grids as necessary.
Once a HP search space has been defined, model selection is fully automatic.

HPs were optimized for each training set size as follows:
For each trial, representation HPs and starting values for regression method HPs (Gaussian kernel length scale and regularization strength) were drawn from the prior. 
The latter were then refined through a randomized local grid search and the resulting HP values used to update the prior.
All optimizations were run for 2 000 steps, and rerun three times, to minimize variance from stochastic optimization.

To reduce computational cost, HPs were optimized on only one outer split.



 §.§ Kernel regression HPs

We used KRR with a single Gaussian kernel, a frequently used combination in the literature.
Note that due to Requirement <ref>, the Gaussian kernel is better suited than less smooth kernels such as the Laplacian kernel. hr2017x

No post-processing of the kernel was performed.
In particular, centering of kernel and labels, which together is equivalent to having an explicit bias termbin the regression, were not performed, as this is not necessary for the Gaussian kernel. pmrrv2001x
Depending on the representation used, labels were normalized for training as needed to either represent values per atom or per entire system. (<ref>)

This setup entails two HPs:
The widthσof the Gaussian kernel, and the regularization strengthλ.
Search spaces for these two HPs (<ref>) were held constant across all representations and learning curves.

See reproducibilityx for HP search spaces and optimized model HPs. 

*2

 §.§ Symmetry function HPs

The five SFs proposed in b2011ex are


    G^1_i =   ∑_j f_c(d_ij)  
    
    	
    	G^2_i =   ∑_jexp(-η (d_ij-μ)^2 )   f_c(d_ij) 
    
    	
    	G^3_i =   ∑_jcos ( κ  d_ij )   f_c(d_ij)  
    
    	
    	G^4_i =    2^1-ζ∑_j,k ≠ i (1 + λcosθ_ijk)^ζ
       exp(-η (d_ij^2 + d_ik^2 + d_jk^2))
        f_c(d_ij)   f_c(d_ik)   f_c(d_jk) 
    
    	
    	G^5_i =    2^1-ζ∑_j,k ≠ i (1 + λcosθ_ijk)^ζ
       exp( -η (d_ij^2 + d_ik^2) )  f_c(d_ij)   f_c(d_ik)  ,


with the cut-off function


    f_c(d_ij) = 
    		0.5 cos(π d_ij / c )    for   d_ij≤ c
    
    			0    for   d_ij > c   .


In <ref>,
indexiis the central atom,j, krun over all atoms in the local environment aroundiwith cut-off radiusc,d_lmindicates pairwise distance,θ_lmnthe angle between three atoms;ηandκare broadening parameters,μa shift,ζdetermines angular resolution.λ= ±1determines whether the angular part ofG_i^4andG_i^5peaks at0or180. 

We utilize the arunnerx,b2015qx,b2017qx software to compute SFs and restrict ourselves to the radial SFsG_i^2and angular SFsG_i^4( functions 2 and 3). 
We use the same SFs for all element combinations to minimize size of HP search space.
This can result in constant or low-variance features,
which are unproblematic for kernel regression with a Gaussian kernel (<ref>) as they enter only through the norm of feature vector differences
(for neural networks, these features could be more problematic).

Similarly, we use an empirical parametrization scheme gsbbm2018x to choose HPsμandηforG_i^2and HPsη, ζ, andλforG_i^4. 

For radial SFs we use two schemes, shifted and centered. 
For shifted,μis chosen on a linear grid whileηis held fixed. 
For centered,μ=0andηis chosen such that the standard deviation of each SF lies on the same grid points. 
Fori ∈{0, 1, …, n }a point on a one-dimensional grid,Δ= c - 1.5/n-1, andr_i = 1 + Δi, 
in the centered scheme,μ_i = 0andη_i = 1/2 r_i^2, 
and in the shifted scheme,μ_i = r_iandη_i = (2 Δ^2)^-1.
In this setting, the only HP is the number of grid pointsn+1, which we allow to vary from 2 to 10 for each scheme.

For angular SFs, we chooseλ=±1andζ= 1, 2, 4.
The only HP remaining is the broadeningη, optimized on alog_2grid between -20 and 1 with spacing 0.5.

The radial SFs and two angular SFs withλ=±1andζ= 1are always included, but the optimizer can enable or disable any of the remainingk=3SFs withλ=±1andζ= 2, 4.

Cut-off radii are varied in integer steps, starting from the integer above the smallest distance found in the dataset, up until at most 9.

The output of  is post-processed to be suitable for KRR, placing all SFs for a given type of central atom in separate blocks of an atomic feature vector with# elements ×# SFscomponents. 
For the Gaussian kernel, this leads to negligible kernel values between representations belonging to different elements.

As SFs are local representations, labels were normalized to (extensive) per-system values.

See reproducibilityx for HP search spaces and optimized model HPs. 



 §.§ Many-body tensor representation HPs

We employed the MBTR implementation in , qmmlpackx 
adding optional normalization byℓ_1orℓ_2norm.
Fork=2,3, representations fork=2andk=3were concatenated.
MBTR exhibits several categorical HPs, with subsequent numerical HPs conditional on prior choices.

We used thek-body functions ,  (k=2), and , ,  (k=3).
No one-body terms were used as atomization and formation energies already contain linear contributions of element counts.

Histogram ranges were chosen based on the whole dataset, as inter-atomic distance ranges are similar for all subsets.
100 discretization bins were used throughout. 
Broadening parameters were restricted to at least a single bin and at most a quarter of the range of the corresponding geometry function.

From the weighting functions, we used , ,  (k=2), and , ,  (k=3). 
The latter two in each set introduce conditional HPs. 
For periodic systems, in particular the  dataset, the ranges of these parameters were manually restricted to avoid excessive computation times (above 30 for one trial). 
The convergence threshold was set to 0.001.

We used the  indexing scheme, which generates all permutations of elements (as opposed to , which does not double-count element combinations, for example, CH and HC).
This seems to lead to more consistent behaviour and higher predictive accuracy for supercells, or unit cells of different sizes, 
and similar accuracy for molecules, at the expense of higher computational cost.

We used per-system energies for the  dataset and per-atom energies for datasets  and .



 §.§ Smooth Overlap of Atomic Positions HPs

We used the  implementation of SOAP with Gaussian-type orbitals, adscribex,hjmfrgrf2020qx 
which we found to provide more accurate predictions at lower computational cost than the aquippyx implementation.

Results are already structured by element types; 
no post-processing was applied.

HPsl_maxandn_maxwere chosen between 2 and 8. 
Cut-off radii were chosen as for SFs, and the broadening adapted to the resulting ranges (53 steps from -20 to 6 on alog_2grid).
We report results for the  basis set, which resulted in lower prediction errors than the  one, and was faster to compute. 
Labels were normalized to per-system values.

See reproducibilityx for HP search spaces and optimized model HPs. 



 §.§ Prediction errors
<Ref> present numerical values underlying the learning curves for RMSE (<ref>) and MAE (<ref>).
For rRMSE (<ref>), standard deviations of 
239.31, 
178.86, 
104.57, 
were used for datasets , , , computed over the whole dataset 
(differences to standard deviations over validation sets were around 1 % or less in all cases).



 §.§ Error metrics

We measure predictive performance by two metrics, an absolute one and a relative one that facilitates comparison across datasets.
In addition, we also provide a metric for qualitative comparison with the literature.

Lety_i,f_i,e_i = f_i - y_idenotei-th observed label, prediction and residual.
Root mean squared error (RMSE) and mean absolute error (MAE) are given by 


    RMSE = √(1/n∑_i=1^n e_i^2 ), 
    	  MAE = 1/n∑_i=1^n |e_i| .


The canonical loss for least-squares regression is RMSE (as it is optimized by the regression).
We also provide MAE since it is often reported in the literature (<ref>).

RMSE and MAE are scale-dependent, and thus not suited for comparison across different datasets.
We therefore also report the scale-independent relative RMSE (rRMSE),


    rRMSE = RMSE/√(1/n∑_i=1^n (y_i - y̅)^2 )
    	= RMSE/σ_y
    	= RMSE/RMSE^∗ ,


wherey̅ = 1/n ∑_i=1^n y_iis the mean of the observed labels andσ_yis their standard deviation.
The rRMSE can be seen as RMSE relative to the RMSE of a baseline model RMSE^∗that always predicts the mean of the labels.

While the latter is more naturally computed using training labels and the former using validation labels, 
as long as as the assumption of independent and identically distributed data holds, the number of samples is more important.

See ctg2017x,b2019x and references therein for an extended discussion of error metrics.



 §.§ Compute times
<Ref> present empirical computational costs, measured by processor wall-time, for calculating representations and kernel matrices, respectively.

Experiments were run on a single core of an Intel Xeon E5-2698v4 2.2 processor.

For <ref>, representations of the 10 k, 1 k, 600 systems (datasets , , ) in the first outer validation set were computed en bloc and the result divided by number of systems; this was repeated three times.

Similarly, for <ref> kernel matrices between the representations of these systems were computed, also over three repetitions. 
The results were divided by the number of entries in the respective kernel matrices, yielding average kernel evaluation times.

<Ref> presents a summary overview of compute times for representations and kernel matrices.



 §.§ Analysis details

Predictive accuracy as measured by rRMSE is worse for solid-state datasets compared to the molecular  one. 
This might indicate that periodic systems pose harder learning tasks than molecules.

MBTR performs worse for solid-state datasets than for the  one, in particular for . 
This might be due to increasing difficulty of the learning problem with system size (see discussion in <ref>) and 
lack of intrinsic scaling with number of atoms, impeding interpolation between unit cells of different size.
The high computational cost of MBTR withk=3for large periodic systems also renders HP optimization more difficult.

For the  dataset at1 600training samples, we observe an increase in RMSE standard deviation compared to neighbouring training set sizes for most methods. 
Comparing to MAE, which exhibits no such effect, and investigating errors individually, revealed that this is due to outliers, that is, few predictions with high error in some, but not all, outer splits. 
The problematic structures are ring molecules, and are not present in the outer training split used for HP optimization. 
This stresses the importance of carefully stratifying benchmark datasets.



 §.§  dataset
<Ref>, and, <ref>
present results for energy predictions on the  dataset, that is, the  dataset with approximate geometries obtained from Vegard's rule. 
In contrast to relaxed structures, such geometries can be obtained at almost no cost, and could be used in virtual screening campaigns.

We observe
(i) a strong increase in prediction errors (14–21 % for rRMSE),  
(ii) collapse of all representations to similar performance,
(iii) large differences between MAE and RMSE, indicating significant outliers.

From this, we conclude that 
the map from unrelaxed structures to ground-state energies is harder to learn than the map from relaxed structures to their energies, and,
that here the representation is not the limiting factor, and other sources of error dominate.



 §.§ Comparison with literature-reported errors

Due to different conditions, such as sampling, regression and HP optimization methods, comparisons with performance estimates reported in the literature must remain qualitative.
Frequently, only MAE is reported, which tends to result in lower absolute values and to de-emphasize outliers.
<ref> presents selected performance estimates from the literature.
Overall, errors in this work appear to be compatible with reported ones.



 §.§ Comparison with DFT and experimental errors

The error of DFT simulations against experimentally measured observations depends on system and property, as well as choice of density functional and other parameters, such as convergence thresholds andk-point density.

For heats of formation and the Becke 3-parameter Lee-Yang-Parr (B3LYP) functional used for the  dataset, 
(systematic) MAEs relative to experiment of≈2.6 have been reported for small organic molecules containing only C, H, N, O. tj2008x

For cohesive energies and the Perdew-Burke-Ernzerhof (PBE) functional used for the  and  datasets,
values of approximately 200 to 300 / have been reported. l2008x,lvvc2014x,zrts2018x

For the PBE functional, reported MAEs in computed energies between different parametrizations of DFT codes and  between 20 different DFT codes on 71 elements in bulk crystalline form were approximately 2 meV/atom and 1.7 meV/atom, respectively; lbbbbbcccdetal2016x
the latter reduces to 0.6 meV/atom for all-electron codes only.

The best models for bulk crystal reported here have RMSEs of 4.6 meV/atom and 3.3 meV/cation on the  and  datasets.

However, the former benchmark values are integrated over a± 6 %interval around the equilibrium volume, 
whereas the values reported here are computed at the minima themselves and therefore measure related but distinct quantities.
This suggests that prediction errors are at least≈2–6 times larger than DFT-intrinsic variations.

section99tocsectionReferencesabbreviations,supplement