Source-linked AI summary
Machine learning a general purpose interatomic potential for silicon
Albert P. Bartok, James Kermode, Noam Bernstein, Gabor Csanyi
TL;DR
Materials simulations need interatomic potentials that avoid first-principles computational limits while remaining accurate across relevant configurations. This paper develops a Gaussian Approximation Potential for silicon by fitting a broad DFT-based dataset, and finds high accuracy across diverse structures and properties, with uncertainty estimates indicating reduced reliability in extrapolative cases. The model is nevertheless scoped mainly to ambient-pressure elemental silicon.
Problem
Computational cost limits first-principles simulation scales, while flexible machine-learning potentials require sufficiently broad yet manageable training databases to remain accurate across relevant configurations.
Method
The authors construct a silicon GAP by fitting a broadly chosen dataset and protocol to DFT energies, forces, and virial stresses, treating the model as an interpolation scheme for DFT.
Results
GAP achieves uniformly high accuracy across a wide range of silicon properties and systems, including bulk structures and point and plane defects, while planar-defect tests remain within 20% of DFT.
Takeaways & Limitations
The study demonstrates that a useful machine-learning interatomic potential can comprehensively describe a material while retaining locality and linear-scaling computational cost.
Takeaways & Limitations
The potential focuses on ambient pressure and elemental silicon, omitting comprehensive high-pressure phases and multielement interactions.
Abstract
from arXiv · showhide
The success of first principles electronic structure calculation for predictive modeling in chemistry, solid state physics, and materials science is constrained by the limitations on simulated length and time scales due to computational cost and its scaling. Techniques based on machine learning ideas for interpolating the Born-Oppenheimer potential energy surface without explicitly describing electrons have recently shown great promise, but accurately and efficiently fitting the physically relevant space of configurations has remained a challenging goal. Here we present a Gaussian Approximation Potential for silicon that achieves this milestone, accurately reproducing density functional theory reference results for a wide range of observable properties, including crystal, liquid, and amorphous bulk phases, as well as point, line, and plane defects. We demonstrate that this new potential enables calculations that would be extremely expensive with a first principles electronic structure method, such as finite temperature phase boundary lines, self-diffusivity in the liquid, formation of the amorphous by slow quench, and dynamic brittle fracture. We show that the uncertainty quantification inherent to the Gaussian process regression framework gives a qualitative estimate of the potential's accuracy for a given atomic configuration. The success of this model shows that it is indeed possible to create a useful machine-learning-based interatomic potential that comprehensively describes a material, and serves as a template for the development of such models in the future.
I. INTRODUCTION
First-principles simulations are limited by computational cost, while conventional potentials struggle to reproduce many properties simultaneously. Machine learning offers flexible high-dimensional fitting, but useful interatomic potentials still require manageable databases covering relevant configurations without sacrificing transferability.
- DFT-based simulations become impractical for thousands of atoms or millions of energy and force evaluations because of computational expense and unfavorable scaling.
- Conventional analytical potentials use simple functional forms that often cannot match many macroscopic and microscopic properties simultaneously.
- Machine learning avoids simply parameterized functional forms, enabling flexible fitting in high-dimensional configuration spaces.
- Machine-learning models are typically accurate near their training data but progressively less reliable for configurations farther away, limiting transferability.
- Accurate machine-learning potentials require tractable databases that adequately sample configurations relevant to a wide range of atomistic simulations.
C. A general potential
The study develops and validates a silicon machine-learning potential against DFT across diverse structures and properties. GAP reproduces bulk, defect, and other tests accurately, while its error estimates identify less reliable extrapolative configurations.
- C. A general potential: The model targets physically important silicon structures and properties spanning bulk phases, surfaces, point and planar defects, and configurations relevant to materials modelling.
- C. A general potential: The silicon GAP is designed as a broadly applicable interpolation of DFT, using a dataset and protocol intended to be minimally specific to material and observable.
- C. A general potential: Comparisons with published analytical potentials are not definitive evaluations because those models were fitted to different property sets and reference sources.
- C. A general potential: GAP reproduces bulk moduli and diamond-cubic elastic constants with fractional errors below 10% relative to DFT.
- C. A general potential: Planar-defect tests outside the training set remain within 20% of DFT, while reduced confidence near these defects corresponds to larger prediction errors.
II. METHODOLOGY
The silicon GAP combines a repulsive pair potential with a SOAP-based many-body Gaussian-process model fitted to DFT energies, forces, and stresses.
- Potential construction: The total energy combines a predefined pair potential with a many-body kernel expansion over representative atomic environments.The pair term captures short-range repulsion, while the many-body fit describes the attractive contribution.
- SOAP representation: The SOAP kernel measures similarity between local neighbourhoods using rotationally and permutationally invariant environment descriptions.Neighbour densities are expanded in radial functions and spherical harmonics, and rotational invariance is constructed through integration over SO(3).
- Fitting procedure: The model coefficients are obtained from a linear system using DFT total energies, forces, and virial stress components.Regularisation weights are interpreted through expected errors in the fitted quantities, including model error from the finite cutoff.
- Fitting procedure: Representative environments are selected with CUR decomposition to construct a low-rank approximation of the full kernel matrix.The selection aims to retain diverse environments while limiting the number of basis environments.
- Uncertainty estimation: Gaussian-process regression supplies a predicted variance interpreted as a one-sigma error bar for atomic energies.The simplified estimate uses σe, typically set to 1 meV, as the regularisation parameter for the uncertainty calculation.
B. Database
The silicon database was assembled from diverse periodic-cell configurations using DFT references, and the resulting GAP achieved substantially lower force errors than the analytical potentials tested.
- Database: The database was assembled over an extended period from configuration types selected to provide coverage for a range of properties.Its final size reflects practical choices influenced partly by the computational cost of electronic-structure calculations.
- Force errors: GAP had a median force-component error of about 0.025 eV/Å, an order of magnitude smaller than the analytical potentials.This comparison was made relative to DFT calculations on the fitting database.
- Force errors: The separate testing database contained a grain boundary, six di-interstitials, stacking-fault paths, and an amorphous configuration.Its force-error distribution was very similar to the fitting-database distribution, although the actual errors depended strongly on geometry proportions.
- Testing design: The testing database represented extrapolation to configurations entirely different from those in training rather than a usual random train-test split.This made it a more stringent test of transfer beyond the fitting configurations.
C. Convergence
SOAP/GAP converges with basis-set size and descriptor resolution, reaching meV-scale energy accuracy, while finite cutoff and database coverage remain important limits.
- Convergence: The SOAP/GAP convergence study varied the number of basis functions and the SOAP descriptor expansion controlling atomic-density resolution.These parameters directly trade computational cost against fitting accuracy.
- Convergence: The model achieved energy-surface precision of order 1 meV for β-Sn and about an order of magnitude better for diamond.The diamond and β-Sn errors were evaluated using the Δ-value against DFT energy-volume curves.
- Limitations: With a 5 Å cutoff, previous investigations indicate an expected diamond force error of about 0.1 eV/Å from locality alone.The finite cutoff means the DFT force is not strictly determined by the model’s finite neighbourhood.
- Limitations: The database is not mathematically complete, so errors can improve for a given application through a larger cutoff only if the database is enlarged accordingly.The practical benefit depends on whether the application’s relevant configurations are represented in the database.
D. Testing
An ASE-based testing framework runs each model and test as an independent Python module against consistent reference DFT calculations.
- Testing framework: Each model and test is implemented as an independent Python module using Atomic Simulation Environment interfaces to QUIP.Reference DFT results use the same tests with an ASE interface to Castep.
- Testing framework: The automated framework keeps starting configurations, minimisation algorithms, and reported test results consistent across models.It also supports rerunning tests after changes to the GAP training database.
III. RESULTS: VERIFICATION
Verification tests assess configurations deliberately included in the fitting database. GAP agrees excellently with DFT across tested crystal structures, while the out-of-database hcp’ structure shows the largest discrepancy.
- Verification tests configurations selected for inclusion in the database to describe the corresponding observables.The authors distinguish verification from validation because the database was designed around these tests.
- The hcp’ structure, absent from the fitting database, has a larger DFT–GAP discrepancy than the other tested structures but remains in good agreement.
- GAP reproduces DFT energy–volume curves across tested crystal structures, including equilibrium volumes, cohesive-energy depths, and bulk-modulus curvatures.
- Only GAP qualitatively reproduces the diamond, β-Sn, and fcc crystal structures together; many other models fail to reproduce β-Sn.
B. Liquid
The liquid and amorphous tests compare GAP with DFT, experiment, and other interatomic models using structural, dynamical, and energetic observables. GAP closely matches reference behavior and produces comparatively stable amorphous silicon structures.
- Liquid: GAP liquid RDFs match DFT peak heights and positions, while its ADF reproduces the narrow 60° and broader 100° peaks.The ADF is described as a more stringent test than the RDF; most other potentials distort the relative peak heights.
- Liquid: Liquid diffusivity was evaluated with 512-atom variable-cell constant-enthalpy simulations from approximately 1700 K to 2200 K against experiment, DFT, and SW-potential results.
- Amorphous phase: Amorphous silicon was generated by cooling a 216-atom liquid sample from 2000 K to 500 K at 10^12 K/s, then relaxing positions and cell degrees of freedom.
- Amorphous phase: Some potentials start from strongly undercooled liquids because they overestimate the melting temperature, limiting direct comparability of their quenched structures with experiment.
- Amorphous phase: GAP, EDIP, and Tersoff best match experimental amorphous RDFs, while GAP and Tersoff have the lowest coordination-defect concentrations, closest to experimental estimates of ≤1%.
- Amorphous phase: GAP gives the closest amorphous–crystal energy difference to experiment, and DFT evaluation shows its quenched structure has the lowest relative energy among compared potentials.DFT relaxation produces only a small energy reduction for GAP, EDIP, and Tersoff, unlike several other potentials.
D. Phase diagram
Phase behavior provides a stringent benchmark because finite-temperature transitions probe the potential energy surface through a balance of energetic and entropic effects. GAP phase calculations used coexistence simulations for several silicon phases and pressures.
- Benchmark and method: Phase behavior tests both model realism and microscopic details of the potential energy surface.Finite-temperature transitions probe relatively high-energy configurations through a delicate energetic–entropic balance.
- Benchmark and method: GAP phase-boundary calculations used coexistence simulations for diamond/liquid and simple-hexagonal/liquid systems at pressures from 0 to 12 GPa.The simulations contained 432 atoms for diamond/liquid and 1024 atoms for simple-hexagonal/liquid systems.
- Additional validation: Point-defect formation energies for vacancies and three interstitial positions were within at most 7% of DFT reference values.The hexagonal, tetrahedral, and dumbbell interstitials were included alongside vacancies in the comparison.
- Additional validation: GAP preserved the DFT ordering of hexagonal and tetrahedral interstitial energies, whereas many other potentials made much larger errors or reversed the ordering.The energy difference between the two interstitial geometries was underestimated by GAP.
2. Surfaces
Surface tests examine bond breaking, reconstruction, and surface energetics across silicon orientations. GAP improves agreement with DFT by representing separation-path configurations and captures key low-energy reconstructions more accurately than analytical potentials.
- Surface significance: Surface properties matter for fracture, environmental reactions, functional layers, and crystal growth because surfaces expose complex bonding environments.Surface bonding reflects an interplay between strain effects and termination-specific structure.
- Decohesion: Including configurations along the separation path corrected GAP’s intermediate decohesion energies while retaining the correct fully separated-surface energy.An earlier GAP trained only on separated surfaces reproduced the endpoint but not the path.
- (100) reconstruction: GAP reproduced the tilted-dimer (100) reconstruction within -2.5° of DFT, whereas analytical potentials showed zero tilt except EDIP at 4°.The DFT dimer tilt is 18° and is associated with a Jahn–Teller effect.
- (111) reconstructions: GAP predicted the (111) DAS reconstruction family below the unreconstructed surface with energy errors below 0.05 J/m2 despite training on one 3 × 3 DAS configuration.Representative environments were selected from atoms unique to the DAS reconstruction family.
- (111) reconstructions: The present GAP minimum occurs at n = 5, within 0.01 J/m2 of the 7 × 7 structure, so the finest DAS-family energy ordering is not robust.More data relevant to these structures would be needed to resolve the smallest relative energies reliably.
F. Crack Propagation
Crack propagation is difficult because realistic fracture requires accurate bond breaking in very large systems. GAP was tested on silicon cleavage and on broader configuration-space probes, reproducing key fracture behavior while revealing imperfect local barriers.
- Motivation and method: Accurate crack-tip bond breaking must be combined with large systems to avoid unrealistic strain gradients, limiting hybrid quantum/classical approaches to DFT-accessible timescales.This makes thermally activated crack-growth studies especially challenging.
- Fracture test: GAP fracture simulations used a 23,496-atom Si(111)[1¯10] cleavage system with molecular dynamics and quasi-static loading.Molecular dynamics was run at 300 K over strain rates from 10^-6 to 10^-4 fs^-1.
- Fracture result: GAP predicted brittle fracture with atomically smooth surfaces and occasional crack-tip reconstructions followed by downward [111] surface steps.The crack-tip reconstruction appeared in a quasi-static simulation at G = 5.13 J/m2.
- Random structure search: Random structure search tested GAP over distorted, high-energy geometries beyond physically reasonable minima, and its energy–volume distribution resembled DFT more closely than those of other potentials.The search used constrained 8-atom cells relaxed with a two-point steepest-descent method.
- Random structure search: DFT relaxation showed that 80% of GAP-minimum volume changes were below 1 Å3/atom, while some GAP minima differed in positions and small barrier heights.Along representative trajectories, DFT and GAP energies remained within 0.05 eV/atom.
B. Phonons and thermal expansion
GAP was evaluated for vibrational, thermal, and stacking-fault properties that are not all explicitly represented in the fitting database. It agrees well with DFT for phonons and several defect energies, but anharmonic quantities remain less resolved.
- Phonons: Phonons and mode Grüneisen parameters probe the potential near diamond-structure minima and influence thermodynamic and transport behavior.These quantities were compared across GAP, DFT, DFTB, and analytical potentials.
- Phonons: GAP phonon frequencies agree excellently with DFT even though frequencies were not explicitly included in the database.Analytical potentials generally overestimate acoustic zone-edge and optical-branch frequencies.
- Anharmonicity: GAP mode Grüneisen parameters are mixed, with a large transverse-acoustic discrepancy indicating insufficient near-equilibrium force data for resolving anharmonicity.This limitation affects the model’s description of the anharmonic potential-energy surface.
- Thermal expansion: GAP thermal-property calculations tested thermal expansion, Grüneisen parameters, and heat capacity using the quasi-harmonic approximation and phonon frequencies.The benchmark included silicon’s negative low-temperature thermal expansion.
- Planar defects: GAP generalized stacking-fault energies showed reasonable agreement with DFT, comparable to the best other interatomic potentials and better than most.The stable stacking-fault energy γsf was positive but underestimated, and its predicted error suggested database extension could improve it.
D. Grain boundary
The silicon GAP is tested on grain boundaries and multiple defect processes, generally matching DFT while exposing specific transferability failures and uncertainty signals.
- Grain boundary: The (112) Σ3 grain boundary was evaluated with DFT and several interatomic potentials despite being absent from the fitting database.Its energy per unit area and relaxed geometry were compared across models.
- Grain boundary: GAP forces and relaxed geometries show small differences from DFT for the tested planar defect, while other potentials vary from close agreement to large energy errors.The DFT-relaxed geometry has small GAP forces, and corresponding displacements are nearly imperceptible.
- Point defects: The GAP does not stabilize silicon’s four-fold coordinated defect because the bond-rotation process and final structure are absent from its database.GAP and DFT energies agree while the predicted error remains below about 3 meV/atom, then strongly deviate.
- Vacancy migration: GAP produces the most accurate vacancy migration minimum-energy path relative to DFT, although it underestimates the barrier.The comparison uses NEB paths for 63-atom diamond cells.
- Vacancy migration: For vacancy migration, predicted error is low where agreement with DFT is good and rises as agreement decreases beyond 3 meV per atom.This makes the uncertainty estimate a qualitative guide to model reliability along the path.
- Di-interstitials: All tested potentials struggle with the TT di-interstitial configuration, with even the best error exceeding 10%, while GAP otherwise stays below about 6%.Stillinger–Weber has competitive energy accuracy but substantially more distorted relaxed geometries than GAP.
V. CONCLUSION
The silicon GAP matches a broad range of target potential-energy-surface properties while retaining locality and linear-scaling cost, but remains limited in scope, database construction, and computational efficiency.
- Conclusion: The non-parametric approach yields uniformly high accuracy across bulk structures, point defects, and plane defects while retaining locality and linear-scaling computational cost.The paper emphasizes accurate energies and forces for fully bonded and bond-breaking configurations.
- Conclusion: Accurate treatment of bond breaking enables reliable descriptions of surface decohesion, unstable stacking-fault paths, and point-defect migration barriers.These processes are identified as especially challenging for conventional interatomic potentials.
- Limitations: The potential was focused on ambient pressure and elemental silicon, excluding comprehensive high-pressure phases, multielement systems, and long-range interactions beyond 5 Å.Including long-range interactions with high-dimensional short-range fits remains an outstanding problem.
- Limitations: The training database was assembled by hand through an ad hoc iterative process, motivating protocols for more automated database construction.The paper suggests using built-in uncertainty quantification to select novel configurations for inclusion.
- Limitations: The implementation costs around 100 ms/atom, and the authors do not claim that the SOAP kernel and basis functions are optimal.Cheaper or better-suited basis functions might achieve the same accuracy.
- Outlook: The potential may be extended with targeted geometries or more accurate reference data, and the authors propose it as a template for models of other materials.The potential is distributed for use with QUIP and LAMMPS.