Performant implementation of the atomic cluster expansion (PACE): Application to copper and silicon

The atomic cluster expansion is a general polynomial expansion of the atomic energy in multi-atom basis functions. Here we implement the atomic cluster expansion in the performant C++ code \verb+PACE+ that is suitable for use in large scale atomistic simulations. We briefly review the atomic cluster expansion and give detailed expressions for energies and forces as well as efficient algorithms for their evaluation. We demonstrate that the atomic cluster expansion as implemented in \verb+PACE+ shifts a previously established Pareto front for machine learning interatomic potentials towards faster and more accurate calculations. Moreover, general purpose parameterizations are presented for copper and silicon and evaluated in detail. We show that the new Cu and Si potentials significantly improve on the best available potentials for highly accurate large-scale atomistic simulations.


I. INTRODUCTION
Atomistic modelling and simulation requires efficient computation of energies and forces. In recent years, machine learning (ML) based interatomic potentials, parameterized to large data sets of reference electronic structure calculations, have provided particularly successful surrogate models of the atomic interaction energy. The ML models construct representations of atomic structure that are used in various regression algorithms to predict energies and forces.
The recently developed atomic cluster expansion (ACE) 1 provides a complete and efficient representation of atomic properties as a function of the local atomic environment in terms of many-body correlation functions. Because of the completeness of the ACE basis 2 these may be employed directly using linear regression for the computation of energies and forces. Furthermore, using simple nonlinear embedding functions, ACE can represent many classical as well as ML interatomic potentials. For example, the widely used family of Embedded Atom Method (EAM) 3 and Finnis-Sinclair (FS) 4 potentials may be viewed as a lowest order ACE. Other properties, for example the moments of the density of states, may also be represented, and recursion or moments-based potentials like the bond-order potentials 5,6 expanded in the form of an ACE.
Moreover, there are deep connections between ACE and several ML representations and formulations. The only other known complete parameterizations, the Moment Tensor Potentials (MTP) 7 and the ML potential of Seko et al. 8 , are both based on a body-ordered invariant polynomial basis and can be exactly represented by ACE by suitable choice of hyperparameters and an explicit linear transformation. In addition, the Spectral Neighbor Analysis Potential (SNAP) 9 , the atomic Permutation Invariant Potentials (aPIPs) 10 and descriptors such as the symmetry functions of Behler 11 and the Smooth Overlap of Atomic Positions (SOAP) 12 can be obtained as special cases or minor variations of the ACE formalism; see Refs. 1,2,13 and the supplementary information (SI) for further details.
Here, we present the Performant Implementation of the Atomic Cluster Expansion (PACE) enabling efficient evaluation of ACE models within the LAMMPS molecular dynamics simulation software package 15,16 . We demonstrate in Figure 1 for two representative elements, Cu and Si, that PACE lowers the Pareto front of accuracy vs. computational cost that was established for several stateof-the-art ML potentials 14 . The details of how these were constructed are provided in the SI. While these benchmarks establish advanced computational performance, we also demonstrate the capacity of the ACE framework to develop highly accurate parameterisations: We present two novel parameterisations of interatomic potentials for Cu and Si that outperform available ML based potentials in terms of performance, accuracy and generalisability.
The fundamental building block from which ACE models are built are atomic properties ϕ i which are expanded in terms of body-ordered functions from the set of neighbors of each atom i, where Φ v are ν-order basis functions (each involving coordinates of ν neighbors) andc v the model parameters.
It appears as if this incurs an O(N ν ) computational cost, where N denotes the number of interacting neighbors; however, ACE exploits a much faster evaluation strategy that makes it possible to compute efficiently high bodyorder terms. This is achieved by (i) projecting the atomic  (top) and Si (bottom) for ACE potentials compared to a recent benchmark study 14 . The timings from Zuo et al. 14 were reduced by constant factors 0.55 (Cu) and 0.60 (Si) to correct for hardware differences and the new ACE timings then overlaid.
on atomic basis functions, φ v (r), resulting in and (ii) choosing a tensor product basis which leads to 1 j1,...,jr We call this reformulation the "density trick" (also used by Bartók et al. 17 and Shapeev 7 in formulating SOAP and MTP, respectively) and it results in the computational cost of an atomic property ϕ i scaling linearly in N (due to evaluating the A ik ) and also linearly in ν (due to evaluating the correlations). Furthermore, in Sec. III C we present an evaluation scheme that avoids the ν-scaling altogether. An ACE model may be defined in terms of several atomic properties ϕ (p) i , p = 1, . . . , P , for each atom i. For the simplest linear model of the potential energy one would use just one property, the atomic energy E i , A more elaborate model may generalize the pairwise repulsion and the pairwise density of the Finnis-Sinclair potential 4 to arbitrary many-atom interactions, i .
In general, a large number of different atomic properties that are regarded as descriptors enter a non-linear function where the non-linearity F may be explicit as in the Finnis-Sinclair model, or represent a general approximator such as artificial neural networks, as used by Behler and Parrinello 18 , or a kernel ridge regression model as used in the Gaussian Approximation Potential (GAP) 17 . Different non-linearities F may be used to incorporate physical or chemical insights in bond formation. Since the d-shell of copper is nearly full, angular contributions are generally small in the bulk, hence copper is modelled well by classical central-force functionals with non-linear EAM or FS type embedding functions that effectively generate high body-order terms 3,4,19 . Our parameterization for copper therefore starts from the FS representation of the energy, as in Eq. (7), but with the two atomic properties not limited to pairwise terms but including many-atom contributions that capture small angular contributions in the bulk and larger angular contributions in small clusters or two-dimensional structures.
On the other hand, the diamond structure of silicon is stabilized by angular contributions over close-packed structures, which highlights the importance of interactions beyond pairwise terms. Many different angledependent potentials have been developed for Si. Perhaps the best known are the Stillinger-Weber potential 20 with a linear three-body term and the Tersoff potential 21 which includes non-linear functions of three-body contributions. The most accurate potential for silicon to date, the SOAP-GAP model of Bartók et al. 22 is an intrinsically high body-order potential. Here, we present a linear ACE for Si, which may be viewed as a generalization of this potential that includes all body-order interactions up to some maximum. In this way, ACE is employed in its basic form shown in Eq. (6), which simplifies the parameterization considerably and avoids implicit assumptions on the form of non-linear terms that are often present in ML frameworks.
We carry out a detailed comparison of both our ACE parameterizations to the most reliable models available from literature. For Cu, we compare to the EAM potential by Mishin et al. 23 , to a recent SNAP 24 parameterisation as well as the GTINV 8 ML potential. For Si, we compare to the GAP that was shown to reproduce a wide range of observable properties for crystalline, liquid and amorphous Si phases 22 .

