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

Yury Lysogorskiy; Cas van der Oord; Anton Bochkarev; Sarath Menon; Matteo Rinaldi; Thomas Hammerschmidt; Matous Mrovec; Aidan Thompson; Gábor Csányi; Christoph Ortner; Ralf Drautz

Source-linked AI summary

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

Yury Lysogorskiy, Cas van der Oord, Anton Bochkarev, Sarath Menon, Matteo Rinaldi, Thomas Hammerschmidt, Matous Mrovec, Aidan Thompson, Gábor Csányi, Christoph Ortner, Ralf Drautz

arXiv:2103.00814v1Mar 1, 2021cond-mat.mtrl-sciphysics.comp-ph

TL;DR

The paper addresses the need for efficient, accurate computation of energies and forces in atomistic simulation by implementing the atomic cluster expansion in PACE. It combines efficient ACE evaluation with general-purpose Cu and Si parameterizations, obtaining a shifted accuracy–cost Pareto front and strong performance against established potentials.

Problem
Atomistic simulations require efficient energy and force computation, while existing machine-learning interatomic potentials must balance accuracy with computational cost.
Method
The paper implements ACE in the performant PACE C++ code with efficient algorithms for evaluating high body-order basis functions and integrates it with LAMMPS.
Results
PACE shifts the accuracy-versus-cost Pareto front toward higher accuracy and faster evaluation, while the Cu and Si parameterizations outperform available machine-learning potentials in reported performance, accuracy, and generalisability.
Takeaways & Limitations
PACE enables large-scale simulations with general-purpose Cu and Si potentials, including phase-melting calculations, while the Si ACE is approximately 30 times faster than GAP.

Abstract

from arXiv · show

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

PACE implements the atomic cluster expansion for efficient energy and force evaluation, connecting a complete many-body basis with practical nonlinear and linear interatomic-potential models. The paper presents Cu and Si parameterizations and evaluates their accuracy and computational cost against established potentials.

ACE represents local atomic environments using complete many-body correlation functions that can be used directly with linear regression for energies and forces.
ACE encompasses or connects to established representations and potentials, including EAM, Finnis-Sinclair, MTP, SNAP, aPIPs, symmetry functions, and SOAP.
PACE enables efficient evaluation of ACE models within the LAMMPS molecular dynamics simulation software package.
The density trick and tensor-product basis avoid the apparent O(N^ν) cost of high body-order terms.
The Cu parameterization starts from Finnis-Sinclair modeling because angular contributions are generally small in bulk copper, whereas Si requires angular interactions beyond pairwise terms.
A linear ACE parameterization for Si generalizes SOAP-GAP by including all body-order interactions up to a chosen maximum, simplifying parameterization and avoiding implicit assumptions about nonlinear terms.
The study compares its Cu and Si parameterizations with reliable literature models, including EAM, SNAP, GTINV, and GAP.

II. RESULTS AND DISCUSSION

The Cu and Si parameterizations were fitted to extensive DFT databases and evaluated across computational efficiency and bulk properties. ACE reproduces key Cu structural, thermal, and energetic behavior while PACE enables fast force evaluation.

The Cu parameterization used about 50000 DFT total-energy calculations covering clusters, bulk structures, surfaces, interfaces, point defects, and modified variants.
The Si parameterization was fitted to a database covering crystalline structures, surfaces, vacancies, interstitials, and liquid phases.
The Cu and Si fits used different strategies: nonlinear optimization for Cu and solution of a linear system for Si.
1.81 meV/atom was the ACE energy error for silicon structures within 1 eV of the ground state, compared with 1.25 meV/atom for GAP.
A single PACE force call took 0.32 ms/atom for Cu and 0.80 ms/atom for Si, supporting large-scale molecular dynamics and Monte Carlo sampling.
ACE matches the DFT structural order fcc → dhcp → hcp → bcc for low-energy Cu phases.
ACE reproduces Cu thermal expansion well up to 600 K, while all evaluated models show minor deviations at higher temperatures.

2. Interfaces and surfaces

The Cu potential was tested on interfaces, surfaces, bond breaking, point defects, clusters, and two-dimensional structures. ACE generally agrees closely with DFT and captures environment-dependent behavior that challenges simpler or shorter-range models.

