Efficient and accurate spatial mixing of machine learned interatomic potentials for materials science

TL;DR

This paper introduces ML-MIX, a tool that combines different interatomic potentials to speed up materials simulations without losing accuracy.

Contribution

ML-MIX enables spatial mixing of machine-learned interatomic potentials for efficient large-scale simulations.

Findings

Speedups of up to 11× were achieved in simulations with ~8000 atoms.

He reflection coefficients matched experimental results up to 80 eV for the first time.

The method works with ACE, UF3, SNAP, and MACE potential architectures.

Abstract

Machine-learned interatomic potentials can offer near first-principles accuracy but are computationally expensive, limiting their application to large-scale molecular dynamics simulations. Inspired by quantum mechanics/molecular mechanics methods, we present ML-MIX, a CPU- and GPU-compatible LAMMPS package to accelerate simulations by spatially mixing interatomic potentials of different complexities, allowing deployment of modern MLIPs even under restricted computational budgets. We demonstrate our method for ACE, UF3, SNAP and MACE potential architectures and demonstrate how linear ‘cheap’ potentials can be distilled from a given ‘expensive’ potential, allowing close matching in relevant regions of configuration space. The functionality of ML-MIX is demonstrated through tests on point defects in Si, Fe and W-He, in which speedups of up to 11× over ~8000 atoms are demonstrated, without sacrificing accuracy. The scientific potential of ML-MIX is demonstrated via two case studies in W, measuring the mobility of \documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$b=\frac{1}{2}\langle 111\rangle$$\end{document}b=12〈111〉 screw dislocations with ACE/ACE mixing and the implantation of He with MACE/SNAP mixing. The latter returns He reflection coefficients which (for the first time) match experimental observations up to an He incident energy of 80 eV—demonstrating the benefits of deploying state-of-the-art models on large, realistic systems.

Figures

Demonstration of convergence of key quantities with hard constraint strength for the Si ACE 2, 10 cheap model.Error is shown compared to the Si ACE 4, 20 expensive model. In each case, the unconstrained case is shown by a dashed horizontal red line and the result at each constraint strength (*α*) is shown by a blue point. a The L2 norm of the error over C_11_, C_12_ and C_44_. b The absolute error in lattice parameter. Note that at a constraint strength of 10^−7^, the lattice parameter error was 0 up to the number of decimal places measured (5).

Region construction schematic.From left to right: starting with a cell containing an initial defect, it is identified (orange atoms in panel two) either using a user-specified LAMMPS group or by the output vector from a user-defined LAMMPS fix. The defect atoms are then iterated over, with atoms that sit less than one core radius out from any defect atom being added to the core region. This construction process is then repeated for the atoms in the blended region.

An illustration of force-mixing with interatomic potentials.A mixed domain is split into two non-interacting subdomains (blue and orange), which are composed of all the atoms in the core and blending regions for each potential (solid line), plus an extra buffer region containing atoms up to the potential cutoff radius out (dotted line). Partial force vectors F_c_ and F_e_ are obtained by evaluating the cheap and expensive subdomains respectively. These partial vectors are then combined into the full force vector on the system F. In this process, forces on buffer atoms are discarded and forces

A schematic overview of the simulation acceleration workflow.Left side: the constrained potential fitting process, whereby an accurate expensive potential (orange) is approximated in local regions of potential energy space by a simpler cheap potential (blue). Right side: schematic showing how simulation is accelerated; cheap potential is used to evaluate force components on atoms in `bulk-like' environments.

The energy barrier associated with vacancy migration in Si, obtained through all-expensive (blue curve) and ML/ML (orange curve) nudged elastic band simulations.For the ML/ML simulation, the final energies were obtained through one-shot energy evaluations of the relaxed structures with the expensive potential. The all-expensive energy barrier (0.0543 eV) agrees with the ML/ML energy barrier (0.0548 eV) within 1 meV. The tangent lines arise from the projected forces.

Average force on stretched bond in silicon over 100 ps at 300 K for different expensive potential radii $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}$$\end{document}$ (blue points).Error bars represent the standard error in bond force measured over this time, with samples taken every 15 fs. For $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{a

Diffusion coefficient as a function of temperature for the diffusion of the iron interstitial dumbbell defect through bulk iron and for helium in bulk tungsten.a The diffusion coefficient of the iron dumbbell for three cases: expensive reference potential only (blue), cheap potential only (orange) and mixed simulation (green). The mixed simulation leads to results that are well matched with the reference. Each point is the average of five diffusion coefficients measured in independent runs of 1 ns each. b The diffusion coefficient for He in W using a helium-tungsten ACE (blue) and a mixed simu

Simulation set-up and results for thermally activated glide of screw dislocations in W at different shear stresses and temperatures.a Simulation set-up, where dimensions in *x*, *y* and *z* are 223Å × 221Å × 68Å(25*b*), respectively, for a total of 195,000 atoms. Stress was applied to the free surfaces at the top and bottom of the cell in 〈111〉 type directions. The expensive atoms around the dislocation (as detected using common neighbor analysis in LAMMPS) are shown in a column, where the colors correspond to the core (red), blending (green) and buffer (white) regions. b The velocity results

Reflection coefficient as a function of incident He energy in eV for normal deposition onto a 1000 K W {100} W surface.Experimental results from Van Gorkum et al. are shown in green^[60](#CR60)^. Results from the all-expensive potential (fine-tuned mpa-0-medium model) are displayed in blue, and results obtained with ML-MIX are shown in orange. Results from a previous MD simulation using a W EAM and W-He pair potential are displayed in purple^[59](#CR59)^. For the all expensive, to limit computational cost, only energies up to 40 eV were tested with 100 repeats each. With ML-MIX, energies up to

Demonstration of convergence of key quantities with hard constraint strength for the Si ACE 2, 10 cheap model.Error is shown compared to the Si ACE 4, 20 expensive model. In each case, the unconstrained case is shown by a dashed horizontal red line and the result at each constraint strength (*α*) is shown by a blue point. a The L2 norm of the error over C_11_, C_12_ and C_44_. b The absolute error in lattice parameter. Note that at a constraint strength of 10^−7^, the lattice parameter error was 0 up to the number of decimal places measured (5).

A schematic representation of the constrained fitting of a cheap potential to an expensive one.a A representation of the local potential energy surface of the expensive reference model (orange curve). Data has been sampled from this potential energy surface (crosses). b The cheap model, fit to this data in order to locally match the expensive model. It approximately matches the `soft' constraint data (pink crosses) and exactly matches the `hard' constraint data (green crosses).

Automated selection of lambda value during a single constrained fit to satisfy a given alpha.The vertical axis shows the energy RMSE per atom. The target alpha value is shown as a dashed red horizontal line. Each blue point represents a step on the search.

Full text

Introduction

Atomistic simulations face two central challenges: modeling quantum-mechanical atomic energies and forces and capturing physical processes at large length and time scales. Addressing these challenges under a fixed computational budget leads to cost-accuracy trade-offs, which have driven methodological development over many decades. The most important (and most successful) approximation for electronic interactions is density functional theory (DFT)—a mean field approach which increases the size of simulations from one or two small atoms to hundreds or even thousands on modern machines^1,2^.

A further level of approximation is to treat electrons implicitly, predicting atomic energies and forces only from nuclear positions and chemical species. For many years, the most accurate models in this category were interatomic potentials with simple functional forms fit to match experimental quantities of interest^3–6^ or DFT-derived properties^7,8^. Many such ‘empirical’ potentials still see considerable use today^9^. In the last three decades a new class of interatomic potentials has emerged, which have complex functional forms and parameters fit to DFT data, with an aim to reproduce large swathes of the potential energy surface^10–15^. These machine learned interatomic potentials (MLIPs) have achieved remarkable accuracy over a range of applications across materials science and chemistry^16–19^. Whilst the evaluation costs of empirical potentials and MLIPs both scale linearly with the number of atoms in the system, the complex functional forms of MLIPs result in much larger scaling pre-factors; MLIPs are typically 2–3 orders of magnitude slower in evaluation time than their simple empirical counterparts^20,21^.

This cost has slowed the uptake of MLIPs in fields where simulation domains are large and/or timescales are long, which encompasses a large range of phenomena in materials science. It is useful to split such simulations into two categories, those where the entire domain is chemically complex, with atoms involved in bond-making/breaking processes throughout (e.g., amorphous materials^22,23^, liquids^24,25^, high entropy alloys^26,27^), and those where the domain contains only isolated regions of complexity surrounded by atoms in simple environments (e.g., materials defects^28^, biochemical simulations^29^, catalysis^30^). In the second case, simulation cost can in principle be reduced through spatial decomposition: using cheaper, approximate models away from the chemically ‘important’ regions of high complexity.