A. Reference data
The parameterization for Cu was obtained by matching to the energies and forces of about 50000 total energy calculations as obtained with density functional theory (DFT) using the PBE 25 functional as implemented in the FHI-aims code 26,27 . The reference data included small clusters, bulk structures, surfaces and interfaces, point defects and their randomly modified variants. Part of the reference data has been briefly described in Ref. 1, but has been extended significantly for the present parameterization. We employed pyiron 28 for generating part of the reference data.
The parameterization for Si was obtained by fitting to the same extensive silicon database GAP was fit to 22 . The database covers a wide range of configurations including crystalline structures, surfaces, vacancies, interstitials and liquid phases. The DFT reference data were generated using the CASTEP 29 software package.

B. Parameterization and timing
We used different parameterization strategies for Cu and Si motivated by their different bond chemistry. In particular, the Si parameterization was obtained from solving a linear system of equations, whereas the Cu fit required non-linear optimization. The Cu potential has a total number of 2072 parameters, of which 756 are expansion coefficients for each of the two densities, and 560 parameters are used for the radial functions. The DFT reference showed that interactions are smaller than 1 meV when atoms are further than the cutoff distance r c = 7.4Å apart, when rigidly separating slabs. Parameter optimization led to a fit with an error of 3.2 meV/atom for structures that are within 1 eV of the ground state. This fit was then fine-tuned towards structures close to the ground state, which further decreased the error to 2.9 meV/atom and slightly increased the error of higher energy structures. For Si we used a total of 6827 basis functions parameterized as a linear model, with a maximum body order corresponding to ν = 4. These basis functions were selected using the construction outlined in Dusson et al. 2 . We show the silicon ACE matches the accuracy of the general-purpose GAP potential introduced by Bartók et al. 22 . More specifically the energy error for the ACE model was found to be 1.81 meV/atom for structures within 1 eV from the ground state. The corresponding errors for the GAP model are 1.25 meV/atom on the silicon database presented in Bartók et al. 22 .
To evaluate the computational efficiency of PACE we carried out molecular dynamics (MD) simulations for face-centered cubic (fcc) Cu and diamond Si structures. We found that a single force call takes 0.32 and 0.80 ms/atom, respectively, for the Cu and Si ACE models 30 . These speeds are sufficiently fast for large scale MD simulations and Monte Carlo sampling, for example, for the computation of phase diagrams. The efficiency of PACE is about two orders of magnitude slower than empirical potentials.