Interfaces and surfaces: ACE predicts Cu stacking-fault energies in very good agreement with DFT, while GTINV gives negative values and SNAP shows larger deviations.
Interfaces and surfaces: ACE gives the best agreement with DFT for all low-index Cu surface energies, whereas SNAP and EAM consistently underestimate them.
Interfaces and surfaces: ACE is the only evaluated model that quantitatively captures how neighboring atoms alter Cu bond-breaking ranges across dimers, adatoms, and bulk slabs.
Point defects, small clusters and 2d structures: ACE reproduces Cu vacancy formation and migration energies well; other potentials overestimate the vacancy reference by 0.1-0.3 eV, while SNAP overestimates the migration barrier.
Point defects, small clusters and 2d structures: ACE and EAM agree best with recent DFT for Cu interstitials despite interstitial configurations being absent from the ACE training set.
Point defects, small clusters and 2d structures: ACE correctly predicts the Cu linear trimer's instability and the metastable bent configuration, unlike the other evaluated models.
Point defects, small clusters and 2d structures: Only ACE and GTINV identify the planar equilateral rhombus as the Cu tetramer's ground state; EAM and SNAP incorrectly favor a close-packed tetrahedron.
Point defects, small clusters and 2d structures: DFT and ACE predict real phonon frequencies for the 2D hcp Cu lattice, whereas EAM and SNAP show out-of-plane dynamic instabilities.

D. Silicon

The silicon ACE potential matches DFT across a broad set of bulk and thermal properties, while improving extrapolation to unseen volumes and remaining comparable in accuracy to GAP. It also reproduces liquid and amorphous structural behavior against reference data.

Scope: The silicon ACE was fit to the extensive database previously used for GAP and evaluated on bulk, surface, liquid, amorphous, and random-structure properties.The benchmark spans crystalline structures, surfaces, vacancies, interstitials, liquid phases, and a random structure search.
Unseen structures: For the hcp′ structure, both ACE and GAP predict the minimum despite its absence from the DFT database; GAP matches the reference energy better, while ACE estimates curvature better.The hcp′ structure has c/a < 1 and was not contained in the reference silicon database.
Extrapolation: ACE extrapolates better than GAP at large volumes, avoiding GAP’s unphysical high-energy minima and remaining close to the DFT energy-volume reference.Both potentials describe the minima near 15 and 20 Å3/atom for bcc and diamond silicon.
Bulk properties: Both ACE and GAP match DFT elastic constants within a few percent and accurately reproduce the silicon phonon spectrum.GAP has a slightly better match to the DFT phonon bandwidth.
Thermal properties: ACE models silicon’s negative low-temperature thermal expansion very well and gives almost perfect agreement for heat capacity, although high-temperature expansion saturates too high.The high-temperature saturation bias is reported for both ACE and GAP.

2. Surfaces

The silicon and copper potentials are tested on surfaces, decohesion, point defects, and random structures. ACE agrees closely with DFT, provides smoother or more reliable behavior in challenging configurations, and demonstrates strong extrapolation beyond the fitting data.

Surface energies: ACE and GAP agree very well with DFT for silicon surface formation energies in the (100), (110), and (111) directions.
Surface decohesion: ACE’s surface-decohesion energy curve is significantly smoother than GAP’s and closer to the DFT reference along the bulk-to-unrelaxed-(110)-surface path.Bulk diamond and the unrelaxed surface were included in the database, along with some configurations along the path.
Point defects: Both ACE and GAP predict silicon point-defect formation energies very well, including vacancy and tetragonal, hexagonal, and dumbbell interstitial defects.
Point defects: ACE stabilizes the fourfold-coordinated defect and gives a significantly better relaxed-defect energy than GAP, despite no nearby defect or transition-state configurations in the fitting database.ACE and GAP make similar errors near the transition state, but ACE better represents the relaxed defect structure.
Random structure search: In random structure search, ACE produces an energy-volume distribution and density of states similar to DFT, while empirical potentials tested previously failed completely.The test uses relaxed eight-atom configurations initialized with interatomic distances greater than 1.7 Å.
Discussion: PACE implements ACE in LAMMPS and shifts the accuracy-versus-cost Pareto front toward faster and more accurate machine-learning potentials.The implementation is presented as enabling efficient evaluation within large-scale molecular-dynamics simulations.
Copper: For copper, ACE improves on EAM for bonding environments requiring angular contributions and reproduces bond breaking and making more accurately than selected SNAP and GTINV comparisons.The longer ACE cutoff supports the bond-breaking and bond-making description relative to DFT.
Silicon: For silicon, ACE is comparable in accuracy to GAP, with smoother behavior at large volumes and surface decohesion, stable fourfold defects, and approximately 30-times-faster evaluation.