Spatial decomposition has been extensively explored in the context of quantum mechanical/ molecular mechanics modeling (QM/MM), a method which couples an expensive first-principles method (typically DFT) with a linear scaling potential for the complex (QM) and simple (MM) regions of a simulation domain, respectively^31–39^. In this paper, we detail a method that is conceptually similar to QM/MM, but instead combines MLIPs of different complexities (ML/ML). The key difference is that ML/ML simulations aim to run for millions of timesteps, compared to only hundreds or thousands in QM/MM, with each timestep taking 3-6 orders of magnitude less wall time. ML/ML, therefore, incurs different software challenges: to minimize overhead, operations that take place each timestep must be closely integrated into the molecular dynamics (MD) software. Whilst some recent studies have explored similar ‘ML/MM’ ideas^40–43^, no generic implementation of ML/ML exists.

The central contribution of this work is ML-MIX, a package created to enable the whole workflow of accelerating a simulation within the popular open-source MD software LAMMPS^44,45^. Given an ‘expensive’ MLIP, we demonstrate how one can distill a ‘cheap’ linear MLIP that locally approximates the potential energy surface and then show how spatial mixing of the cheap and expensive MLIPs can be used to run accurate simulations at a fraction of the cost. The process is schematically depicted in Fig. 1. ML-MIX aims to be as generic as possible; it is able to perform ML/ML and ML/MM simulations with many native LAMMPS potentials out of the box, in parallel, on both CPU and GPU. A non-exhaustive list of compatible MLIPs includes SNAP^46^, ACE^13^, UF3^21^ and MACE^15^. An up-to-date full list of tested compatible potentials can be found in the ML-MIX GitHub repository^47^.Fig. 1A schematic overview of the simulation acceleration workflow.Left side: the constrained potential fitting process, whereby an accurate expensive potential (orange) is approximated in local regions of potential energy space by a simpler cheap potential (blue). Right side: schematic showing how simulation is accelerated; cheap potential is used to evaluate force components on atoms in `bulk-like' environments.

Here we present the results of six ML/ML case-studies, targeting one static and five dynamic quantities of interest in Si, Fe and W, demonstrating that ML-MIX can be used to generate results that agree with those from fully expensive simulations (and in the final case-study, experiment) at a fraction of the cost. The first case studies are designed primarily to demonstrate different features of the ML-MIX package. These particularly aim to highlight how ML/ML methods can produce accurate results for simple quantities, even for very small (~50 basis functions) linear, cheap potentials. We then present two more case studies, which were primarily conducted for scientific interest. Both studies concern the application of W as a diverter material in nuclear fusion reactors. In the first of these final two case-studies, we consider the thermally activated glide of screw dislocations in W at different stresses and temperatures. In the second, we consider the normally incident deposition of He at different energies into W at 1000 K. In this final case-study, we were able to reproduce (for the first time) experimental observations of the reflection coefficient of He into W for deposition energies up to 80 eV.

Results

Potential fitting

The starting point for our study is five ‘expensive’ potentials; four linear atomic cluster expansion (ACE) potentials (Si, Fe, W, W–He) fit using the ACEpotentials.jl package^48^, and one MACE potential created from fine-tuning the mpa-0-medium foundation model to W–He data using multihead replay fine-tuning^19^. Table 1 shows the fitting parameters of each expensive potential and where to find the corresponding DFT data and parameters. Additional hyperparameters for the Fe, W and W–He ACE potentials are specified in Section S1 of the Supplementary Information. The root-mean-squared per-atom-energy and force errors (RMSEs) for all potentials on their fitting datasets are shown in the upper left quadrant of Table 2.Table 1A summary of expensive potential fitting parameters and source of corresponding DFT dataExpensive potentialKey fitting parametersDFTSi ACEν = 4*, p* = 20, rcut = 6.0 Å^16^Fe ACEν = 3, p = 20, rcut = 5.5 Å^82^W–He ACEν = 3, p = 20, rcut = 5.0 Å^56,73^W ACEν = 3, p = 21, rcut = 5.0 ÅData:^81^ Params:^56,73^W–He MACEmpa-0-medium^19^ fine tuned, 128 L = 1 channels, 128 L = 0 channels, 2 layers, rcut = 6.0 ÅPretraining data: MPtraj^83^, sAlex^84,85^. Fine tuning: Data:^56,73,81^ Params:^56,73^For the ACE potentials, ν refers to correlation order and p refers to maximum total polynomial degree. For the MACE potential, L refers to the message rank. In all cases, rcut is the cutoff radius for the local environment.Table 2. Training RMSE table for both expensive and cheap potentials on the full DFT datasets used for the creation of the expensive potentials and the small synthetic MD datasets used for the creation of the cheap potentialsRMSEsRMSEsExpensive potentialsEnergy (meV)Force (eV/Å)Cheap potentialsEnergy (meV)Force (eV/Å)Full DFT dataSi ACE (4, 20)2.140.075Si ACE (2, 10)106.50.27Fe ACE (3, 20)2.320.083Fe ACE (2, 10)18.480.156W–He ACE (3, 20)2.150.070W UF365490.16W ACE (3, 21)1.900.066W ACE (3, 14)4.4890.116W–He MACE4.100.075W SNAP**Synthetic MD dataSi ACE (4, 20)––Si ACE (2, 10)0.210.045Fe ACE (3, 20)––Fe ACE (2, 10)0.810.090W–He ACE (3, 21)––W UF31.020.053For the W–He potential, only configurations containing W are used (as the cheap UF3 potential cannot describe He interactions). Note that the synthetic data was generated by the expensive potentials, hence the lack of RMSE values in the lower left quadrant. Energy errors are per-atom. *The full DFT data used for the fitting of the W SNAP is not public, so these RMSEs are not stated.

The following cheap potentials were generated: (i–ii) Si and Fe linear ACE potentials (fitted with ACEpotentials.jl) with correlation orders of 2 and maximum total polynomial degrees of 10 and cutoffs of 6.0 Å and 5.5 Å, respectively, (iii) a W ultra-fast (UF3^21^) potential with a two-body cutoff of 6.0 Å and a 3-body cutoff of 3.5 Å and (iv) a W linear ACE potential (fitted with ACEpotentials.jl) with a correlation order of 3 and a maximum polynomial degree of 14 and a cutoff of 5.0 Å. We also use (v) the pure W SNAP potential generated by Wood and Thompson^49^. Whilst all these potentials have significantly lower evaluation cost than their expensive counterparts, (iv) and (v) are more sophisticated than (i–iii). This is due to (i–iii) only needing to describe relatively simple bulk interactions, whilst (iv) and (v) need to also describe surfaces and longer-range strain fields. In each case, cutoff radii were selected for the cheap potential that were the same or smaller than the expensive potential.

The correlation order 3, maximum total polynomial degree 14 W ACE was fit to the same DFT dataset as its expensive counterpart (the correlation order 3, maximum polynomial degree 21 W ACE). RMSE values for this fit are shown in the upper right quadrant of Table 2. Whilst the W SNAP potential was also fit to much of the same underlying DFT data^49^, the full fitting dataset is unpublished and therefore the RMSEs are left out of Table 2.

For ~5000 atoms, the speedup for the 3, 14 W ACE over the 3, 21 W ACE was ~4.4×, whilst the speedup of the pure W SNAP potential over the fine-tuned mpa-0-medium model was ~29×. Note that the MACE/SNAP speedup was measured on GPU, whilst all other quoted speedups were measured on CPU.

To generate the small cheap potentials (Si and Fe ACE, W UF3), the expensive potentials were used to generate ‘synthetic’ training data focused around a suitable local region of the potential energy surface. In each case, the synthetic fitting data were divided into two types: ‘soft constraint’ data which consisted of bulk-crystal high-temperature MD trajectories, and ‘hard constraint’ data which comprised a set of zero-temperature homogeneous lattice deformations. Constrained linear fitting was used to ensure that the elastic constants and lattice parameters closely matched between the cheap and expensive potentials. The constrained fitting scheme is described in detail in the section “Constrained potential fitting”.

The training energy and force root mean squared error (RMSE) values for the synthetically trained potentials are displayed in the lower right quadrant of Table 2. For comparison, potentials with identical hyperparameters were also fit directly to the same DFT data as the expensive reference. The training RMSE for these fits is shown in the upper right quadrant of Table 2. Compared to the corresponding expensive potentials, the speedup for force evaluations with the Si ACE potential was ~30×, the W UF3 potential was ~75×, and the Fe ACE potential was ~8×. Measured speedups were the same for systems of both 8000 and 250,000 atoms.

Table 3 shows the elastic constants predicted by the expensive and cheap potentials. For the potentials fit to synthetic data, it also includes rows demonstrating the impact of applying the hard constraint on the elastic constants and lattice parameters. The elastic constants were computed using the matscipy package^50^. In the case of the Si and Fe ACE potentials, adding the hard constraint data led to near-perfect matching of the reference elastic constants. For the pure W UF3 potential, the C_11_ and C_44_ elastic constants got closer to the underlying reference, but C_12_ was worse, something that can be attributed to the simplicity of the UF3 functional form. In all cases, the equilibrium lattice parameter of the model became closer to the expensive reference under constrained fitting.Table 3. Elastic constants and lattice parameters according to each potential (expensive, cheap unconstrained, cheap constrained) for Si, Fe and W–HeC_11_ (GPa)C_12_ (GPa)C_44_ (GPa)a (Å)Si ACE (4, 20)Expensive150.8 ± 1.056 ± 270.5 ± 0.45.461Si ACE (2, 10)Cheap (unconstrained)139 ± 381.29 ± 0.1252.2 ± 1.95.661Si ACE (2, 10)Cheap (constrained)148 ± 255 ± 369.8 ± 0.35.461Fe ACE (3, 20)Expensive283 ± 12154 ± 6105 ± 62.834Fe ACE (2, 10)Cheap (unconstrained)426 ± 16364 ± 2239 ± 23.360Fe ACE (2, 10)Cheap (constrained)284 ± 11154 ± 7106 ± 32.834W–He ACE (3, 20)Expensive520 ± 14199 ± 5143 ± 53.181W UF3Cheap (unconstrained)604 ± 12209 ± 12180 ± 63.136W UF3Cheap (constrained)527 ± 11157 ± 4157 ± 33.180W ACE (3, 21)Expensive520 ± 14199 ± 5143 ± 53.180W ACE (3, 14)Cheap (unconstrained)520 ± 14196 ± 5144 ± 53.180W–He MACEExpensive540 ± 15215 ± 8147 ± 43.179W SNAPCheap (unconstrained)518 ± 14196 ± 5144 ± 53.180Elastic constants are fit to stress-strain data using Bayesian linear regression, with the error (σ) in each value found from the maximum likelihood estimate. For Si and Fe, the addition of the hard constraints to the cheap potential fitting produces elastic constants that match the expensive within 1σ. For the W–He potential, the constrained fitting improves C_11_ (matching within 1σ) and C_44_ (matching within 2σ), but leads to a significantly worse C_12_ than the unconstrained case. In all cases, the constraint dramatically improves the lattice parameter.

Energy barrier for Si vacancy migration

We first selected a simple zero-temperature (static) quantity of interest for a materials defect. We measured the energy barrier for vacancy migration in a 10 × 10 × 10 (8000 atom) block of Si for both a fully-expensive simulation and an ML/ML simulation where the expensive potential was localized around the vacancy. For the ML/ML NEB, the combined relaxation was performed using ML-MIX. Figure 2 shows the results from this investigation. The energy barrier obtained through the ML/ML simulation (0.0548 eV) agrees with the all-expensive reference energy barrier (0.0543 eV) to within 1 meV. The NEB relaxation ran ~5× faster, close to the theoretical ‘zero-overhead’ value of 5.90× expected from the fraction of expensive atoms in the ML/ML simulation.Fig. 2. The energy barrier associated with vacancy migration in Si, obtained through all-expensive (blue curve) and ML/ML (orange curve) nudged elastic band simulations.For the ML/ML simulation, the final energies were obtained through one-shot energy evaluations of the relaxed structures with the expensive potential. The all-expensive energy barrier (0.0543 eV) agrees with the ML/ML energy barrier (0.0548 eV) within 1 meV. The tangent lines arise from the projected forces.

Stretched bond in silicon

Turning to a dynamical example, we investigated the quantity of interest proposed as a test for the quality of different QM/MM techniques by Bernstein et al.^35^—the average force on a single stretched and rigidly fixed Si bond in a block of bulk Si at 300 K.

In a 10 × 10 × 10 (8000 atom) supercell of bulk Si, a bond was stretched an additional 0.1 Å along its length and fixed. The average force on this stretched bond recorded over 100 ps was then found, both for all-expensive reference and a number of ML/ML simulations, where the impact of changing the radius of the expensive potential region around the stretched bond atoms \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}$$\end{document} on the measured average force was investigated. The first simulation, at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}=0$$\end{document}  Å, corresponded to an all-cheap simulation. Further simulations were conducted for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}=4$$\end{document}  Å, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}=6$$\end{document}  Å, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}=8$$\end{document}  Å and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}=10$$\end{document}  Å. A core region size of 2 Å was not tested as a core radius smaller than the Si-Si bond length of 2.35 Å results in no expensive atoms, and thus an identical result to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}=0$$\end{document}  Å case.

The results for these simulations are presented in Fig. 3. It can be seen that as the expensive potential replaces the cheap potential in the core ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core} > 0$$\end{document} ), the forces jump to closely matching the NVE reference. Images of the simulation domain highlighting the atoms treated expensively around the fixed bond are shown in the insets in Fig. 3.Fig. 3. Average force on stretched bond in silicon over 100 ps at 300 K for different expensive potential radii \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}$$\end{document} (blue points).Error bars represent the standard error in bond force measured over this time, with samples taken every 15 fs. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}=0$$\end{document} , which corresponds to using the cheap potential everywhere, the average bond force measured does not match the all-expensive NVE reference (red dashed line). Once the expensive potential is introduced at and around the stretched bond, this difference vanishes, and agreement between the mixed simulation and the NVE reference is within statistical error.

The theoretical and measured speedups are shown in Table 4. It can be seen that in serial, the measured speedups match the theoretical, whilst in parallel, the values are lower due to imperfect load balancing and parallel overheads. As \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}$$\end{document} increases, the number of expensive atoms in the simulation increases, and the overall speedup decreases. Load balancing is further discussed in Section S2 of the Supplementary Information.Table 4. Speedup table for mixed simulations in SiSimulationTheoretical speedupSerial speedupParallel speedupSi: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}=4$$\end{document} 12.9312.898.89Si: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}=6$$\end{document} 10.7910.7310.10Si: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}=8$$\end{document} 8.558.558.89Si: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}=10$$\end{document} 6.726.726.88The measured serial speedup closely tracks the theoretically predicted value, whilst the parallel speedup is worse than expected for small rcore due to imperfect load balancing.

Diffusion coefficients

Next, we moved to a less artificial dynamical quantity of interest. Here, we consider the diffusion of a dumbbell self-interstitial in Fe and an interstitial He impurity in W, both of which have been extensively studied^51–54^. For each of these quantities, an atomistic domain of (~8000 atoms) was used. Whilst domains of this size are not necessary for accurate quantities of interest in these simple cells (we found one can attain converged diffusion coefficients in cells of ~500 atoms and smaller, depending on the defect), we selected this size primarily for demonstration purposes; it is around the minimum cell size beyond which we believe ML/ML methods are worth implementing.

For both target defects, reference simulations were carried out using the expensive ACE potentials. These were followed by ML/ML simulations where the expensive potentials were limited to a small region around the diffusing defect, which was tracked as it moved. Figure 4 presents the results from these investigations side-by-side. For the Fe system, measurements of the diffusion coefficient as a function of temperature were taken with only the cheap potential, only the expensive potential, and using the ML/ML approach. For W-He, only the ML/ML and expensive-only simulation results are presented, as the UF3 potential was fit only for pure W and could not model W-He interactions. For both quantities of interest, the ML/ML simulations agreed well with the all-expensive reference calculations at multiple temperatures, despite the inaccuracy of the all-cheap simulation in the Fe case. The theoretical and measured speedups of the mixed simulations over the reference for both simulation types are shown in Table 5.Fig. 4. Diffusion coefficient as a function of temperature for the diffusion of the iron interstitial dumbbell defect through bulk iron and for helium in bulk tungsten.a The diffusion coefficient of the iron dumbbell for three cases: expensive reference potential only (blue), cheap potential only (orange) and mixed simulation (green). The mixed simulation leads to results that are well matched with the reference. Each point is the average of five diffusion coefficients measured in independent runs of 1 ns each. b The diffusion coefficient for He in W using a helium-tungsten ACE (blue) and a mixed simulation with both a helium-tungsten ACE and a UF3 tungsten potential (orange). Again, agreement is good. Note that there is no 'cheap-only' reference in this case because the UF3 potential cannot be used to model tungsten-helium interactions. Each point is the average of 25 diffusion coefficients measured in independent runs of 60 ps each. The error bars on the points in both a and b represent the standard error in each value.Table 5. Speedup table for mixed simulations of diffusion coefficients in Fe and W–HeSimulationTheoretical speedupSerial speedupParallel speedupFe3.793.802.23W–He10.5510.467.39The measured serial speedup closely tracks the theoretically predicted value, whilst the parallel speedup is lower due to imperfect load balancing.

Velocity of screw dislocation in W

The thermally activated glide of screw dislocations is known to play an important role in the deformation of W, with molecular dynamics simulations being a key tool used to understand their mechanics^55^. Following Bacon et al.^55^, we have used a periodic array of dislocations with Burgers vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=\frac{1}{2}\langle 111\rangle$$\end{document} to study the thermally activated dislocation glide velocities at a range of applied shear stresses (between 600 and 2000 MPa) and temperatures (300, 600 and 900 K). In each case, all-expensive simulations were conducted using the 3, 21 W ACE potential^56^ at a limited range of stresses before ML/ML simulations were conducted using ML-MIX for a greater range, with the expensive region covering only the dislocation.

The simulation domain is presented in Fig. 5a. The ML-MIX expensive region is highlighted. Figure 5b presents the measured velocity-stress profiles obtained in both the all-expensive and mixed simulations. Each data point is the mean of an ensemble of three repeats, with the error bars representing the standard error. The mixed and all-expensive curves agree with each other within error. The all-expensive simulations in this case were approximately 4× the cost of the mixed simulations.Fig. 5. Simulation set-up and results for thermally activated glide of screw dislocations in W at different shear stresses and temperatures.a Simulation set-up, where dimensions in x, y and z are 223Å × 221Å × 68Å(25b), respectively, for a total of 195,000 atoms. Stress was applied to the free surfaces at the top and bottom of the cell in 〈111〉 type directions. The expensive atoms around the dislocation (as detected using common neighbor analysis in LAMMPS) are shown in a column, where the colors correspond to the core (red), blending (green) and buffer (white) regions. b The velocity results at a range of temperatures and shear stresses. Green, orange and blue lines represent 300, 600 and 900 K, respectively. Results from all-expensive simulations are shown with dashed lines, whilst those from mixed simulations are shown with solid lines. For each point, three separate simulations were conducted and the results were averaged. Results from the mixed simulations agree well with those obtained from the all-expensive reference. The error bars on the points represent the standard error in each value.

He implantation into W

Tungsten is a well-known candidate for plasma-facing materials within modern fusion reactors. Proposed as both a diverter and first-wall material^57,58^, one of its key properties is the ability to withstand high-energy flux from ions and neutrons. Part of this flux includes impaction from He, the fusion product. In order to predict the damage level and subsequently the longevity of the W material, it is crucial to know (1) how much of the He is expected to enter the W for different incident energies and (2) how deep the He penetrates. In a previous molecular dynamics study, it was found that these quantities are highly dependent on the choice of interatomic potential^59^, and so far no reconciliation has been made between molecular dynamics results and the limited experimental data that exists^59,60^.

To compare with experiment, we specifically targeted simulations of He implantation into {100} W surfaces at 1000 K. Initially, we performed all-expensive simulations using a fine-tuned version of the MACE-mpa-0 foundation model^19^. We then conducted ML/ML simulations in which the expensive potential interactions were limited to the immediate vicinity of the He atom. In both all-expensive and ML/ML simulations, all computation was conducted on GPU. Figure 6 shows the measured reflection coefficient as a function of incident He energy. The all-expensive (mpa-0-medium-finetuned) data (blue points), ML-MIX data (orange points) and experimental data (green points) all broadly agree within error. Average implantation depths were also measured for each energy. These are shown in Fig. 7, where it can be seen that results from the ML/ML simulation are in broadly good agreement with the all-expensive. The reduced cost of the ML/ML simulations allowed extension to larger system sizes and thus greater He energies and implantation depths.Fig. 6. Reflection coefficient as a function of incident He energy in eV for normal deposition onto a 1000 K W {100} W surface.Experimental results from Van Gorkum et al. are shown in green^60^. Results from the all-expensive potential (fine-tuned mpa-0-medium model) are displayed in blue, and results obtained with ML-MIX are shown in orange. Results from a previous MD simulation using a W EAM and W-He pair potential are displayed in purple^59^. For the all expensive, to limit computational cost, only energies up to 40 eV were tested with 100 repeats each. With ML-MIX, energies up to 100 eV were repeated 500 times each. Error bars are obtained from the 95% confidence interval of a fitted Beta distribution.Fig. 7. Implantation depth as a function of incident He energy in eV for normal deposition onto a 1000 K W {100} W surface.Results from the all-expensive potential (fine-tuned mpa-0-medium model) are displayed in blue, and results obtained with ML-MIX are shown in orange. For all expensive, to limit computational cost, only energies up to 40 eV were tested with a 100 repeats each. With ML-MIX, energies up to 100 eV were repeated 500 times each. Error bars represent the standard error in the mean values.

The speedup of the mixed simulation over the expensive simulation was approximately 10× for simulations where a reflection takes place and approximately 5× for simulations where the He implanted.

Discussion

We have presented ML-MIX, a hybrid simulation package that aims to break existing cost-accuracy tradeoffs in simulations by localizing accuracy requirements, enabling simulation speed-ups of 4−10 × with minimal loss in accuracy on quantities of interest. ML-MIX has the following advantages that set it apart from previously published similar ML/MM schemes: (1) ML-MIX is an easy to install, direct plugin to LAMMPS, which is fast and has no external dependencies, and aims to be compatible with a large number of existing pair_styles. (2) ML-MIX functions as an external wrapper around LAMMPS pair_styles, meaning there is no necessary editing of C++ code, (3) ML-MIX works with existing LAMMPS parallelism and load-balancing, and (4) ML-MIX (though our Kokkos implementation) is compatible with pair_styles that execute on both CPU and GPU.

We have described two ways in which cheap potentials can be generated; either by fitting directly to the same underlying DFT data as the expensive reference, or by fitting to synthetic data generated by the expensive potential in a targeted region of the potential energy surface. We found that fitting to synthetic data was particularly important for very lightweight potentials, as for these we observed that attempting to fit to the diverse reference DFT data directly lead to unacceptably high RMSEs. To ensure close matching of specific properties, we developed a constrained fitting scheme, the details of which are described in section “Constrained potential fitting”. For linear ACE potentials, we found that imposing constraints allowed for exact matching of elastic constants between the cheap and expensive potentials. For UF3 potentials, a less flexible architecture, we have shown that under hard constraint elastic constants improve up to some limit, beyond which any further improvements were limited by the misspecification of the functional form. We stress that whilst in this study we have specifically discussed constrained fitting in the context of distilling one MLIP from another, it would be easy to apply this methodology to constrain the fitting of a large, flexible MLIP to DFT data directly, extending recent retraining approaches^38^.

We found that constrained fitting is particularly important for small potentials of <50 basis functions; without it, errors are far more pronounced at small strains. This can lead to erroneous Eshelby-inclusion^61^ like strain fields in the vicinity of the ML/ML potential boundary under small deformations. This is explored in more detail in Section S3 of the Supplementary Information.

The first case-study we present in this paper was a nudged elastic band energy barrier calculation for the migration of an Si vacancy. We showed that, using an ML/ML simulation, a value for the migration barrier can be attained which is within 1 meV of the all-expensive reference at a fifth of the computational cost. This speedup, whilst already considerable, is not as great as what is seen for ML/ML MD simulations on the same 8000 atom domain—this is because a larger fixed set of expensive atoms is necessary for a NEB to converge. It is important to point out that total energies are undefined in ML-MIX ML/ML simulations; the final energies were computed by performing one-shot energy evaluations of the relaxed structures using only the expensive potential.

Next, we presented an investigation into the average force on a stretched and rigidly fixed bond in an Si bulk structure. This was selected to mirror previous tests performed in QM/MM studies^35^. It should be noted that for an \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}$$\end{document} of 0 Å (cheap potential only) the average force measured on the stretched bond is within 5% of the correct value, which is a significantly smaller error than one typically sees in QM/MM simulations. Adding even a small amount of expensive potential region around the stretched bond led to an immediate improvement in the measured average bond force, leading to a value which was within 1 standard error of the NVE reference value. This behavior persisted for larger \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}$$\end{document} values. The measured serial speed-up closely matched theoretical predictions, showing that our implementation incurs minimal overhead due to the mixing process. However, in parallel testing (on 48 cores) the measured speedup is below theoretical predictions, in particular for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{core}=4$$\end{document}  Å case. We attribute this performance to the limited strong scalability of the expensive interatomic potential at low atom count, leading to significant load balancing issues in LAMMPS when there are few expensive atoms. The limited speed-up is thus a property of the expensive interatomic potential rather than our load balancing protocol; near-theoretical parallel speed-ups are achievable on larger systems with simple load balancing requirements. This is demonstrated in Section S2 of the Supplementary Information.

We also presented case studies on defect diffusion coefficients in Fe and W. The measured diffusion coefficients consistently agreed with the all-expensive reference simulations within statistical error. The close matching of the measured diffusion coefficients with the reference values suggest that one only requires spatially local accuracy to attain correct defect diffusion dynamics in materials; corroborating findings from earlier QM/MM studies^62^.

Next, we considered some more complex and scientifically relevant case-studies, aiming to demonstrate the type of simulations where we believe ML-MIX will find application. These were (1) a study of the thermal glide of screw dislocations in W at a range of temperatures and pressures, and (2) a study of normal deposition of He in W at 1000 K.

We considered the thermally activated glide of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=\frac{1}{2}\langle 111\rangle$$\end{document} screw dislocations in W. We found that in a mixed ACE/ACE simulation, we were able to reproduce the dislocation velocities measured in the all-expensive reference, meaning that the mixed simulation was able to reproduce the correct dynamics of the kink-nucleation-propagation mechanism that controls the glide of screw dislocations in W. Whilst there was scope to study longer dislocation lines, we found increasing the length of the dislocation increased the probability of cross-kinking—a phenomenon that can lead to debris formation. This can cause nonphysical interference with subsequent dislocation motion when the dislocation re-encounters the debris due to periodicity. As such, we leave an exploration of larger dislocation studies with ML-MIX to future work.

Finally, we simulated the normal implantation of He into (100) W surfaces at He energies up to 70 eV at 1000 K. This led to a novel result; for the first time, experimental observations of reflection coefficient were reproduced up to 40 eV with the mpa-0-medium foundation model fine-tuned to W–He DFT data. Following this, we carried out the same simulations with ML-MIX, where we limited interactions with the expensive potential to directly around the He and used a pure W SNAP potential for the remainder of the interactions. In order to ensure efficient evaluation of the expensive MACE potential, ML-MIX was ported to GPU using Kokkos^63,64^. Using ML-MIX, we showed that (1) up to 40 eV results agreed well with the all expensive simulation and experiment, and (2) above 40 eV results continue to line up closely with experiment, with deviations only occuring at relatively high energy (>80 eV). Measurements of implantation depth were also taken in each case, and those from the mixed simulations aligned well with those from the all-expensive simulations. Due to the relative low cost of the mixed simulations, we were able to perform more repeats and create the larger cells necessary for high energy simulations without needing to parallelize over more GPUs.

It is important to note that the experiment we have drawn comparisons with involved incident He^+^ ions^60^, whilst in this study we are considering incident He atoms. One possible explanation for our close agreement despite this is that in reality, as the He^+^ approaches the surface, it gains an electron from the W conduction band on timescales faster than the movement of the nuclei. A recent time-dependent DFT study by Ito et al.^65^ studying plasma-wall interactions observed electron transfer in the case of He^2+^, but was inconclusive on He^+^. We believe that reconciling this and performing a wider exploration of the implantation parameter space are both beyond the scope of this primarily method-focused paper, and therefore leave them as interesting directions for future work.

A key limitation in ML/ML (and indeed, QM/MM) simulations arises from the lack of energy conservation; if a force mixing method is used to combine potentials (as it is here), simulations do not conserve energy. This means that a thermostat may be required to keep simulations from substantially heating up as they progress. The energy drift in well-matched ML/ML simulations is relatively small ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{O}}(10\,K/ps)$$\end{document} ), meaning thermostats can be weak enough to not impact dynamics—but this should always be checked. Force mixing is discussed further in section “Spatial potential mixing”, and a full investigation into the factors affecting energy conservation in ML/ML simulations is presented in Section S4 of the Supplementary Information.

A possible solution to this problem is the use of so-called ‘energy-mixing’^35^, where a Hamiltonian is constructed for the mixed system, for example by writing the total energy as a sum of local site energies from each potential. Upon testing this scheme however, we found the forces derived had errors of >50% in the vicinity of the ML/ML boundary, even when using potentials that nominally are well matched in force-mixing. This backs up previous work on energy-mixing for QM/MM conducted by Bernstein et al.^35^ and expanded on by Chen and Ortner^33^ where it was found that eliminating these spurious ‘ghost forces’ in energy-mixed QM/MM for infinitesimal displacements away from equilibrium required the MM potential have an exactly matching force-constant matrix (second order derivatives of the site energies with respect to atomic positions). This requirement is stringent and as such would require more flexible (and thus expensive) cheap potentials, and even if satisfied only would eliminate ghost forces for atoms near equilibrium. Additionally, if the regions change over time as atoms move (a common feature of long-timescale ML/ML simulations), energy would not be conserved anyway; as regions update, potential energy would be continually added or removed from atoms. Conversely, a buffered force-mixing approach naturally has no spurious ghost-forces and can be made to work with far simpler cheap potentials. In view of this, we believe (despite the lack of energy conservation) that force-mixing is the clear choice for ML/ML simulations.

A further limitation is that we have made no attempt to quantify or estimate how the error associated with the size of the expensive region selected may propagate through to quantities of interest. Uncertainty quantification (UQ) for MLIPs and QM/MM is a very active area of research^66–69^, and we believe that applications of UQ to ML/ML simulations is a very exciting possible direction for future work.

Through these studies, we have demonstrated the flexibility and efficiency of the ML-MIX package. In our initial, simple case-studies, we show that ML-MIX can be used to accelerate simulations without compromising accuracy both in simple cases where the same atoms are always involved in the defects (i.e., for He in W) and for more complex cases where different atoms are involved as the defects move (i.e., for Fe dumbbell interstitials in Fe). In two further case-studies, we have shown that ML-MIX can be used to attain a significant speedup without any meaningful loss in accuracy in applications of genuine scientific interest. In the first of these, we present results for the speed of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=\frac{1}{2}\langle 111\rangle$$\end{document} screw dislocations in W at a range of temperatures and shear stresses. In the second, we provide estimates of reflection coefficients and implantation depths for normally deposited He into 1000 K {100} W surfaces. Our measured reflection coefficients match experimental observations up to a deposition energy of 80 eV. Looking forwards, we anticipate that ML-MIX will find application in a broad range of atomic simulations of heterogeneous systems, including the calculation of defect free energy barriers^70^, interface migration^71^ and nucleation phenomena^72^.

Methods

Constrained potential fitting

The constrained potential fitting process is shown schematically in Fig. 8. By fitting a cheap potential to highly localized regions of the potential energy surface of the expensive potential, it is possible to attain high accuracy in local regions whilst maintaining a low complexity (and thus low evaluation cost).Fig. 8A schematic representation of the constrained fitting of a cheap potential to an expensive one.a A representation of the local potential energy surface of the expensive reference model (orange curve). Data has been sampled from this potential energy surface (crosses). b The cheap model, fit to this data in order to locally match the expensive model. It approximately matches the soft' constraint data (pink crosses) and exactly matches the hard' constraint data (green crosses).

We first consider the challenge of selecting suitable constraints. The requirements of the cheap potential are split into two types: soft constraints (loose matching), and hard constraints (tight matching). Archetypal hard constraints are the elastic constants, to ensure seamless matching of the long range elastic stress fields between the cheap and expensive potentials in a solid state simulation^35^.

Writing the design matrices for the hard and soft constraints as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{A}}}_{H}\in {{\mathbb{R}}}^{{N}_{H}\times {N}_{D}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{A}}}_{S}\in {{\mathbb{R}}}^{{N}_{S}\times {N}_{D}}$$\end{document} , and the fitting data as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{y}}}_{H}\in {{\mathbb{R}}}^{{N}_{H}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{y}}}_{S}\in {{\mathbb{R}}}^{{N}_{S}}$$\end{document} , we can express the potential parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{c}}\in {{\mathbb{R}}}^{{N}_{D}}$$\end{document} as a solution to the constrained optimization problem:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop{\min }\limits_{{\bf{c}},| | {{\bf{A}}}_{{\mathrm{H}}}{\bf{c}}-{{\bf{y}}}_{\mathrm{H}}| {| }^{2} < \alpha }(| | {{\bf{y}}}_{{\mathrm{S}}}-{{\bf{A}}}_{{\mathrm{S}}}{\bf{c}}| {| }^{2})$$\end{document}

where regularization terms have been excluded for simplicity. The hard constraint is imposed by constraining the fit to the subspace ∣∣AHc − yH∣∣^2^ < α, where α is the maximum allowed error for the model on the hard constraint data. Expressing this constraint through a Lagrange multiplier argument gives us the augmented Lagrangian:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{L}}({\bf{c}},\lambda )=| | {{\bf{y}}}_{{\mathrm{S}}}-{{\bf{A}}}_{\mathrm{S}}{\bf{c}}| {| }^{2}+\lambda (| | {{\bf{y}}}_{{\mathrm{H}}}-{{\bf{A}}}_{{\mathrm{H}}}{\bf{c}}| {| }^{2}-\alpha ).$$\end{document}

By requiring that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\nabla }_{{\bf{c}}}{\mathcal{L}}({\bf{c}},\lambda )=0$$\end{document} , we can express optimal constrained model parameters c at fixed λ as a solution to the unconstrained problem:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop{\min }\limits_{{\bf{c}}}(| | {{\bf{y}}}_{{\mathrm{S}}}-{{\bf{A}}}_{{\mathrm{S}}}{\bf{c}}| {| }^{2}+\lambda | | {{\bf{y}}}_{{\mathrm{H}}}-{{\bf{A}}}_{\mathrm{H}}{\bf{c}}| {| }^{2}).$$\end{document}

From this expression, it is clear that the quadratic constraint is equivalent to adding the hard constraint configurations to the fit with weights scaled by λ. Obtaining a solution to the full constrained problem therefore requires finding the minimum λ such that the solution to (3) lies on the surface of the constraint subspace ∣∣AHc − yH∣∣^2^ = α.

Scripts to perform constrained linear fitting for ACE and UF3 potentials are provided in the accompanying code^47^. To perform the λ hyperparameter search, the linear problem is solved repeatedly. Initially, λ is increased logarithmically until ∣∣AHc − yH∣∣^2^ < α. Interval bisection is then used to find λ such that ∣∣AHc − yH∣∣^2^ − α < tol, where tol is set by default as α/10. An example λ search to satisfy a constraint of α = 10^−7^ is shown in Fig. 9. AH and AS only need to be assembled once, and as the cheap potentials are small (often ≲ 100 parameters), the whole constrained fitting process takes approximately 5 min on a single CPU core.Fig. 9. Automated selection of lambda value during a single constrained fit to satisfy a given alpha.The vertical axis shows the energy RMSE per atom. The target alpha value is shown as a dashed red horizontal line. Each blue point represents a step on the search.

Generation of synthetic data

Simple high-temperature bulk equilibrium dynamics was selected as the soft constraint data. For Si, a 4096-atom diamond structure bulk cell was simulated at 500 K for 10 ps, with snapshots saved every 0.5 ps. For W and Fe, 2000-atom body centered cubic (BCC) bulk cells were simulated at 1200 K, with snapshots saved every 0.1 ps. The soft constraints consisted of the total energy and forces from each snapshot.

For the hard constraints, following^38^, data were generated by applying multiple shear and expansion deformations to 5 × 5 × 5 bulk supercell of Si, Fe and W lattices and saving the resulting energies. Seventeen configurations were generated between ±0.5% strain for all of the following strain states: hydrostatic, uniaxial [100], uniaxial [111], shear ([100], (010)), shear ([110], (001)) and shear ([110], (1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{1}$$\end{document} 2)). As the energy of weak, homogeneously deformed lattice structures are described fully by the elastic constants^38^, fitting to the selected hard constraint enforced elastic constant matching between cheap and expensive potentials.

For the Si and Fe ACE potentials, an α value was selected that corresponded to an enforced energy error of less than 10^−7^ eV/atom on the hard constraint data. For Si, the corresponding λ value was 3.25 × 10^9^. For Fe, λ was found to be 5.5 × 10^9^. To a reader familiar with MLIPs, a target per-atom-energy error of α = 10^−7^ eV/atom may seem very small (a standard aim when fitting MLIPs is a per-atom-energy error of around 1 meV). However, we found that pushing α to such small values was necessary to achieve the tight property matching that we were aiming for. This is demonstrated for the Si cheap potential in Fig. 10; such small α values were necessary to match elastic constants within a fitting error of ~ 2 GPa and lattice parameter within a fitting error of 10^−5^ Å. For the W potential the constrained fitting process failed due to the UF3 potential not being capable of simultaneously fitting the hard constraint data and maintain a low RMSE on the soft constraints. By searching the space of λ values, it was found that above approximately λ = 4.3 × 10^5^, ∣∣AHc − yH∣∣^2^ stopped decreasing. λ was therefore set to 4.3 × 10^5^.Fig. 10. Demonstration of convergence of key quantities with hard constraint strength for the Si ACE 2, 10 cheap model.Error is shown compared to the Si ACE 4, 20 expensive model. In each case, the unconstrained case is shown by a dashed horizontal red line and the result at each constraint strength (α) is shown by a blue point. a The L2 norm of the error over C_11_, C_12_ and C_44_. b The absolute error in lattice parameter. Note that at a constraint strength of 10^−7^, the lattice parameter error was 0 up to the number of decimal places measured (5).

Additional DFT data

In the case of the cheap and expensive potential fit for the W dislocation and He–W implantation studies, additional DFT data was acquired. The DFT parameters used matched those described in ref. ^73^. For the dislocation mobility simulations, the dataset in ref. ^73^ was extended to include 1000 K surface and dislocation kink configurations. Forty structures for both (110) and (112) surfaces were generated, each containing 45 and 48 atoms, respectively. Twelve structures for each of the two kink types (left and right) were generated, containing 330 and 345 atoms, respectively. For the He–W implantation study, 129 snapshots of small 91-atom implantation cells were generated, from simulations of He implantations at 10 eV and 100 eV. These were added to the full He–W dataset provided by Nutter et al.^56^, which was then used for generation of the expensive and cheap potentials for the He–W implantation simulation.

Region tracking

When running an ML/ML simulation, it is crucial that the process of tracking defects and building the expensive potential region(s) is: (1) Accurate—if the cheap potential is erroneously used in a region that requires the expensive potential, it could impact on the accuracy of the simulation. (2) Flexible—any defect type should be trackable. (3) Fast—region building should not add considerable overhead to a simulation.

In order to satisfy all three of these simultaneously, the region building and tracking were implemented directly in LAMMPS through a custom fix. By re-using the same neighbor lists that are built for the pair style evaluation step, this fix can quickly rebuild regions. Additionally, the fix is MPI local; whilst a limited amount of communication between neighboring processes is necessary, no expensive global communication is required.

The fix builds regions in two distinct stages: (1) Identify the ‘seed’ atoms through user-defined criteria. (2) Construct the regions around these atoms.

Seed atoms are identified in one of two ways; through a LAMMPS group or through querying the output vector of a separately defined LAMMPS fix. Tracking through a predefined LAMMPS group is for situations where seed atoms are unchanging, and are known at the start of the simulation, as is the case with He atoms in W. Tracking through the output vector of a user-supplied fix is useful for tracking defects which are not tied to any one particular set of atoms, for example a vacancy, dislocation or an Fe dumbbell interstitial. Seed atom identification is schematically represented in the first two panels of Fig. 11.Fig. 11. Region construction schematic.From left to right: starting with a cell containing an initial defect, it is identified (orange atoms in panel two) either using a user-specified LAMMPS group or by the output vector from a user-defined LAMMPS fix. The defect atoms are then iterated over, with atoms that sit less than one core radius out from any defect atom being added to the core region. This construction process is then repeated for the atoms in the blended region.

Once the seed atoms are identified, regions are constructed around these seed atoms, using a method that is schematically identified in the last two panels of Fig. 11.

During region construction two arrays are populated, i2_potential and d2_eval. Both arrays are of size number of potentials × number of atoms. The i2_potential array holds integers and is used to determine which atoms are to be evaluated with each potential; i2_potential[1][10] = 1 indicates that atom 10 should be evaluated with potential 1, whilst i2_potential[1][10] = 0 indicates that it should not be. d2_eval holds double-precision real numbers and describes how much each potential contributes to forces on atoms. By allowing d2_eval to hold fractional values, it is possible to introduce a blending region in which atoms move under a mixture of the forces from the two separate potentials; d2_eval[0][10] = 0.4 and d2_eval[1][10] = 0.6 indicates that atom 10 should get 40% of its’ forces from potential 0 and 60% from potential 1. Two possible blending functions are implemented in ML-MIX—linear blending and cubic blending. A brief discussion of the impact on energy conservation of these two blending schemes is presented in Section S4 of the Supplementary Information.

Three regions are built in the region construction process, the core region, in which the force on all the atoms is given entirely by the expensive potential, the blending region, in which the forces are a mixture between the cheap and expensive potentials and the buffer region, which is crucial for the accuracy force-mixing scheme, and is discussed in the next section. The size of each region is specified by user controlled parameters: r_core, r_blend and r_buffer which are passed to the fix.

Algorithms 1, 2 and 3 describe how these arrays are populated in the region construction process for one of the two potentials (potential 0). Note that in algorithm 2, a linear blending function to mix together potentials in the blending region is used. This is not mandatory and in principle could be replaced by any blending function. So long as there are only two potentials, it is trivial to populate i2_potential and d2_eval for the other potential once the regions have been constructed around the seed atoms. Currently, ML-MIX is limited to two potentials only.

Algorithm 1

Build core region

1: Initialize i2_potential to 0

2: Initialize d2_eval to 0.0

3: for i ← seed atoms do

4:  i2_potential[i][0] = 1

5:   d2_eval[i][0] = 1.0

6:   for j ← neighbors of atom i do

7:    if ∣∣pos(i) − pos(j)∣∣ < r_core then

8:      i2_potential[j][0] = 1

9:     d2_eval[j][0] = 1.0

10:   end if

11:  end for

12: end for

13: Communicate()

Algorithm 2

Build blended region

1: for i ← core atoms do

2:   for j ← neighbors of atom i do

3:     if ∣∣pos(i) − pos(j)∣∣ < r_blend then

4:      i2_potential[j][0] = 1

5:      r = ∣∣pos(i) − pos(j)∣∣

6:      v = max(d2_eval[j][0], 1.0-(r/r_blend))

7:      d2_eval[j][0] = v

8:     end if

9:   end for

10: end for

11: Communicate()

Algorithm 3

Build buffer region

1: for i ← blended atoms do

2:   for j ← neighbors of atom i do

3:     if ∣∣pos(i) − pos(j)∣∣ < r_buffer then

4:       i2_potential[j][0] = 1

5:     end if

6:   end for

7: end for

8: Communicate()

Spatial potential mixing

The mixing scheme implemented in ML-MIX is force-mixing, which is a scheme commonly used in QM/MM. For in-depth reviews of force-mixing in the QM/MM context, please see Bernstein et al.^35^. Here we briefly review the force-mixing method, and discuss simplifications that can be made when using two local potentials.

The process of force mixing is schematically shown in Fig. 12. To obtain the forces, the domain is split into separate segments where an additional buffer region is included to stop individual regions seeing an artificial surface. The forces on these regions are then evaluated to create separate parts of the overall force vector. The full force vector is constructed as a combination of these separate vectors. In QM/MM, as DFT is a non-local method, the size of the buffer region around the QM domain is a crucial convergence parameter—one needs to balance the error due to the separation of the region of interest from the artificial surface against additional cost incurred as atoms are added.Fig. 12. An illustration of force-mixing with interatomic potentials.A mixed domain is split into two non-interacting subdomains (blue and orange), which are composed of all the atoms in the core and blending regions for each potential (solid line), plus an extra buffer region containing atoms up to the potential cutoff radius out (dotted line). Partial force vectors Fc and Fe are obtained by evaluating the cheap and expensive subdomains respectively. These partial vectors are then combined into the full force vector on the system F. In this process, forces on buffer atoms are discarded and forces on blending atoms are mixed between potentials.

Local potentials simplify this. Consider a set of atoms Λk, which is a subset within a larger domain Λk ⊆ Λ. We want to attain correct forces on atoms in Λk with a local potential Φk that has a defined cutoff radius rcutoff. The force on atom i ∈ Λk can be written in terms of local energies as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${F}_{i}=\frac{\partial E}{\partial {{\bf{x}}}_{i}}=\mathop{\sum }\limits_{j}^{j\in {\mathrm{neigh}}(i)}\frac{\partial {E}_{j}^{k}}{\partial {{\bf{x}}}_{i}}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{j}^{k}$$\end{document} represents the local energy of atom j as evaluated by Φk. Due to the local nature of the potential, we only need to evaluate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{j}^{k}$$\end{document} for atoms that lie within rcutoff of atom i. For an atom that lies on the edge of Λk, it is clear that attaining the correct force requires evaluation of atoms up to rcutoff outside Λk. Defining the distance between two atoms i, j as d(i, j), the buffer region Λk-buffer is given by:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{l}{\Lambda }_{k-{\rm{b}}{\rm{u}}\mathrm{ff}{\rm{e}}{\rm{r}}}=\{i\in \Lambda \backslash {\Lambda }_{k}\,| \,\exists j\in {\Lambda }_{k}\\ \,{\mathrm{and}}\,d(i,j)\le {r}_{{\mathrm{cutoff}}}\}.\end{array}$$\end{document}

Note that through the dependence of local energies on local atomic positions, there is an implicit dependence on atomic positions up to two cutoff radii out. Hence, r_buffer need not be any greater than the potential receptive field, as at this radius the forces match what they would be in non-mixed simulations on each domain. For MLIP architectures with a large receptive field due to message passing (e.g., MACE), one often does not need a buffer that captures it in its entirety; force errors reduce in size rapidly with buffer radius.

Force mixing is implemented in ML-MIX through a wrapper pair_style which evaluates multiple sub-pair_styles in turn. Algorithm 4 demonstrates this process.

Algorithm 4

Force mixing

1: Initialize forces to 0

2: for i ← potentials do

3:   for j ← atoms do

4:    if i2_potential[i][j] == 0 then

5:       Remove j from neighbor list

6:    end if

7:   end for

8:   temp_forces = potentials(i) → compute()

9:   temp_forces = temp_forces * d2_eval[i][:]

10:   forces = forces + temp_forces

11:   Restore neighbor list

12: end for

Si vacancy migration

To generate the vacancy structure, a 10 × 10 × 10 (8000) atom Si supercell was generated, and the atom nearest the center was removed. Initial guesses for the minimum energy pathway were formed by linearly interpolating the migration of a neighboring atom into this vacancy. Including the start and end states, 9 images were generated along this pathway, which were then relaxed using the nudged-elastic-band (NEB) method^74^ as implemented in LAMMPS. Further analysis was conducted in Python using ASE^75^. For the ML/ML simulation, the seed atoms for the expensive potential region were a union of all atoms immediately surrounding the vacancy in the final and end states. These neighbors were identified using coordination analysis with a cutoff of 4.0 Å, with all atoms that have 15 neighbors being selected. An additional 4.0 Å of core region was added around these atoms, as well as 4.0 Å of blending region. 6.0 Å of buffer region was then included, which matches the cutoff of the cheap and expensive potentials. Region rebuilding was switched off; provided expensive atoms do not move far this allows relaxations to progress smoother and faster by stopping atoms changing the proportion of forces they gain from each potential as the minimization progresses. For both the all-expensive and ML/ML NEB simulations, no energy tolerance was set (as total energies are not present in the ML/ML simulation) and the force tolerance was set such that ∥F∥∞ < 0.0005 eV/Å. Tighter tolerances than this would frequently fail to converge for the ML/ML simulation. If tighter tolerances are necessary, then we recommend switching to the all-expensive regime for the last steps of a relaxation. The final energies were obtained through one-shot energy evaluations of the relaxed structures with the expensive potential.

Si stretched bond

To generate the structure, a 10 × 10 × 10 (8000 atom) periodic supercell of bulk Si was constructed and a single bond was stretched an additional 0.1 Å along its length and fixed using the rigid package in LAMMPS^76^. This block was then pre-thermalized to 300 K in a 2 ps simulation with a Langevin thermostat using only the expensive potential. This thermalized initial structure was used as the starting point for all subsequent comparison simulations, in which the average force on this bond was recorded over 100 ps. For the all-expensive simulations, an NVE ensemble was used. For the ML/ML simulations, the seed atoms were selected to be the atoms involved in the stretched bond, no blending region was used (abrupt force mixing), as well as a weak thermostat of time constant 2.0 ps applied to only the cheap potential region in order to stop temperature drift due to the ML/ML boundary. For all simulations a buffer region size of 6.0 Å was used (which matches the cutoff of both potentials). The size of the core region was varied between simulations, as described in section “Stretched bond in silicon”. The expensive potential region was rebuilt around this seed atom every timestep. To compute the average bond force, samples of atomic forces were taken every 15 fs, long enough to avoid significant correlation between samples.

Diffusion coefficients

To calculate the diffusion coefficients, individual defects were placed into periodic cells of size 16 × 16 × 16 unit cells (~8000 atoms). Cells were heated to the target simulation temperature over 2 ps with a Langevin thermostat. For the Fe simulation, diffusion coefficients were measured at 800 K, 900 K, 1000 K and 1100 K. For W–He, diffusion coefficients were measured at 400 K, 600 K and 800 K. The reference simulations were then switched to NVE, whilst the ML/ML simulations had a weak Langevin thermostat (damping parameter - 2.0 ps) applied only to the cheap potential region. For the Fe simulation, each data point consisted of the average of five diffusion coefficients derived from the mean square displacements from non-overlapping time-lags of 5 × 1 ns simulations. For the W–He simulations, each point was the average of twenty-five diffusion coefficients derived from the mean square displacements from non-overlapping time-lags from 25 × 60 ps simulations. In all cases, the error was estimated from the standard error in the averaged diffusion coefficients. For the W–He ML/ML simulation, the seed atom was selected to be the single He present in the simulation, and regions were rebuilt around this seed atom every timestep. For the Fe ML/ML simulation, seed atoms were selected by finding the most highly coordinated atoms (on average) in the simulation every 100 timesteps. For this, measurements of atomic coordination Z (with a cutoff of 2.0 Å) were taken every 10 timesteps, and these were averaged using fix ave/atom. New seed atoms were selected according to the criteria that 〈Z〉 > 0.4. Regions were rebuilt around these seed atoms every timestep. For the first 100 timesteps in the Fe ML/ML simulation, only the expensive potential was used. Both simulations used a core potential region of 6.0 Å around the defect atoms, as well as a blending region of 4.0 Å. For the Fe potential, a buffer width of 5.5 Å was used, whilst for the W–He potential, a width of 6.0 Å was used, matching the UF3 potential cutoff.

W dislocation glide

Simulating the thermally activated glide of W screw dislocations is known to require ≈ 10^5^ atoms due to finite-size effects^55,77,78^. These are (1) in x, the cell length must be large enough to prevent self-interaction and to ensure that the work done can dissipate without substantial heating, (2) in y, the surface must be far enough away to prevent attraction and dislocation annihilation, and (3) in z, to ensure that the dislocation moves by the correct mechanism of the nucleation and propagation of kink pairs. In this work, the dimensions of the simulation cell in the glide direction (along \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x=[10\overline{1}]$$\end{document} ), glide plane normal (along \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y=[\overline{1}2\overline{1}]$$\end{document} ), and dislocation line (along z = [111]) were chosen to be 223 Å × 221 Å × 68 Å (25b), respectively, for a total of 195,000 atoms. At each temperature simulated, the thermally expanded lattice parameter was measured and the cell was rescaled to match. In each simulation, the dislocation was placed at the center of the simulation cell, where the atomic displacement field was generated using the matscipy dislocation module^50^. Periodicity was imposed in x and z, with free surfaces in y. Shear stresses σyz were applied to these free surfaces, simulating {112} loading. During glide, the dislocation moved on alternating {110} planes.

Using a 2 fs timestep, the cell was equilibrated to the chosen temperature for 50 ps with a Langevin thermostat. Following this, a 250 ps production run was conducted, from which dislocation positions were extracted (using OVITO dislocation analysis^79^) and an average dislocation velocity was calculated. In the all-expensive simulations, an NVE ensemble was used, whilst in the mixed simulations, a weak Nosé-Hoover thermostat was applied to the whole system (damping time 10 ps) in order to counteract the gradual heating that occurred due to energy flux from the ML/ML boundary. A short initial test was conducted to ensure that this thermostat did not impact dynamics.

The dislocation core was identified as the expensive region using Common Neighbor Analysis (CNA) in LAMMPS. The radius of the core, blending and buffer regions were 5.0, 4.0 and 5.0 Å, respectively. The expensive region corresponded to ≈ 1% of the total atoms in the simulation cell.

He implantation into W

For the study of He implantation into W, a modified version of the procedure outlined in Borovikov et al.^59^ was followed. We targeted normal deposition into 1000 K (001) W surfaces. A slab-type simulation domain was selected, which was non-periodic in z [001] and periodic in x [100] and y [010]. First, preliminary studies were run to determine the necessary size of the slab. It was found that high energy He atoms (> 40 eV) can penetrate very deeply into W samples, in some cases up to 4× the mean depth. It was also found that cells needed to be made wide enough that a deflected, obliquely penetrating He atom that passed through periodic boundaries did not pass close to atoms it had previously disturbed. Table 6 contains the cell sizes and number of atoms used for each He implantation energy. All simulations were carried out in LAMMPS, with the following procedure. Firstly, the slabs were thermalized to 1000 K over 10 ps of MD, with a 1 fs timestep using a Langevin thermostat with time constant of 0.01 ps. For computational efficiency, the same initial thermalized slab configuration was used for sets of implantations carried out in parallel. In all following simulations (both all-expensive and mixed), the thermostat was removed and the molecular dynamics was run using an NVE ensemble.Table 6. The cell size and number of atoms used for He implantation simulations at energies up to 100 eVHe energy (eV)Cell size (Å)Number of atoms5–3039.8 × 39.8 × 46.1507035–4046.1 × 46.1 × 62.0900050–6046.1 × 46.1 × 93.813,5007046.1 × 46.1 × 125.618,0008046.1 × 46.1 × 157.422,5009046.1 × 46.1 × 189.227,00010046.1 × 46.1 × 221.031,500

To generate the initial configurations for each implantation, a He atom was placed in a uniformly sampled random position 5.0 Å above the W surface and given a velocity in − z corresponding to the investigated kinetic energy (between 5 and 100 eV). An initial simulation timestep was selected corresponding to an initial He per-timestep displacement of 0.02 Å. The initial surface impact simulation was then run for 750 timesteps, which would correspond to a total straight-line He distance of 15.0 Å. At the end of this, the He atom position was measured. If it was found to be more than 2.0 Å above the surface, and with a positive z velocity, the simulation was stopped with the result being that the atom reflected. If not, the timestep was updated to either 0.02 Å divided by average He velocity or 0.1 fs, whichever was smaller, and a further 0.05 ps was run. This was repeated until either the He escaped the surface or a total time of 1 ps passed. For all energies, 1 ps was enough time for the He atom to thermalize.

All expensive simulations were carried out with He energies up to 40 eV. For each energy, 4 sets of 25 implantation attempts were carried out, with each set of 25 sharing the same thermalized slab structure.

For the mixed simulations, the expensive region was set to track the He atom, with a core radius of 6.0 Å a blending radius of 4.0 Å and a buffer radius of 8.0 Å. The thermalization was entirely conducted using the cheap potential. For all energies, 5 sets of 100 repeats were used.

The potentials were evaluated in LAMMPS on 40 GB A100 GPUs using the Kokkos^63,64^ accelerated variants of the symmetrix package for MACE^80^ and ML_SNAP package for SNAP^46^. In order to conduct mixed simulations where both the expensive and cheap potentials were evaluated on the same GPU, a Kokkos version of ML-MIX was developed. This can be found in v0.3 of ML-MIX onwards, and allows for force-mixing on GPU of many Kokkos-compatible LAMMPS pair_styles.

Supplementary information

Supplementary Information