C. Copper
The ACE for Cu has been comprehensively validated against DFT and available experimental data and compared to three other Cu potentials. The potentials we chose for the comparison were (i) the EAM potential of Mishin et al. 23 , which exhibits an excellent overall accuracy and is considered as the reference Cu EAM potential, (ii) the SNAP model of Li et al. 24 , which was trained to strained crystalline as well as melted Cu phases obtained by ab-initio MD, and (iii) the ML interatomic potential Cu-gtinv-934 (GTINV) of Seko et al. 8 , which was fitted to an extended DFT database of 10 4 structures and reached RMSE values of 8.2 meV/atom. The EAM and SNAP potentials were computed through the OpenKIM API 31 .
We evaluate the models for a broad range of structures and properties that not only exceed beyond the reference data but are also relevant for observable macroscopic behavior of Cu. Fig. 2 gives an overview of the binding energy over large volume changes. ACE provides a very good match to the reference data at all distances, while the shorter range of EAM means that interatomic interactions are cut off too early when the atoms are separated. The even shorter cutoff of SNAP leads to abrupt bond breaking, illustrating that the cohesive energy was not fitted in the construction of the potential and therefore one cannot apply the potential, for example, for gas phase condensation simulations. GTINV shows significant oscillations at larger interatomic distances. These observations also apply to the dimer shown in Fig. 9.

Bulk properties
a. Structural energy differences A detailed analysis of structures that are energetically close to the fcc ground state is presented in Fig. 3. All potentials reproduce correctly the structural order of fcc → dhcp → hcp → bcc. The energy minima predicted by EAM and GTINV potentials are shifted to smaller volumes, which may be due to different DFT reference data. EAM and SNAP also show larger discrepancies for the fcc-bcc energy difference. b. Elastic moduli The elastic moduli for the ground state fcc structure are summarized in Tab. I. ACE and EAM reproduce the DFT reference very well, while small deviations are observed for SNAP and slightly larger for GTINV. Similar outcomes are obtained for other bulk phases (see SI) with ACE and EAM giving consistently the best agreement with DFT.
c. Phonons Figure 4 shows a comparison of phonon band structures and densities of states (DOS) for fcc Cu. Despite not having fitted any phonon frequencies explicitly, ACE provides the best match to the reference DFT data. The EAM and GTINV potentials overestimate the width of the DOS , while SNAP underestimates it. These conclusions apply also for the phonon DOS of other crystal structures that are shown in the SI.
d. Structural transformations Transformations between different crystal structures present a sensitive test for any interatomic potential as both bond distances and bond angles are changed simultaneously. In addition, the associated changes in atomic coordination effectively scrutinize the screening of pairwise terms by many-atom contributions. As shown in Fig. 5, all potentials agree well with the reference DFT data for the tetragonal, trigonal and hexagonal paths. However, only ACE provides an excellent quantitative interpolation for all structures along all considered transformation paths. Especially, the orthorhombic transformation, which can be regarded as a generalization of the Bain path 39 , is challenging for the other potentials.
e. Melting transition and thermal expansion We used thermodynamic integration to evaluate the free energy of the solid and liquid phases of Cu. The free energies intersect at T = 1272 K, about 20 K above the 1251 ± 15 K predicted by DFT 37 . The EAM and SNAP melting temperature at T = 1325 K and 1372 K are close to the experimental value of 1358 K. The prediction of the melting point with GTINV was not possible due to long evaluation times and the lack of a parallel implementation. Fig. 6 shows the thermal expansion as obtained from MD simulations in the NPT ensemble. All models agree well with the experimental data for temperatures up to 600 K and exhibit minor deviations at high temperatures.