1. Energy

PACE evaluates ACE atomic properties, energies, and forces through symmetry-aware many-body basis functions and efficient algorithms. The formulation supports linear ACE models as well as supplied nonlinear embedding functions.

Energy model: The energy uses a general nonlinear function F of ACE atomic properties, allowing nonlinear embedding functions alongside the linear ACE expansion.The paper gives distinct ACE forms and parameterization strategies for copper and silicon.
Atomic basis: PACE represents each atomic property with an ACE expansion built from neighbor positions, radial functions, spherical harmonics, and permutation-invariant products.The resulting product basis has body order ν + 1 and is indexed by species, radial-function, and angular-momentum labels.
Symmetry: ACE obtains an isometry-invariant basis by coupling the atomic basis through generalized Clebsch–Gordan coefficients, with free coefficients optimized during fitting.The constraints enforce invariance under rotation and inversion.
Forces: PACE computes forces from the energy using pairwise force contributions evaluated with an adjoint method.The paper provides explicit expressions for the force on atom k and the pairwise forces.

3. Additional Symmetries

Rotational invariance and the real-valued nature of the expansion provide symmetry reductions for product-basis and force evaluations. These reductions substantially lower computational effort.

Additional Symmetries: Combining coefficients related by m and −m reduces product-basis evaluation effort by nearly 50%.This exploits the identity t_mt = 0 required for rotational invariance.
Additional Symmetries: Force evaluation needs only the real part because the imaginary part sums to zero.The resulting restriction removes unnecessary complex-term evaluations.
Additional Symmetries: 75% of the multiplications are saved compared with fully evaluating all complex terms.This reduction applies to the force evaluation.

B. Algorithms

PACE models are specified by radial and many-body basis information together with expansion coefficients, then evaluated through a five-step energy-and-force pipeline. A compressed basis representation organizes the implementation around interaction-order tuples.

B. Algorithms: A PACE model uses a radial basis, a list of basis functions for each order, and corresponding expansion coefficients.The radial basis may use splines or polynomial recursion.
B. Algorithms: The compressed specification indexes one-particle basis functions by v and represents many-body functions as tuples v = (v1, . . . , vν).The tuple length ν specifies interaction order.
B. Algorithms: The compressed format condenses notation and simplifies specification into one-particle and many-body basis-function lists.Many-body functions are represented by interaction-order tuples.
B. Algorithms: The algorithms are designed for efficient implementation after the basis specification has been reorganized into compressed form.The paper introduces the algorithms following the five-step evaluation sequence.
B. Algorithms: Energy and force evaluation proceeds through five stages from atomic-base evaluation to energy derivatives, adjoints, product-basis derivatives, and force assembly.The listed steps include evaluating A, obtaining E_i and derivatives, computing adjoints and derivatives, then assembling f_ji.

4. Optimisations

PACE optimizations combine recursive basis evaluation, contiguous basis storage, and real-valued symmetry reductions. The implementation also uses splines for efficient radial-function evaluation and avoids storing temporary product-basis values.

4. Optimisations: Performance relies on recursive basis evaluation, contiguous memory layout, recursive many-body evaluation, and real-expansion reductions.These are identified as the four most important performance-oriented ingredients.
4. Optimisations: Radial functions, spherical harmonics, and gradients are computed before evaluating the atomic base.Spherical harmonics are computed directly in Cartesian coordinates.
4. Optimisations: Algorithm 1 evaluates atomic-base contributions from radial functions and spherical harmonics, including nonnegative m terms.The atom type μ indexes the atomic-base contributions.
4. Optimisations: Representing pairwise radial functions as splines with thousands of interpolation points enables efficient use of arbitrary radial basis functions.The spline representation is used for numerical efficiency.
4. Optimisations: Product-basis functions and their derivatives are constructed after the atomic base, with temporary product-basis values omitted from storage.The product basis contributes to the properties used to obtain the energy.
4. Optimisations: The energy is evaluated as a function of the constructed properties after product-basis contributions are accumulated.The implementation then proceeds to the corresponding derivatives and adjoints.