Interfaces and surfaces
Planar defects include internal interfaces, such as stacking faults (SF) and grain boundaries (GBs), where 0 500 1000 T (K) the local atomic density does not vary significantly but bond angles change compared to bulk. In contrast, at free surfaces the bond angles remain mostly unaltered but the surface atoms loose about half of their neighbors. Typically, central-force models such as EAM provide a good description of structures and energies of GBs but cannot capture well the large local density changes at surfaces which usually leads to underestimation of surface energies.
a. Stacking faults The small energy differences between the close-packed fcc, hcp and related structures in Cu imply small SF energies. ACE predicts the SF energies in very good agreement with DFT reference data, as shown in Tab. I, with comparable predictions from EAM. GTINV predicts negative stacking fault energies, hinting at a different ground state. SNAP provides stacking fault energies with slightly larger deviations from the reference data.
b. Grain boundaries The energies of several twin and twist symmetric GBs (Σ = 3, 5, 9) are compared in Fig. 7 to reference DFT data from the Materials Project database 41 (see SI for more details). As expected, all potentials predict the GB energies very accurately, which suggests good transferability of all models for environments with small local density variations.
c. Surface energies As noted above, surfaces present a much more stringent test than GBs. ACE provides the best agreement with DFT reference data for all low-index surfaces, as shown in Fig. 8, while both SNAP and EAM consistently underestimate the surface energies. For GT-INV we observed a detachment of the top surface layers during relaxation which resulted in unphysically high surface energies that were excluded from the comparison.
d. Bond breaking In addition to the energetics of surfaces, we examined bond breaking in various atomic environments. Such tests have practical relevance as they are related to fracture, surface adsorption or vaporization. We designed three distinct decohesion tests that  are schematically shown in Fig. 9. These tests compare bond dissociation in the Cu dimer, detachment of a Cu adatom from the (111) surface, and an ideal rigid decohesion of bulk Cu slabs that leads to the formation of two (111) free surfaces. As can be seen from Fig. 9, ACE is the only model that is able to describe quantitatively accurately the impact of the atomic environment on bond breaking. The presence of neighboring atoms leads to an effective screening of the interatomic bonds and their interaction ranges 42 . The dimer and the adatom have no neighbors so that their interaction range is longer than the interaction between two surfaces whose atoms are surrounded by bulk. Given the simplicity of EAM, it provides a surprisingly good account of bond breaking in the very different environments, while GTINV and SNAP have problems with this test. configurations were part of the reference data, ACE reproduces the vacancy formation energy very well while the other potentials overestimate the DFT reference by 0.1-0.3 eV, see Tab. I. The migration barrier is reproduced well, too, by the models, apart from SNAP that overestimates the barrier.
None of the interstitial configurations were included in the ACE training set, but ACE results are consistent with those of the other potentials and together with EAM agree best with recent DFT results. 35 The 100 dumbell is predicted to have the lowest energy, followed by the octahedral and tetrahedral configurations. These predictions are consistent with those for other fcc metals. 43 b. Small clusters Small metallic clusters, important for catalysis and nanotechnology, usually form a large number of isomers with energies and structures often governed by subtle electronic structure contributions. For this reason, the predictions of the detailed energetics and structural stability is very challenging for interatomic potentials that are typically aimed at the description of bulk systems. We compared the predictions of ACE and the other models for three-and four-atomic clusters.
For the Cu trimer, the ground state structure is an isosceles triangle configuration while the linear trimer corresponds to an energy saddle point and is not dynamically stable 44 . In fact, the linear trimer transforms to a metastable configuration of a bent molecule with an obtuse angle of 130 • . ACE is the only model that correctly reproduces the instability of the linear trimer and the existence of the metastable bent configuration. EAM predicts the equilateral triangle as the only stable configuration while for SNAP and GTINV both the linear trimer and the equilateral triangle are stable configurations. The energy differences between the configurations are also reproduced most accurately by ACE while the other models either significantly underestimate (GTINV) or overestimate (EAM, SNAP) the DFT values.
In the case of the tetramer, only ACE and GTINV give correctly the planar equilateral rhombus 44 as the ground state, albeit GTINV shows also additional metastable configurations. Both EAM and SNAP favor incorrectly the close-packed tetrahedron which may originate from the lack of or weak angular contributions.
c. 2D structures Planar 2D structures belong to a family of structures that is usually not included in the validation of interatomic potentials for bulk metals. It has been found recently that Cu is the only metal whose free-standing monolayers arranged in honeycomb, square and hexagonal close-packed lattices are dynamically stable 45 . We investigated in detail the 2D hcp lattice; both EAM and SNAP potentials show dynamic instabilities related to out-of-plane atomic displacements for the 3×3×1 supercell that we used in our calculations (Fig. 10). In contrast, DFT and ACE predict real phonon frequencies that confirm excellent transferability of ACE once more. We note that the 2D hcp structure could not be stabilized using the GTINV potential.

D. Silicon
The ACE for silicon was created by fitting to an extensive database first introduced to create a general-purpose GAP model for silicon 22 . This GAP was shown to describe silicon accurately and to also be a qualitatively better interatomic potential than all other models tested, each best in their class: Stillinger-Weber 20,46 , EDIP 47 , Tersoff 21 , MEAM 48 , DFTB 49 and ReaxFF 50 . In this paper we show that the silicon ACE potential achieves the same accuracy as the GAP model, while being around 30 times faster in evaluation time and also better at extrapolating to unseen configurations.
The following section presents a benchmarking of the ACE silicon potential on a wide range of properties including bulk, surface, liquid and amorphous properties as well as a random structure search (RSS) 51 test.

Bulk properties
a. Structural energy differences The energies of the diamond, hexagonal diamond, β-Sn, bc8, st12, bcc, fcc, simple hexagonal (sh), hcp and hcp' are compared to DFT in Fig. 11. Excellent agreement with the DFT reference is observed for all structures apart from hcp'. Si hcp has two minima 22 , the conventional hcp with c/a ≈ 3/2, and hcp' with c/a < 1. The hcp' crystal structure is not contained in the DFT reference silicon database, however, both ACE and GAP predict the minimum. The GAP predicts the DFT reference energy at the minimum more accurately than ACE, while the latter gives a better estimate of the curvature.
The energy versus volume curves for the silicon diamond and bcc are extended over a wide volume range in Fig. 12. Both potentials accurately describe the minima around 15 and 20Å 3 /atom for bcc and diamond, as previously shown in Fig. 11. At larger volumes GAP exhibits unphysical high-energy minima. ACE does not show these minima and is close to the DFT reference, demonstrating better extrapolation compared to GAP. This extrapolative behavior is remarkable since there is no reference data at these large volumes as shown by the data density in the lower panel.
b. Elastic moduli The elastic constants for Si in the diamond structure are summarized in Tab. II. Both ACE and GAP match the DFT reference within a few percent.
c. Phonons The phonon dispersion for ACE and GAP is compared against the DFT reference in Fig. 13. Both GAP and ACE accurately describe the phonon spectrum in comparison to the DFT reference, with the band width of GAP showing a better match to DFT. More silicon phonon spectra are found in the SI.
d. Thermal expansion We investigated the thermal expansion, Grüneisen parameter and heat capacity of ACE and GAP in the quasi-harmonic approximation as shown in Fig 14. Diamond silicon displays negative thermal expansion at low temperatures 52 , which ACE models very well compared to the DFT reference, and more accurately than GAP. The heat capacity is described with almost perfect agreement, whereas the thermal expansion saturates at a slightly too high value for high temperature (for both GAP and ACE).
e. Liquid phase ACE was also tested in a liquid simulation on an eight-atom 2x2x2 supercell (64 atom) and  compared to GAP and the DFT reference. The radial distribution function (RDF) and angular distribution function (ADF) where averaged over 20000 MD steps (0.25 fs timestep). The DFT reference data was generated using CASTEP averaging over 9700 MD steps (0.25 fs timestep) taken at T =2000K, using a 200 eV plane-wave energy cutoff and 0.05Å −1 k-point spacing. The results are shown in Fig. 15   ing from the melt. Here we quench a 216-atom sample of liquid Si from 2000 K to 500 K at a rate of 10 12 K/s with a 1 fs time step (1.5 × 10 6 steps) using the LAMMPS software 16 . After the MD steps the final configuration was relaxed to a local minimum with respect to cell size and shape. The radial distribution functions (RDF) of both GAP and ACE are compared to experimental results 53 (since DFT results are not computationally feasible) in Fig.16. Both GAP and ACE accurately describe the first and second neighbor peaks, and no atoms in the range (2.5Å≤ r ≤ 3.25Å).  Radial distribution function (RDF) for a 216-atom configuration generated by cooling liquid Si from 2000K to 500K comparing ACE, GAP and experimental data unrelaxed surface and bulk crystal diamond were part of the database and accurately fitted, as well as some configurations along the path. The ACE energy is significantly smoother than GAP and closer to DFT reference.