7. Adjoints ω

Adjoint and force evaluation is optimized through local temporary quantities, backward differentiation, and real-valued force expressions. The resulting computational cost scales linearly with key basis and neighbor counts.

7. Adjoints ω: Adjoints are computed after the energy derivatives, using locally stored temporary quantities for the required intermediate values.The adjoints ω_iμnlm follow Eq. (17).
7. Adjoints ω: Backward differentiation computes product-basis derivatives with cost linear in ν instead of O(ν^2) for a naive implementation.Forward and backward passes provide the derivative factors.
7. Adjoints ω: 3ν − 5 multiplications are required for ν ≥ 2 after removing multiplications by one.This is the optimized multiplication count for computing dA.
7. Adjoints ω: Force gradients are obtained from the optimized derivative expressions and assembled over neighboring atoms.Algorithm 5 initializes forces and accumulates neighbor contributions.
7. Adjoints ω: The overall cost scales linearly with neighbor count N, maximum correlation order νmax, and the number of properties ϕ(p).Atomic-base evaluation contributes O(N · #A), while correlation evaluation contributes O((νmax + #ϕ) · #A).
7. Adjoints ω: An alternative evaluation scheme is introduced to further reduce dependence on ν.The paper states this reduction follows the initial cost analysis.

C. Recursive evaluator

The recursive evaluator represents higher-order ACE basis functions as products of lower-order functions arranged in a directed acyclic graph. Forward and reverse traversals reduce arithmetic costs while accommodating basis-function constraints through zero-coefficient artificial nodes.

C. Recursive evaluator: The recursive evaluator reduces higher-order basis-function evaluation to products of two lower-correlation-order basis functions.This replaces the conventional sequence of ν − 1 products with a single product for each decomposed basis function.
C. Recursive evaluator: Basis-function constraints can prevent valid decompositions, so the method inserts artificial basis functions with zero coefficients.The resulting construction is a directed acyclic graph whose nodes represent basis functions.
C. Recursive evaluator: The graph is built from atomic-base root nodes and recursively inserted nodes whose parent decompositions are already present.Evaluation proceeds in increasing correlation order so each node’s parents are available.
C. Recursive evaluator: The heuristic graph construction achieves excellent performance, although further optimization could reduce artificial nodes and memory access.Its effectiveness depends on the graph structure and the number of inserted auxiliary nodes.
C. Recursive evaluator: The forward algorithm computes basis-function values and property accumulations by traversing graph decompositions.Only interior-node values need persistent storage; leaf values are used locally.
C. Recursive evaluator: Reverse-mode differentiation propagates adjoints from each graph node to its two parents, with root adjoints used to assemble forces.The graph must be traversed in reverse order for this propagation.

3. Computational cost of the recursive evaluator

The recursive evaluator’s cost is linear in neighbor count and graph size rather than directly in basis-function correlation order. In practice, it is substantially faster for large basis sets and high correlation order, with measured speedups for both Cu and Si.

3. Computational cost of the recursive evaluator: The forward and backward passes require fixed per-node operation counts, making their cost seemingly independent of each basis function’s correlation order.The forward pass requires 1 + P multiplications and 5 + P memory accesses; the backward pass requires 2 + P and 7 + P, respectively.
3. Computational cost of the recursive evaluator: O(N · #A) + O(#G · P) is the overall cost, scaling linearly with neighbors N and graph nodes #G.The atomic-base setup contributes the first term and is unaffected by recursive evaluation.
3. Computational cost of the recursive evaluator: 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.This behavior reflects the observed predominance of leaf nodes over interior nodes and the resulting limited number of artificial nodes.
3. Computational cost of the recursive evaluator: 0.43 to 0.32 ms/atom/MD-step is the Cu timing change when the recursive evaluator is used.The Cu potential uses a smaller number of basis functions than the Si potential.
3. Computational cost of the recursive evaluator: 1.84 ms/atom/MD-step to 0.80 ms/atom/MD-step is the Si timing change with the recursive evaluator.The more pronounced Si speed-up occurs for a potential employing more basis functions.