Point defects
The diamond vacancy formation energy and interstitial formation energies including tetragonal, hexagonal and dumbbell are shown in Table. II. Both ACE and GAP predict point defect formation energies very well. a. Fourfold defect The lowest formation energy point defect is the "fourfold-coordinated defect" which consists of a bond rotation followed by a reconnection of all broken bonds 54 . We performed the following test using the ACE model: optimise the defect structure (using a 64 atom cell) with DFT, then re-optimise it with ACE, and finally compute the minimum energy transition path to the perfect crystal. When this test was performed with GAP in Bartók et al. 22 , no local minimum was found corresponding to the defect. With ACE however, we do find a local minimum, as shown in Fig. 18, where we also show the energy of the path evaluated with DFT and GAP. Remarkably, while both GAP and ACE make a similar error near the transition state, the ACE energy is significantly better for the relaxed defect structure, leading to the stabilisation of the defect. Note that there are no configurations in the fitting database (which is identical for GAP and ACE) near the defect structure and the transition state, so again this shows the extrapolation power of the ACE model.
b. Random structure search In the random structure search (RSS) 51,55 test, randomized atomic configurations are relaxed, providing a view of the fitted potential energy surface for higher energy configurations. The RSS tests were performed on eight-atom configurations with close to cubic initial shapes and initial interatomic distances > 1.7Å. These configurations were then relaxed using the two-point steepest-descent method 56 . The resulting energy per atom versus volume per atom distribution is shown in Fig. 19. ACE shows a similar distribution compared to DFT, with the diamond structure at the correct volume and a few structures up to 0.2 eV per atom higher at comparable or somewhat larger volumes. A larger group of configurations is found at higher energies over a wider distribution of volumes. The density of states on the right panel of Fig. 19 shows excellent agreement with DFT, as does the GAP model. This test is a strong discriminator between potentials, with all empirical potentials tested in Bartók et al. 22 failing completely.

E. Discussion
We present a performant implementation of ACE in the form of the PACE code. We demonstrate that ACE, as implemented in PACE, shifts the Pareto front to higher accuracy and faster evaluation times, as compared to a number of machine learning potentials from Ref. 14. For our general purpose parameterizations of Cu and Si the CPU time per atomic force call is below 1 ms. As our implementation is fully compatible with LAMMPS, large scale simulations become possible, which we demonstrate through the computation of the free energies of liquid and solid phases for evaluating the melting temperature. PACE provides a simple interface for implementing nonlinear functions F (Eq. 8) as well as arbitrary radial functions which enables to adapt quickly to future ACE parameterizations.
We choose two distinct elements to illustrate parameterizations of ACE. Copper, for which classical potentials such as EAM are known to provide a good description of the interatomic interaction, and Si. For Si many different potentials were published to date and a recent GAP was shown to perform significantly better than other potentials 22 . We compare our Cu parameterization to a very good EAM potential and recent GTINV and SNAP potentials. The ACE for Si is compared to GAP.
For copper, EAM provides a very good description of the energy and ACE improves on this in particular for bonding environments that require angular contributions. Excellent extrapolation to new atomic environments is demonstrated, for example, for the phonons in a free standing Cu monolayer in Fig. 10. Furthermore, the longer cutoff of ACE enables us to reproduce bond breaking and making in accurate agreement with the DFT reference data. It appears that SNAP and GTINV were parameterized to selected reference data, which leads to deviations from DFT in several of our tests.
The Si ACE is comparable in accuracy to GAP, with a few key improvements. The ACE hypersurface is smoother than GAP, which is important in particular for extrapolation to large volumes as shown in Fig. 12. The improved smoothness can also be seen in the surface decohesion curve in Fig. 17 showing behavior closely matching the DFT reference. Another example of the ACE extrapolation is the fourfold defect which was highlighted in the original GAP paper, predicted erroneously to be unstable. However, ACE was fitted to the exact same DFT database and does predict a stable fourfold defect. Furthermore, it is notable that this Si ACE potential is approximately 30 times faster than the GAP of Bartók et al. 22 .

III. METHODS
We give detailed expressions for energies and forces and efficient algorithms for their evaluation in PACE in the following. For Cu and Si we employ distinct ACE forms, and different parameterization strategies follow for the two elements. The details of the parameterization strategy are provided in the SI.
A. Expressions for energy and forces

Energy
The energy of atom i is given by where F is a general non-linear function that may be supplied. Each atomic property ϕ (p) i is given by an ACE expansion, which is obtained as follows: given the relative neighbor positions r r r ji = r r r j −r r r i , r ji = |r r r ji | and directionŝ r r r ji = r r r ji /r ji , we first evaluate the atomic base where the one-particle basis φ is given in terms of spherical harmonics Y lm (r r r ji ) and radial functions R µj µi nl (r ji ) by Permutation-invariant many-body basis functions are obtained by forming products, The body order of the products is ν + 1 and the species of atom i is µ i . The vectors µ µ µ, n n n, l l l and m m m have length ν and contain atomic species, radial function indices, and spherical harmonics indices, respectively. The ACE expansion of an atomic property ϕ with expansion coefficientsc µiµnlm µnlm µnlm and lexicographically ordered indices µnlm µnlm µnlm.
The coefficientsc µµnlm µnlm µnlm are not free model parameters to be fitted since the A basis does not satisfy all required symmetries. An isometry invariant basis B is obtained by coupling elements of the A basis through the generalized Clebsch-Gordan coefficients, B = CA, which yields a linear model

Forces
The force on atom k is written as (14) and the pairwise forces f f f ki obtained using an adjoint method, 1,13 where the adjoints ω iµnlm are given by

Additional Symmetries
A straightforward opportunity for optimisation arises due to the fact that the product basis functions fulfill and t m t = 0 for rotational invariance. As we are interested in a real-valued expansion, this identity is exploited by combining the A µnlm µnlm µnlm and A µnl µnl µnl−m m m and thus reducing the computational effort for evaluating the product basis by nearly 50%. Similarly, when evaluating the forces only the real part needs to be evaluated, as the imaginary part has to add up to zero. Since ∇φ µj µinl−m (r r r ji ) = (−1) m ∇φ µj µinlm (r r r ji ) * , (21) and therefore one can limit the force evaluation to Re(ω iµ k nl0 ) Re(∇ k φ µ k µil0 (r r r ki )) Re(ω iµ k nlm ∇ k φ µ k µilm (r r r ki )) , which saves about 75% of the multiplications compared to fully evaluating all complex terms.
B. Algorithms

Model specification (PACE Input Parameters)
A PACE model is specified through four ingredients:

Compressed basis representation
To formulate the evaluation algorithms it is convenient to reorganize the basis specification into a "compressed" format. First, we enumerate the list of one-particle basis functions and the atomic base A by identifying v ≡ (µ, n, l, m), and A iv ≡ A iµnlm .
A tuple v = (v 1 , . . . , v ν ) can then be identified with µnlm and specifies a corresponding many-body basis function An atomic property ϕ i can now be written more succinctly as This format condenses the notation as well as simplifies the basis specification, now given by (i) a list of oneparticle basis functions indexed by v; and (ii) a list of many-body basis functions, each represented by a tuple v = (v 1 , . . . , v ν ) where the length ν specifies the interaction order.

Force and energy evaluation : Summary
The energy and force for a given atom i are obtained in five steps: 1. Evaluate atomic base A iµnlm , Eq.(10).
2. Evaluate product basis A iv , Eq. (12) and properties ϕ In the following we summarize algorithms for an efficient implementation.

Optimisations
Although we will not go into details of performance oriented code optimisations, we briefly mention the four most important ingredients: (1) recursive algorithms to evaluate the polynomial, radial and spherical basis sets; (2) contiguous memory layout for the many-body basis specification; (3) recursive evaluation of the many-body basis (cf. § III C); and (4) reducing the basis size and force evaluation by exploiting that the expansions are real; cf. Eqs. (20) and (23).

Stage 1: Atomic base A
First the radial functions, spherical harmonics and their respective gradients are obtained. The spherical harmonics are computed in cartesian coordinates directly 13 . Then the atomic base is evaluated. For −l ≤ m ≤ 0 we exploit Note that only for evaluations of the atomic base we need to work in µnlm notation.

Algorithm 1 Atomic base A
A iµnlm = 0 for j ← neighbors of atom i do µ = type of atom j compute R nl , dR nl , Y lm , dY lm for n, l, (m ≥ 0) do A iµnlm += R nl · Y lm end for end for for n, l, (m > 0) do For numerical efficiency all pairwise radial functions can be represented as splines with several thousands interpolation points, which makes it possible to implement arbitrary radial basis functions efficiently.  The atomic properties ϕ (p) i are computed following Eq. (13). Because A iv is used only to construct the ϕ (p) i , it need not be stored. Next, the energy E i = F(ϕ for t ← 1, . . . , ν do ωiv t += Θiv · Re(dAivt) end for end for Here, Θ iv and dA ivt are only required locally and stored in a temporary variable. The derivatives dA ivt can be computed via backward differentiation with cost that scales linearly in ν instead of the O(ν 2 ) scaling for a naive implementation. This important optimisation is implemented as follows: Algorithm 4 Compute dA ivt , t = 1, . . . , ν 1: A fwd = 1; dAiv1 = 1 2: for t ← 1, 2, . . . , ν − 1 do
Only the values A iv corresponding to interior nodes v (i.e. nodes that have at least one child), must be stored, while those corresponding to leaf nodes (i.e. nodes without any child) are only required locally to update ϕ (p) i . To evaluate the adjoints ω iv ≡ ω iµnlm we use the observation that where ∂ is a differential operator, and hence ∂F(ϕ (1) i , . . . , ϕ This shows that the adjoint of A iv can be propagated to the adjoints of the two parents. This immediately leads to the following reverse mode differentiation algorithm, which computes adjoints ω iv for all basis functions A iv (or at least those corresponding to interior nodes of the graph). However, only the adjoints for root nodes, ω iv = ω iµnlm , are eventually used to assemble the forces (Eq. (15)). The traversal of the graph must now be done in reverse order, that is, a node v may only be visited once all of its children have been visited, for example, by traversing in reverse correlation order.

Computational cost of the recursive evaluator
The forward pass, Algorithm 2R, requires 1 + P multiplications and 5 + P memory access operations at each iteration. The backward pass, Algorithm 3R, requires 2 + P multiplications and 7 + P memory access operations. In particular, the cost is (seemingly) independent of the correlation order of each basis functions. The overall cost is given by i.e., it scales linearly in the number of neighbors N and the number of nodes in the graph. The first part for setting up the atomic base is unaffected by the recursive evaluator.
Comparing the cost between the two approaches is difficult since we have no estimates on the number of artificial nodes that must be inserted into the graph. In practise we observe that there are always more leaf nodes than interior nodes, which means that relatively few artificial nodes are inserted and hence the recursive algorithm is significantly faster for large basis sets and high correlation order, but roughly comparable for small basis sets and low correlation order.
For the Cu potential with a smaller number of basis functions the timing decreases from 0.43 to 0.32 ms/atom/MD-step when the recursive evaluator is used. The Si potential, employing more basis functions, has a more pronounced speed-up from 1.84 ms/atom/MD-step to 0.80 ms/atom/MD-step.