Deep potentials for materials science

1.   Introduction

Atomistic simulations are playing an increasingly important role in materials science and changing how research in this heavily-experimental field is conducted [1]. All atomistic simulations in materials modelling require the input of some form of potential energy surface (PES) to describe how atoms interact; from which atomic forces are determined. The most accurate way to obtain the PES, within the Born-Oppenheimer approximation [2], is by solving the Schrdinger equation based on a quantum mechanical treatment of the electronic structure for fixed atomic nuclei coordinates [3]. However, in the most-widely applied electronic structure approach, density function theory (DFT) [4], the scaling is normally $\propto N^3$ where N is the number of atoms. This scaling makes DFT very costly in applications to large materials systems (over 1000 atoms) and for long simulation times (e.g. nanoseconds in molecular dynamics, MD). A widely applied solution to this efficiency problem is to develop empirical interatomic potentials [5-10], which describe the relationship between atom positions and system energy by assuming an analytical functional relationship, often based upon physical and chemical insights. Although large-scale, long-time atomistic simulations may readily be performed using empirical interatomic potentials, the accuracy of the atomistic simulations in describing real materials systems is often limited by the assumptions inherent to these empirical descriptions. In this context, one faces the dilemma that quantum mechanics methods are highly accurate but extremely inefficient for such atomistic simulations while empirical interatomic potentials are efficient, but commonly of limited accuracy.

Many approaches have been proposed to strike a balance between accuracy and efficiency in atomistic simulations. One approach to overcoming the low efficiency of widely-used DFT methods was the development of the ONESTEP program [11] in which plane-wave DFT calculations are performed on parallel computers, leading to computational costs that are linear in the number of atoms. DFT calculations have also been implemented on GPU machines [12-15] leading to accelerations by a factor of more than 20 compared to CPU machines [15]. The accuracy of empirical interatomic potentials can be improved, to some extent, by developing analytical functional forms that include many additional parameters to account for additional physical insights. A typical example is the modified embedded atom method (MEAM) potential [10] which extends the embedded atom method (EAM) [9] by considering the angular nature of electron density distributions. However, this tradeoff leads to a decrease in computational speed of the MEAM potentials compared to the simpler EAM potentials; such angular potentials are also more difficult to accelerate than simpler potentials. Nonetheless, the progress from these two approaches are encouraging and have led to many applications. Nevertheless, as materials systems of interest become larger, more complex and demands on predictability more severe, obtaining a better balance between accuracy and efficiency is both an urgent and challenging problem. This challenge requires fundamentally new approaches, rather than tweaking existing techniques. New insights may be garnered from recent progress in other disciplines.

Machine learning (ML) is well known for its surprising successes in, for example, pattern recognition and, as a result many different ML methods were developed in recent years [16-18]. ML potentials may be viewed as versatile descriptions of the PES parameterised using a flexible ML-based analytical form. The flexibility of ML potential increases the representability (accuracy) compared to empirical interatomic potentials and the analytical form significantly improves the efficiency relative to DFT calculations. In this perspective, ML potentials, fit to DFT results, have the potential for achieving DFT accuracy and empirical interatomic potential efficiency. Developing such ML potentials is a challenge.

Since the pioneering work of Blank et al [19] where ML neural network methods for describing the PES were first introduced, a variety of ML potentials have been proposed. Behler and Parrinello (BP) [20] introduced ML neural network potential (NNP) approach in which radial and angular symmetry functions are used as atomic environment descriptors; NNP has found applications in bulk silicon [20], carbon [21], TiO₂ [22] and many other materials [23]. BP-NNPs were categorised into four generations [23] for which there are several recent reviews [23-25]. Schtt et al [26, 27] developed SchNet and the SchNetPack package based upon a neural network framework to model the chemical properties and PES of molecular materials. Many other NNPs exist and are popular for different material systems [28-31]. Apart from NNP and SchNet, several other types of ML potentials were introduced; Gaussian approximation potentials [32-35], moment tensor potentials (MTP) [36-38], spectral neighbour analysis potentials [39-41], gradient-domain ML [42, 43], etc [44]. A comprehensive comparison of the major ML potentials in terms of accuracy and efficiency can be found in a recent paper [45].

Deep potentials (DP) [46-49] are of the NNP type that were first published in 2017. This approach has been used extensively for different material systems. The underlying theory has also developed continuously pushing these potentials to increasingly favourable combinations of accuracy and efficiency. Recently, DP has been applied to MD simulations of more than 100 million atoms with ab initio accuracy on a state-of-the-art supercomputer [50]. This is a good example of the power of integrating physical modelling, ML and high-performance computing. In this review, we focus on the application of DP in materials science and discuss a vision for future DPs. The paper is organised as follows. In section 2, we review the basic theory underlying the DP method (section 2.1), demonstrate the steps for developing DPs and their application for atomistic simulation (section 2.2), introduce the extant software for DP development (section 2.3) and the DP Library (section 2.4), and discuss how to make DP more practical for atomistic simulations (section 2.5). Section 3 lists many examples of DP applications in materials science, covering elemental bulk systems (section 3.1), multi-element bulk systems (section 3.2), aqueous systems (section 3.3), and other applications (section 3.4). Next, we discuss the efficiency and accuracy of DPs in practice (section 4.1) and a comparison of the computational speed of DP versus other approaches (section 4.2). We conclude section 5 with an assessment of where DP is going in the near future.

2.   Deep potential

2.1   Theory

The local (neighbourhood) dependence of the atomic energy is an assumption. There are, however, non-local (long-range) interactions, arising mainly from Coulombic interactions within the electron density distribution. For metallic systems, the local dependence assumption is reasonable as a result of shielding effects. For homogeneous materials, the long-range interactions attenuates rapidly with increasing atomic separation such that a sufficiently large r_c can always satisfy the local dependence assumption. For materials where long-range interactions dominate, these must be explicitly considered in the model construction. Although there are research reports which introduce long-range interaction in the model construction [28, 30, 51-59], there is still no widely accepted method which can handle this interaction appropriately. In the DP method, we focus on the most common situation, in which $E_i^w$ is well-described by the local interaction assumption. Also note, that when the system is extended (rather than molecules), we explicitly involve periodic boundary conditions in the usual way; for more details, see, e.g. [60].

Different types of ML potentials employ different descriptors. Some examples are the symmetry functions from BP [20], the Smooth Overlap of Atomic Positions from Bartok et al [34], the SchNet descriptors’ from Schtt et al [26, 27], and the moment tensor from Shapeev [36]. In the following, we introduce the DNN and then focus on the construction and physical meanings of the descriptors employed in DPs.

2.1.1   Deep neural network.

Deep neural networks (DNN) and kernel methods are the two dominant methods employed in constructing ML potentials. Kernel-based methods usually construct descriptors explicitly, while DNN methods can be used to learn’ descriptors. These are both broad areas, so we refer interested readers to a standard reference [18]. Due to the robust fitting ability of DNN for high-dimensional and nonlinear properties, the DP uses DNNs both as a fitting method for ${F}$ and also to construct descriptors of local atomic environments (as described below, section 2.1.3)

Quantification of the representability of a DNN is an active area of research. Barron [62, 63] proved that a neural network with only 1 hidden layer (P = 1) and an arbitrary number of neurons can approximate a class of functions with arbitrary precision. Many researchers have subsequently investigated the approximation ability of DNNs with P>1 [64-67] and explained why the DNN is more successful (and widely used) than wide neural networks (many neurons in one hidden layer). E et al [68, 69] explained and demonstrated why DNNs are suitable for high dimensional problems.

2.1.2   Descriptors for the non-smooth DP.

A schematic of the descriptor for the non-smooth DP is shown for a water molecule in figure 1. The red and blue spheres denote oxygen and hydrogen atoms, (e_x, e_y, e_z) is the local frame of atom i and its neighbour atom j, and $\left(x_{ij}^{\mathrm{\prime }},y_{ij}^{\mathrm{\prime }},z_{ij}^{\mathrm{\prime }}\right)$ are the Cartesian components of ${\mathit{r}}_{\mathit{i}\mathit{j}}$ in this local frame. The descriptor ${{D}}_{ij}$ may contain only radial (distance) information and/or angles. In the water molecule example, neighbours of atom i are first sorted according to their chemical species (oxygen first, then hydrogen). Then within the same species, the inverse distances to atom i are used to sort the atoms. Finally, ${{D}}_{ij}$ is applied to the sorted input data for atom i.

Figure  1.  Schematic of the non-smooth DP descriptor for a water molecule. The red and blue spheres denote oxygen and hydrogen atoms, respectively.

DownLoad: Full-Size Img PowerPoint

The advantage of non-smooth descriptor is that all neighbour information is preserved. However, due to uncertainty in the choice of neighbour atoms a(i) and b(i), the descriptor is non-smooth. In practice, a(i) is picked as the nearest neighbour and b(i) as the second nearest neighbour. Continuous change in atom positions can thus result in a discontinuous change of the atom number, the local frame, and the local coordinates. In addition, the sorting operation for the other neighbours introduces additional discontinuities in the descriptor and its derivatives (see equation (18)).

2.1.3   Descriptors for the smooth DP.

Figure  2.  Workflow for constructing the smooth DP descriptors.

DownLoad: Full-Size Img PowerPoint

2.2   Developing and applying DPs

After illustrating the descriptors and DNN fitting in the DP method, we now outline the main steps we apply to develop a DP: preparing training datasets, the training process, and model validation. We also explain how to apply DPs in atomistic simulations.

2.2.1   Preparing training datasets.

The preparation of the training data for a DP has two parts: (1) providing atomic configurations, (i.e. the coordinates of atoms and the cell shape tensor) and (2) labelling (i.e. calculating the energy, atomic forces and virial tensor for the configuration). Labelling is done through a DFT calculation; e.g. well-established DFT packages such as the Vienna Ab initio Simulation Package (VASP) [71, 72], Quantum ESPRESSO (QE) [73], and Atomic-orbital Based Ab-initio Computation at UStc (ABACUS) [74]. Since DPs are trained using a DFT training set, the accuracy of a DP will never exceed that of its training data. The errors in DFT calculations have two main sources. The first is the error introduced by the approximate form of the exchange-correlation functional. This type of error may be reduced by moving up the Jacob’s ladder [75] of increasingly accurate exchange correlation functionals (usually at higher computational cost). It is also possible to use the high order post-Hartree Fock methods like the Mller-Plesset perturbation [76], coupled cluster [77], and configuration interaction [78] for labelling. The other source of error in DFT calculations is numerical, i.e. error introduced by the numerical discretisation of wave functions in real and k spaces and convergence. This type of error can be systematically controlled with the use of more complete basis sets (increasing the energy cutoff in plane-wave approximations), reducing k-space mesh spacing and using stricter stop criteria for self-consistent field iterations. Generally speaking, increasing label quality implies larger computational demands. One usually seeks for a balance between quality and cost.

The training dataset is another critical issue for generating a DP; here, the two main issues are (1) completeness and (2) compactness. By completeness, we mean the training datasets need to sample the relevant configuration space as completely as possible. Increasing the diversity in the members of the training dataset helps increase the transferability of the DP. By compactness, we mean the training data should be the minimal subset of the sampled configurations from which a model with uniform accuracy on the sampled configurations is trained. This is important for minimising DFT computational time. Different approaches have been used to sample the configuration space including MD simulations, genetic algorithms [79, 80], enhanced sampling methods [81], active learning [82] and concurrent learning schemes [83, 84]. Among these, the concurrent learning scheme we refer to as the deep potential generator (DP-GEN) is found to be particularly effective at generating training datasets that satisfy both the completeness and compactness conditions (see section 2.3.3).

2.2.2   Training process.

We employ the DeePMD-kit package [49] to train a DP. A python package, dpdata (web link 1 in section S2 of the supplementary material (SM) available online at stacks.iop.org/MF/1/022601/mmedia) is used together with the DeePMD-kit to transform the DFT labels from DFT software to the data format accepted by the DeePMD-kit package. DeePMD-kit and dpdata package details are discussed in section 2.3.

There are several important issues in the training process. First, the DP should be trained such as to avoid under-fitting and over-fitting. Under-fitting implies that the DP performs poorly for both the training and validating datasets. Increasing the number of fitting parameters (wider or deeper DNN) or adjusting the ML algorithms can help avoid under-fitting. Under-fitting is easily detected in the training process. On the other hand, over-fitting suggests that a DP is very good at reproducing the training data, but poor at prediction. This may be solved by increasing the training dataset size, decreasing the number of fitting parameters, or using the force and virial as labels (see below). Over-fitting is more difficult to identify; hence, comprehensive testing is necessary.

Second, while training may involve the energy, force and virial labels (normally from DFT calculations), not all of these are necessary. However, it is strongly recommended to use force labels because (a) force has more information: 3N vs. 1 compared to the energy label and (b) training with gradient information helps avoid over-fitting. The schematic plot in figure 3 demonstrates that the gradients of an over-fit model deviate from the ground truth, thus training with forces (gradients of energy with respect to coordinates) and virials (gradients of energy with respect to the cell tensor) as labels helps avoid over-fitting.

Figure  3.  Schematic depiction of an over-fit ML model (a) and a properly fit ML model (b). The ground truth and ML model are denoted by the red and blue lines. The training data are denoted by red . The over-fit values are close to the ground truth for the training data, while the gradients deviate from the ground truth (red and blue arrows). Both the values and gradients of the properly fit model are close to the ground truth for the training data.

DownLoad: Full-Size Img PowerPoint

Third, there are many tunable hyper-parameters in the training process in the DeePMD-kit, including neural network size, learning rate, prefactors in the loss functions (equation (3)). In practice, we observe that the quality of the DP is not very sensitive to these hyper-parameters and the default settings in the DeePMD-kit routinely provide reasonable accuracy. In some cases, high accuracy for energy is required and fine tuning’ with larger energy prefactors helps. A more detailed discussion of the hyper-parameters in DeePMD-kit is in section 2.3.2.

2.2.3   Model validation.

After the training process is complete, it is advisable to validate the obtained DP to decide whether additional training datasets are required prior to atomistic simulation applications. There are two main approaches to test the DP performance. (1) We can construct a small set of datasets which are unrelated to the training datasets. The DP is used to predict the total energy, atomic forces, and virial tensors of the constructed datasets and the results are compared with DFT results for the same atomic configuration. Experience shows that the root mean square error (RMSE) for energy and forces for a good DP should be smaller than 10 meVatom^-1 and 100 meV Å^-1. It is not uncommon to achieve 1 meVatom^-1 and $<$50 meV Å^-1 errors for energy and forces. (2) The DP can be applied to calculate properties of interest to be compared with DFT. For example, elastic constants are important properties for structural applications that can be obtained both from DFT and DP MD. Application of DPs in atomistic simulations are discussed below and the Autotest package in DP-GEN conveniently calculates several different properties (see section 2.3.3).

2.2.4   Model inference.

Model inference (a commonly used term in the ML community), is the process of providing live data to the ML model to obtain an output. In the context of ML models for the potential energy surface, the inference means taking the configuration (atom coordinates and cell tensor) as input to calculate the energy, force and virial tensor. With the interfaces provided by the DeePMD-kit package, one can easily make inferences about the DP in a Python or C++ programming environment. This makes it possible to use DPs in various molecular simulation tasks such as MD, Monte Carlo, geometric optimisation, by interfacing DeePMD-kit with molecular simulation packages, such as the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [85], the atomic simulation environment (ASE) [86], i-PI [87], and GROMACS [88].

2.3   Software for DP development

There are three main packages for DP development: (1) dpdata for converting the output of DFT software to the data format accepted by DeePMD-kit; (2) DeePMD-kit for DP training and inference; (3) DP-GEN for efficient sampling and labelling of the training data. The DP-GEN package integrates with DeePMD-kit and DFT calculation packages to automatically generate and test DPs.

2.3.1   dpdata.

dpdata is a python package which converts and manipulates training data in different DFT package formats to the compressed format used by DeePMD-kit. A typical dpdata workflow is as follows. (1) Load data from data files. Data files can be written in one of the following package formats: VASP [71, 72], LAMMPS [85], Gaussian [89], SIESTA [90], CP2K [91], QE [73], FHI-aims [92], QUIP (web link 2 in section S2 of the SM), PWmat [14, 15], AMBER [93], GROMACS [88], and ABACUS [74]. (2) Data may be manipulated through operations including replication of atom configurations in a supercell, perturbation of the cell vector and atom positions, and replacement of a number of one type of atoms with others. (3) Output data in a format of one of the aforementioned software. Here we give a simple example of transforming the OUTCAR file from VASP into a training dataset that DeePMD-kit can read:

1 import dpdata

2 d_outcar = dpdata.LabeledSystem("OUTCAR")

3 d_outcar.to("deepmd/npy", "dpmd_raw")

After execution of these commands, a directory named dpmd_raw is created that stores the numpy compressed data appropriate for DeePMD-kit use. The dpdata package source code and manual is available at web link 1 in section S2 of the SM.

2.3.2   DeePMD-kit.

The DeePMD-kit package [49] (first made publicly in November 2018) is available on GitHub (web link 3 in section S2 of the SM) under the GNU Lesser General Public License (LGPL). DeePMD-kit has continued to evolve since first becoming available. The DeePMD-kit package interfaces with TensorFlow [94] to make the training and inference codes more efficient and automated. The training is integrated and made available via a command line interface. Model inference is provided through C++ and Python interfaces, which accept atom positions and the cell tensor and return the energy, force and virial for this configuration. These interfaces may be used by MD and molecular simulation packages written in C, C++ or Python, as discussed in section 2.2.4. The DeePMD-kit supports GPU accelerated training and inference. When interfaced with the LAMMPS MD package, parallel and distributed computations accelerated by GPUs are available. (On Summit, one of the most powerful supercomputers in the world, DeePMD-kit has pushed the limit of MD with ab initio accuracy to 100 million atoms and achieved a peak performance of 91 PetaFLOPS in double precision (45.5% of the theoretical peak) [50]). The DeePMD-kit package is now well-developed, been installed over 30000 times (GitHub) and received over 650 GitHub stars.

DeePMD-kit can now be easily installed via off-line packages and package managers such as Conda and Docker. A detailed introduction can be found at web link 4 in section S2 of the SM. The brief release history and key milestones of DeePMD-kit can be found in the SM.

2.3.2.1   Training and testing of DPs.

After installation of DeePMD-kit package, we follow the path in section 2.2 to develop a DP. First, dpdata is used to generate the training datasets to be read by DeePMD-kit package (see section 2.3.1). Then, DP training is started with the following command:

1 dp train input.json where input.json is the input script controlling the training process. If the training program stops, it can be restarted with

1 dp train-restart model.ckpt input.json where model.ckpt is the checkpoint file storing the existing model and training status. A detailed explanation of the output files and training (validation) errors can be accessed at web link 5 in section S2 of the SM.

After the training process, the architecture and the parameters of the DP can be abstracted from the checkpoint file and saved in the DP file in the Google’s protobuf format

1 dp freeze -o graph.pb where graph.pb is the DP file. The DP file can be used for testing and model inference via the Python and C++ interfaces. For example, the DP can be used to give energy, forces, and virial errors on a designated datasets:

1 dp test -m graph.pb -s /path/to/system -n 30 where the tested model follows -m tag, -s gives the path to the tested system, and -n is the number of frames in the tested system.

2.3.2.2   Fine-tuning a DP.

The adjustment of hyper-parameters in ML is never a trivial issue. Fortunately, the efficiency of training a DP is not very sensitive to hyper-parameters and the default settings in DeePMD-kit usually yield a DP of reasonable accuracy. Detailed documentation of all available training parameters can be accessed at web link 6 in section S2 of the SM. On occasion, the hyper-parameters settings may not be satisfactory; e.g. if one wants a DP with higher energy and virial accuracy without a large penalty in force accuracy, fine-tuning of the hyper-parameters may be necessary.

Model fine-tuning restarts training. In the fine-tuning, the model parameter are initialised from the DP trained with the default setting, the start learning rate is decreased to 10^-4 while the stop learning rate is kept, and the loss prefactors set to $p_e^{\text{start}}=10$, $p_e^{\text{limit}}=100$, $p_f^{\text{start}}=1$, $p_f^{\text{limit}}=1$, $p_v^{\text{start}}=10$, $p_v^{\text{limit}}=100$. The idea behind this setting is to inherit model parameters from a DP with high force accuracy, and to train the model with significantly larger energy and virial prefactors to focus the training on the energy and virial accuracy. The learning rate is decreased to slow the forgetting of the model parameters of the original model and to preserve the force accuracy. If the numerical error in the labels is low, then the fine-tuning usually leads to a significant improvement on the energy and virial accuracy with little decrement to the force accuracy [95].

Another hyper-parameter critical to the accuracy of the DP is the number of training steps. Different settings may be employed for different DP purposes. In the concurrent learning procedure, (for example DP-GEN see section 2.3.3), 0.4-2 million training steps should suffice. If the DP is intended for high-accuracy atomistic simulations, over 8 million training steps are often employed for from-scratch training or fine-tuning.

Finally, the atom_ener parameter may be used to specify the energy of an atom in the vacuum for each atom. For metals and alloys, the cohesive energy calculated from DFT does not always match experiment. For example, for Ti, the cohesive energy calculated from DFT is 5.34 eVatom^-1 [96], while experiment gives 4.85 eVatom^-1 [97]. Because DP is fitted to DFT data, this command helps to correct the DFT errors’ for the isolated atom. Accurate cohesive energy relative to experiment is critical for DP in some applications, such as the fracture behaviour of metals and alloys.

2.3.2.3   Interfacing with third-party packages.

After a DP is trained and tested by DeePMD-kit, application requires linking with other software. DeePMD-kit provides Python and C++ interfaces for model inference, which is helpful for calculating energy, atomic forces, and virial tensors with input atomic coordinates and the cell tensor. DeePMD-kit also interfaces with ASE [86], LAMMPS [85], i-PI [87], and GROMACS [88] for DP-based atomistic simulations. The detailed use of DP for these can be found at web link 7 in section S2 of the SM.

Here, we explicitly discuss using DPs in LAMMPS. The installation package for DeePMD-kit already incorporates the latest stable version of LAMMPS. Running MD or other molecular simulations with DP in LAMMPS is very simple. The user only needs to add two lines of commands specifying the interatomic interaction with the LAMMPS command pair_style and pair_coeff to the LAMMPS input script:

1 pair_style deepmd graph.pb

2 pair_coeff where graph.pb is the DP file. The user need not change any other part of the input script to conduct DP-based molecular simulations in LAMMPS.

2.3.3   Deep potential generator.

Deep Potential GENerator (DP-GEN) is a software package that implements the concurrent learning framework for generating high quality DPs.

2.3.3.1   Basic concepts and framework.

As discussed in section 2.2.1, generating a complete and compact dataset is critical for training high quality DPs. A straightforward way to generate data is to run finite-temperature ab initio MD (AIMD) simulations and use the configurations along the trajectories with labels (energy, forces, and virial tensors) as training datasets. However, this method of data generation is not particularly efficient. First, AIMD simulations are computationally expensive since labels are calculated at each time step by DFT. Second, computationally affordable AIMD simulations are typically very short and do not effectively explore phase space. Many important phenomena, such as phase transformations are difficult to be observed on AIMD time scales. Also configurations generated at successive MD time steps are extremely similar to one the other.

The DP-GEN concurrent learning scheme is designed to overcome these difficulties by generating a much more complete and compact training dataset for DP training. The DP-GEN is an iterative scheme. In each iteration it performs exploration, labelling and training. The training dataset is gradually enriched, and the quality of the DP improves on each iteration. The DP-GEN scheme is deemed converged when all relevant configurations are explored and the DP is uniformly accurate on the explored configurations. The three steps in each DP-GEN iteration are as follows.

(a) Exploration. We start with an ensemble of DPs. One of these DPs is used in the exploration step (figure 4(b)) to efficiently sample the relevant configuration space. The sampler typically performs several DP MD simulations (different initial configurations) for a set of thermodynamics conditions. In principle, sampling can also be performed using Monte Carlo simulations, enhanced sampling MD simulations [98], and any molecular simulation method that explores configuration space.

Figure  4.  Schematic illustration of the DP-GEN scheme. (a) The DP-GEN scheduler runs iteratively, performing three steps: exploration, labelling, and training. (b) In the exploration step, different structures are sampled using DP MD and an error indication is applied to choose the candidates for labelling. (c) In the labelling step, DFT calculations are performed on the candidates to obtain the energy, forces, and virial tensor. (d) In the training process, DeePMD-kit package is used to train a new DP based on the initial datasets and candidates from each iteration. Reprinted from [84], Copyright (2020), with permission from Elsevier.

DownLoad: Full-Size Img PowerPoint

The ML potential models cannot be extrapolated’ to the configurations that they are not trained on [34], thus the sampling should generate training datasets that explore the relevant configuration space as completely as possible. In the context of the DP-GEN scheme, the configurations in the training dataset are proposed in the exploration step, thus it should explore the relevant configurations as complete as possible by exploiting the high efficiency of the DP itself. The design of the exploration strategy depends on the applications for which the DP is intended. For example, a DP for liquid water properties need not fully explore ice configurations. A DP for chemical reactions may need enhanced sampling techniques to explore reaction pathways usually not available by standard MD simulations [99].

How the model deviation indicates the error is shown in the schematic figure 5. Four deep learning (DL) models (potentials), DL 1 to 4, are trained from the same training data (blue crosses) with random initialisation of the model parameters. At times, DL1 to DL4 all show small training error in the region covered’ by the training data, but the predicting error is large in the region far away from the training data. The standard deviation of the DL model predictions follows exactly the same trend. In the region covered’ by the training data, all models/potentials are trained against the same target, thus their predictions agree with each other, while in the region far from the dataset the different DL potentials (training from different initialisations lead to different potentials) yield different predictions. Therefore, the model deviation is a good indicator of the DL prediction error.

Figure  5.  Visualisation of the model deviation. DL 1 to 4 indicates 4 Deep Learning models trained from the existing data (blue crosses).

DownLoad: Full-Size Img PowerPoint

Upper and lower trust levels, _hi and _lo, are used to select configurations for labelling. When the error indicator () is smaller than _lo, it suggests that the atomic forces of all atoms in the configuration are accurately predicted and there is no need to do DFT labelling for this configuration. When $ϵ⩾{\sigma }_{hi}$, the potential behaves poorly and DFT labelling will not be performed, because such explored configurations may be nonphysical (e.g. overlapping atoms) due to poor potential quality. Only when ${\sigma }_{lo}⩽ϵ<{\sigma }_{hi}$, will the configuration be selected as candidates for DFT labelling. The user may also set the maximal number of structures in the candidate set to perform DFT labelling; this leads to a randomly sampling to form the labelled dataset.

(b) Labelling. In the labelling step (figure 4(c)), a DFT software package is called to calculate the energy, force and virial tensor on the selected configurations. The labelled data are added to the existing training dataset.

(c) Training. In the training step (figure 4(d)), an ensemble of new potentials is trained based on the same training dataset but with different and random initialisation of the potential parameters. Relatively short trainings are performed at this step, because the configurations are selected for labelling by the force criteria, and reasonably accurate DP forces can be obtained with relatively short trainings. Experience shows that 400000 to 2000000 training steps are appropriate. One may also initialise the potential parameters in the ensemble of potentials from the previous iteration, and start the training in the next iteration with a relatively small initial learning rate (e.g. 10^-4), and with a bias to the new training data.

(d) Initial datasets for starting the DP-GEN workflow. The DP-GEN workflow introduced above is automatic, but human intervention is required at the beginning when no DP exists. Generally, we generate some initial datasets and train an ensemble of initial DPs to kick off the workflow. In metals and alloys, the generation of initial bulk datasets usually consists of the following steps. (a) Start from a supercell ($2\times 2\times 2$) of crystal unit cells (BCC, FCC, HCP). (b) Compress and dilate the supercell uniformly to cover a range of densities. (c) Randomly perturb the atom coordinates and supercell vectors. (d) Randomly replace atom elements to generate alloys. (e) Start from different structures and do a short AIMD at low temperature (100 K). (f) Collect the DFT results and generate initial training datasets. The initial datasets are not particularly important as long as they reasonably cover the starting structures of DP MD in the exploration step to avoid breakdown of the DP MD within the first few steps. This generating initial dataset method is simply one possibility; other methods such as generating initial structures according to crystal symmetry in genetic algorithm [80] are also applicable.

2.3.3.2   A practical guide.

Here we provide a practical guide for generating a general-purpose DP for a bulk metal or alloy using the DP-GEN scheme.

(a) Generation of initial datasets. Initial datasets are needed to generate the first ensemble of DPs. Because initial datasets are only used to start DP-GEN, datasets with only several hundred different configurations should suffice. The accuracy of the DPs is not sensitive to the choice of the initial datasets and would be gradually improved as iterations increase assuming that the DP-GEN run parameters are set appropriately.

(b) Exploration. The exploration strategy is not unique and is determined by the user according to the applications of interest. Here, we focus on the general exploration strategy using, again, bulk metals or alloys as our example. DP-based MD simulations at different temperatures and pressures are employed to explore the configurational space. For efficient sampling, we increase the temperature during the exploration. If the highest melting point of the element in an alloy is T_m, we usually divide the temperature range from 50 K to 2T_m into 4 regions. (a) [50 K, (0.1 , 0.2 , 0.3 , 0.4)T_m] (b) [0.5, 0.6, 0.7, 0.8, 0.9]T_m; (c) [1.0, 1.1, 1.2, 1.3, 1.4]T_m; (d) [1.5, 1.6, 1.7, 1.8, 1.9]T_m. For each temperature region, the pressure range is varied over [0.001, 0.01, 0.1, 1, 5, 10, 20, 50] kBar. In this case, there are $5\times 8=40$ different MD conditions in each temperature region. Other parameters for the DP MD simulations is the MD simulation duration and the number of different initial structures. At the beginning of each new temperature region, the number of MD step should be small (500 or 1000) and the number of starting structures should also be small (e.g. 5 each for different distorted crystal supercells). This is because the DPs have little knowledge’ of this new temperature region such that the expected accuracy is low. As the iteration number increases within a temperature region, the MD runs are longer and the number of starting structures is increased.

After the MD runs, candidate systems are selected based upon the trust levels _lo and _hi. The key idea is to make DP-GEN converge within each temperature region. If _lo and _hi are too strict, for example _lo is too low, the accurate ratio ($ϵ<{\sigma }_{lo}$) stays low (smaller than 90%) in the temperature region. If _hi is too high, too many structures are selected as candidates for labelling and the ratio of the failed structure ($ϵ⩾{\sigma }_{hi}$) is very low. The values of _lo and _hi used in our previous work for different systems (see table 1); _hi is normally 0.15-0.30 eV Å^-1 higher than _lo. In principle, _lo should be higher than the lowest error that a DP can achieve (error introduced by data, fitting ability of optimiser, representability of the DP). _hi should be set so as to eliminate unphysical configurations from the dataset, since the convergence of the DFT calculations for such configurations is difficult to obtain. See [84] for additional guidance.

Table  1.  Trust levels _lo and _hi employed in DP-GEN for several systems.

System	[lo,hi]
Mg [83]	[0.03,0.13]
Al [83] & Al-Mg [83]	[0.05,0.15]
Cu [84]	[0.05,0.20]
Mg-Al-Cu [100]	[0.05,0.20]
Ti [101]	[0.10,0.25] at $T<1.5T_m$a for bulk
	exploration and [0.15,0.30] elsewhere
W [102]	[0.20,0.35]
Ag-Au [103]	[0.05,0.20]
water [95]	[0.15,0.25] in first 24 iterations
	[0.18,0.32] in iterations 25 to 32
	[0.20,0.35] in iterations 33 to 36
SiC [104]	[0.15,0.30]
Li10(Ge,Si, or Sn)P2S12 [105]	[0.12,0.25]
Tm is 1941 K, which is the experimental melting point for Ti.

 | Show Table

DownLoad: CSV

(c) Labelling. After candidate structures are selected from the exploration step, users can choose the maximal number of structures sent to DFT labelling in each iteration. On the one hand, if there is a large number of structures for DFT labelling in one iteration, the computational cost for labelling would be high and this dataset of configurations may be redundant, leading to wasted computational resources. On the other hand, if the number of structures for DFT labelling is too small (e.g. one) in one iteration, the datasets after DP-GEN run would be very compact, requiring many iterations, which again would be computationally expensive for training in each iteration. In this case, the maximal number of structures for DFT labelling in one iteration should be set appropriately, we typically set this number to 100.

(d) Training. The prefactor of the virial tensor in the loss function (equation (3)) is usually set to zero in the DP-GEN training loop for two reasons. First, only the model deviation of the forces is used to select candidates in the exploration step; this is irrelevant for the virial tensors. Second, longer training is needed to obtain a DP with good accuracy on both forces and virial tensors. The training step in each iteration of DP-GEN should be small, so we do not initially train on the virial tensors.

(e) DP-GEN convergence. In principle, DP-GEN is considered converged when the model deviation () on each structure is smaller than the trust level _lo. However, there are often some (very low probability) configurations with undesirable accuracy in long enough MD simulations. These configurations generally do not improve the DP performance. In this case, DP-GEN allows users to set the convergence criterion at which level the DP is regarded as converged. For example, the user can skip the labelling step when the percentage of configurations with model deviation smaller than _lo is larger than 99%. After convergence, long training steps (usually over 8000000 steps) are performed on all DP-GEN loop datasets by including labelled energy, forces, and virial tensors. Comprehensive tests of the resulting model are then performed and the Autotest in DP-GEN package is applied to perform different property tests (see below). The DP obtained after careful DP-GEN run and long training is normally good for general purposes, but does not guarantee good performance for subtle properties. Specialisation for these properties is required and the details of specialisation are presented in section 2.5.2.

2.3.3.3   Autotest.

After one or a series of DPs is trained, we employ the Autotest package, to calculate a simple set of properties and/or perform tests for comparison with DFT and/or empirical interatomic potentials (EAM, MEAM, and etc). Because the DP was obtained by fitting DFT calculation results, the degree of agreement between DP and DFT for a series of properties is usually excellent; if not, this comparison provides a guide for further training and parameter settings. Autotest, as part of the DP-GEN package, has standardised the calculation of some critical physical properties and provides a series of reliable benchmark testing to better evaluate the performance of the DPs. As of this writing, Autotest includes the calculation of the following set of properties (additional properties are added continuously): (a) equilibrium structural parameters (relaxation), (b) equation of state (eos), (c) elastic constants (elastic), (d) vacancy formation energy (vacancy), (e) interstitial formation energy (interstitial), and (f) surface formation energy (surface).

Autotest can also use LAMMPS or VASP to reproduce and refine previous calculation results. The current package is mainly targeted for simple crystal structures (metals and alloys), but is extensible for users to implement new features. The detailed manual and framework of Autotest package is available at web link 8 in section S2 of the SM.

2.4   The DP Library

The DP-GEN package provides a relatively automatic routine for generating DPs, and the DP Library project is a place for sharing and publishing the DPs and the training datasets. First, DP Library is a place for model sharing, which is similar to other interatomic potential repositories like the NIST (web link 9 in section S2 of the SM) and OpenKIM (web link 10 in section S2 of the SM). Second, DP Library provides an opportunity for data sharing. One can contribute and download DFT datasets used to train a published DP. If a DP needs to be refined or improved, one can first add new data to the downloaded training datasets from DP Library and then perform training using DeePMD-kit. For example, one may develop a DP for A-B alloys based on the training data for element A and element B, instead of generating the datasets by him/herself. The settings (usually settings to use a DFT software) used to generate the dataset are asked to provide, so the new data can be generated in the same way as the downloaded dataset. The shared training data makes the published DPs reproducible and improvable. DP Library can be accessed at web link 11 in section S2 of the SM and an image of the website is shown in figure 6. In the periodic table, the available DP for elements is in black and details about the DP and DFT data can be found by just several clicks. For more details, please refer to the manual of DP Library at web link 12 in section S2 of the SM.

Figure  6.  Elements for which DPs are currently available in the DP Library (web link 11 in section S2 of the SM) are indicated in black.

DownLoad: Full-Size Img PowerPoint

2.5   Efficiency and accuracy of DPs for applications

To make DP more practical for applications in different materials systems, efficiency and accuracy must be balanced.

Although DP is very much faster than DFT, it is still much slower than empirical interatomic potentials like EAM and MEAM. This is not surprising given the vast number of parameters in the neural net of the DP (compared to simple empirical potentials such as EAM). Hence, optimisation of the efficiency in use of the DP can greatly enhance the practical usability. Recently, a highly optimised GPU version of DeePMD-kit pushed the limit of MD with ab initio accuracy to 100 million atoms [50, 106] (recognised by the 2020 Gordon Bell Prize). Further optimisations are still possible based on advances in neural networks.

DPs generated by the example exploration protocol (section 2.3.3) may not be accurate for all applications. For example, the Ti DP yields screw $⟨\mathbf{\text{a}}⟩$ dislocation properties that are not consistent with experiment (or DFT). For some particular applications, the DP must be tuned or specialised for important, yet subtle properties. This may be accomplished through specialisation, as discussed below [101].

We first introduce DP model compression which can easily accelerate DP by a factor of 4 to 18, based on experience. Next, we discuss DP specialisation and consider the example of Ti, mentioned above.

2.5.1   DP compression.

The most computationally intensive part of using the DP is the evaluation of the embedding net (equation (23) or (24)) and the assembly of the descriptor by equations (26)-(28) [107]. The goal of DP compression is to reduce the computational and memory overhead associated with the embedding net with little loss of accuracy. The evaluation of the DP involved the mapping from a scalar to a vector of dimension M, each dimension of which can be approximated by a piece-wise fifth order interpolating polynomials. The range of the embedding net is first discretised by nodes $x_1,x_2,\dots ,x_l,\dots ,x_{L+1}$, and the lengths of intervals (tabulation step) are assumed to be the same. On each interval, e.g. $[x_l,x_{l+1})$, the 6 fitting parameters of the fifth-order polynomials are uniquely determined by matching the value and the first and second derivatives of the polynomial to those of the embedding net on nodes x_l and $x_{l+1}$. This approximation is referred to as tabulation, because the parameters of the polynomials are stored in a table, and use involves simply a look-up the table for the polynomial when evaluating the embedding net. The accuracy of the tabulation is controlled by the tabulation step; a size of 0.01 leads to errors in energy and force smaller than 10^-7 eVatom^-1 and 10^-6 eV Å^-1. After tabulation, the multiplication between the embedding and the environment matrices becomes the bottleneck. The embedding matrix is stored in the memory after it is computed and then loaded from memory to the register for the matrix multiplication. This requires significant I/O with the memory. This can be eliminated by merging the tabulation and the multiplication; i.e. once one component of the embedding net is computed by tabulation, it is immediately multiplied with the environment matrix pre-loaded in the register, and accumulated to the result. With this optimisation, the multiplication between the embedding net and the redundant zeros in the environment matrix can be avoided.

The DP compression was benchmarked for Cu, H₂O and Al-Cu-Mg ternary alloy DPs [107]. The model inference was accelerated by 9.7, 4.3 and 18.0 times on a CPU and by 9.7, 3.7 and 16.2 times on an Nvidia V100 GPU for Cu, H₂O and Al-Cu-Mg, respectively. The maximum number of atoms handled by one GPU also increased from 12, 49 and 5 10³ to 129, 246 and 61 10³, respectively.

Figures 7(a) and (b) display the speed comparison of the compressed Ti DP [101] with an EAM [108], and an MEAM potential [109] on a CPU and GPU machine. Note that the DP has a larger radius cutoff distance than EAM and MEAM in this case. On CPUs, the compressed DP is 200-300 times slower than EAM potentials and 30-40 times slower than the MEAM potential. On GPUs, the compressed DP is 20-30 times slower than the EAM potential (MEAM is currently not ported to GPU in LAMMPS). All potentials show a linear scaling with the number of atoms. Because of this linearity and speed, the compressed DP can be used to perform large scale MD simulations to investigate a wide range of properties with ab initio accuracy; e.g. defect properties, phase transformations, and transport coefficients.

Figure  7.  A comparison of the speed of MD simulations using the compressed Ti DP, an EAM, and/or an MEAM potential on (a) CPU and (b) GPU systems [101].

DownLoad: Full-Size Img PowerPoint

DP compression is supported in the DeePMD-kit package (releases beyond 2.0.0) and the compressed model can be easily generated using the following command:

1 dp compress -i graph.pb -o graph-compress.pb

Using the optional -s flag, followed by the tabulation step with the default value of 10^-2 typically gives very accurate compressed DPs (compared with the original). It is recommended that the user compare the values of a few key properties from the original and compressed DPs.

2.5.2   Specialisation.

The DP generated from the DP-GEN scheme described above may not be sufficiently accurate for some complex phenomena and specialisation of this general-purpose DP may be required. This is not surprising since the exploration strategy may not provide a sufficient sampling of the relevant local structures that are inherent to the complex phenomena of interest. Therefore, some special structures should be added to the training process to better represent the requisite subtle properties. The workflow for specialising DP is shown in figure 8. The initialisation and DP-GEN loop steps were discussed in section 2.3.3.

Figure  8.  The workflow for specialising a general DP. R_i is the atomic coordinate of atom i, E is the total energy of one configuration, V is the virial (stress) tensor of one configuration, f_i is the force on atom i, n is the number of atoms in one configuration, DP¹ is the first ensemble of trial DPs, labels the ^th DP in the ensemble, is the standard deviation, and _lo and _hi are two thresholds in DP-GEN. The specialisation step is shown in the green box. Reprinted by permission from Springer Nature Customer Service Centre GmbH: [Nature] [npj Computational Materials] [101].

DownLoad: Full-Size Img PowerPoint

In the specialisation step, Special’ structures are first created based on the special properties/structures of interest. For example, if the DP elastic constants are not sufficiently accurate, Special’ structures may be created corresponding to a range of different crystal deformations. If the specialised DP is to describe dislocation properties in complex crystal, then Special’ structures that include sheared structures akin to those used in determining the generalised stacking fault energy -lines [110] may be of use. DFT calculations are then performed on these Special’ structures and the energies, atomic forces, and virial tensors along with the configuration form the Special’ training sets. The Special’ training sets then combine with the Classic’ training sets from DP-GEN loop and weights are chosen for the Special’ training sets. The default weight is 1 but this weighting of the Special’ training sets should be increased because (1) the number of Special’ training sets are usually much smaller than those in the Classic’ training sets, and (2) experience shows this may be necessary to properly reproduce the special properties. Next, DeePMD-kit package is applied to retrain the DP on all of the training sets. Finally, the specialised DP is further tested to insure the best overall performance.

Here, we show an example of DP specialisation for the mechanical response of Ti [101]. Figures 9(a) and (b) present a comparison of different empirical interatomic potentials with DFT on the -lines of Basal and Prism planes in HCP Ti. Illustrations of the Basal and Prism planes of HCP Ti are shown in the insets of figures 9(c) and (d). All of these empirical potentials yield inaccurate -line profiles and/or generally predict low stacking fault energy (the first minimum at 0.33 along the x-axis) on the Basal plane with respect to DFT. These empirical potentials are not specialised for dislocation properties. When general-purpose DP is obtained from DP-GEN scheme, it also shows similar systematic inadequacy as the other empirical potentials: low Basal and high Prism stacking fault energy, which prohibits the use of a generalised DP for investigating dislocation properties in HCP Ti. Specialisation is then performed by adding the structures in the red dashed box as well as the origin in figures 9(c), (d) and the structures along -lines in figure 9(e), all with a weight of 100. The DP -lines in figures 9(c)-(e) are overall in good agreement with DFT. In particular, the stable stacking fault energy order Basal $>$ Prism $>$ Pyramidal I narrow (the stable stacking fault energy on this plane will decrease significantly after full relaxation compared to that from -line) is different from previous empirical potentials in figures 9(a) and (b) and follows the trends observed in the DFT data. Based on this, the relative energies of the screw $⟨\mathbf{a}⟩$ dislocation on the Prism and Pyramidal I planes are shown in figure 9(f). The screw $⟨\mathbf{a}⟩$ dislocation is more stable on Pyramidal I plane than on Prism plane for DP, which agrees with previous DFT calculations and experimental measurements [113]. We refer readers to [101] for more details. As an additional example where the difference between DP and empirical potentials is particularly vivid; we consider the case of the screw $⟨c+a⟩$ dislocation core on the Pyramidal I plane of HCP Ti, as determined using DFT, DP [101], and the most widely used interatomic potential [109] in figure 10. The DP-determined core-structure shares the key features of the DFT-determined core structure; a long extended region and an asymmetric compact, localised core. The DP- and empirical interatomic potential-determined dislocation move very differently. This example of Ti DP demonstrates the importance of specialisation of the DP and the superior flexibility of DP over other empirical potentials.

Figure  9.  The generalised stacking fault energy (-lines) on the (a), (c) Basal, (b), (d) Prism, and (e) Pyramidal I plane of HCP Ti calculated using DFT, empirical potentials (Hennig [109], Ko [111], Dickel [112], and Mendelev (Ti3’ in [108])), and a specialised DP. The configurations in the dashed red box and at zero slip (origin) are included in the training sets. All configurations on the Pyramidal I narrow -line (e) are included in the training dataset. (f) Each point denotes the energy of a simulation cell containing a screw $⟨\mathbf{a}⟩$ dislocation core from a 600 K MD simulation after quenches to zero temperature (each point is one picosecond apart in the MD simulation). Two types of screw dislocation cores are observed: one delocalised onto a Prism plane (above the blue line) and one on Pyramidal I plane (below the blue line). Some of the cores are slightly distorted leading to small variations in the energy. See [101] for details.

DownLoad: Full-Size Img PowerPoint

Figure  10.  The screw $⟨c+a⟩$ dislocation core structure on the Pyramidal I plane of HCP Ti determined via DFT, DP [101], and the Hennig MEAM potential [109]. The red, blue, and white shaded atoms denote local HCP, BCC, and indeterminate atomic environments.

DownLoad: Full-Size Img PowerPoint

3.   DP applications in materials science

In the past three years, DPs have been applied in a number of systems in materials science including (1) elemental bulk systems, (2) multi-element bulk systems, (3) aqueous systems, (4) molecular systems and clusters, and (5) surfaces and low-dimensional systems. Table 2 shows a list of the material systems to which DPs have been applied (as of the writing of this paper). We choose several examples from each category to briefly discuss the corresponding DP application and how DP aids materials science research.

Table  2.  Example applications of DP in materials science.

System	Reference
Elemental bulk systems
Al	[83, 114-118]
Mg	[83]
Cu	[84]
Ti, W	[101, 102]
Ag, Au	[103, 119, 120]
Li	[121]
Be	[122]
Ga	[123]
Sb	[124]
C	[125]
Si	[126, 127]
P	[98]
Multi-element bulk systems
Al-Mg, Al-Cu-Mg	[83, 100, 128, 129]
Al-Cu, Al-Zn-Mg	[130, 131]
Al-Cu-Ni	[132]
Ag-Au	[103, 119]
Pd-Si, Nb5Si3, Zr77Rh23, Bi2Te3	[133-136]
Al90X10 (X = Tb, Cr, or Ce)	[137-140]
(Pd, Pt)${}_{\mathrm{x}}$(Ge, Sn, Pb)${}_{\mathrm{y}}$	[141]
P2Sn5	[142]
Silica, silicate	[143-146]
SiC	[104, 147]
B4C	[148]
Molten salt LiF, FLiBe, and chloride	[149-157]
Li or Na-based battery materials	[105, 158-162]
TiO2	[163]
-Ga2O3	[164]
Ferroelectrics HfO2	[165]
Ag2S	[166]
MoS2	[167]
SnSe	[168]
Zr${}_{1-x}$WxB2	[169]
(Hf0.2Zr0.2Ta0.2Nb0.2Ti0.2)X (X = C or B2)	[170, 171]
Aqueous systems
Water	[95, 172-183]
Zinc ion in water	[184]
Water-vapour interface	[185, 186]
Water-TiO2 interface	[187]
Ice	[188, 189]
Molecular systems and clusters
Organic molecules	[99, 190-195]
Metal and alloy clusters	[119, 196]
Surfaces and low-dimensional systems
Metal and alloy surfaces	[103, 119, 129]
Graphane	[125, 197]
Monolayer In2Se3	[198]
2D Co-Fe-B	[199]

 | Show Table

DownLoad: CSV

3.1   Elemental bulk systems

To date, DPs have been applied to a wide-range of pure systems, including Al, Mg, Cu, Ti, W, Ga, C, Si, as shown in table 2. Al was the first metal system to which DP was applied and a general-purpose DP developed [83]. It accurately reproduces the lattice parameter, elastic constants, vacancy and interstitial formation energies, surface energies, twin and stacking fault energies, melting point, enthalpy of fusion, diffusion coefficient. The general-purpose DP for other metal elements is also accurate for the same properties. For other properties, not included in the training datasets, DP is in better agreement with DFT than MEAM [10] for phonon dispersion relations, equations of state, and the liquid state radial distribution function [83]. Based on this, Wang et al [114] smoothly interpolated the Ziegler-Biersack-Littmark (ZBL) screened nuclear repulsion potential with a DP to obtained a DP-ZBL model for irradiation damage simulations, surpassing the widely adopted ZBL MEAM [200] or EAM [201] potentials. Later, a DP was developed for warm dense Al to simulate ion dynamics near the hydrodynamic limit [115], structural and dynamic properties [116], and electronic and ionic thermal conductivities [117]. A DP was also developed for high temperature and high pressure liquid Al to calculate shear velocity [118].

Dislocation properties play important roles in the plastic response of most structural materials, including Ti and W. The specialised Ti DP accurately depicts the -lines on different planes (figures 9(c)-(e)) and the screw $⟨\mathbf{a}⟩$ dislocation core energy ordering between the Prism and Pyramidal I planes (figure 9(f)). In addition, the screw $⟨\mathbf{a}⟩$ dislocation core structures on Prism and Pyramidal I planes are also in surprisingly good agreement with DFT results [113]. Because the dislocation core structures could not be explicitly included into DFT training datasets, the example of Ti DP shows that surrogate properties can be used to optimise DPs for dislocation properties of complex HCP systems. While the dislocation core structures of BCC W are not as complex as in HCP Ti (because it only adopts a compact core), the Peierls barrier for BCC W is difficult to reproduce with other empirical potentials (cf DFT results [102]). The DP-SE2 model with only a two-body embedding descriptor yields a very poor prediction on the Peierls barrier, but the DP-HYB that hybridises descriptors with two-body and three-body embeddings reproduces this property very accurately [102]. The Ti and W examples demonstrate that DPs can accurately describe dislocation core structures and Peierls barrier of BCC and HCP metals.

The DP approach has also been applied to Ag and Au (widely used in catalytic applications). Andolina et al [119] and Wang et al [103] developed DP for Ag and Au that is accurate at lattice parameters, elastic constants, surface formation energies, interstitial and vacancy formation energies, etc. Furthermore, Andolina et al [119] got accurate adsorption energy and diffusion barriers for adatoms on {100}, {110}, and {111} compared to DFT results. Wang et al [103] presented a comprehensive study of the Au {111} surface reconstruction using a DP that yields excellent agreement with DFT results. From another perspective, Chen et al [120] used DP MD to illustrate the dynamics compression process of Au. The developed DP could accurately reproduce the experimentally determined phase boundaries and the short-to-medium range orders are proposed to reduce the Gibbs free energies of the shocked structures. The examples above for Ag and Au validate the applications of DP in both catalytic and shock compression areas.

In addition, DP has also been applied in many other elemental bulk systems and here we only list some of the examples in table 2. For Li, which is an important element for battery, Jiao et al [121] developed DP to reveal self-healing mechanisms in a large Li-metal system. For Ga, Niu et al [123] used DP to construct the phase diagram of liquid Ga, -Ga, -Ga, and Ga-II, in good agreement with experimental results. In addition, the local structure of liquid Ga and the nucleation process into -Ga, -Ga were also studied [123]. For carbon, a DP was developed by Wang et al [125] to simulate the structural properties of 12 different bulk and low-dimensional carbon structures. For Si, which has both covalent and metallic bonding behaviour, the first DP was trained based on datasets generated by classical metadynamics simulations [126]. The DP was then applied to study the crystallisation and the free energy surface between liquid and solid. Many thermodynamics properties near the critical point were found to be close to experimental data. Li et al [127] trained a DP on DFT data of silicon in the crystalline, liquid, and amorphous phases and thermal conductivity was accurately reproduced. DPs were used to study a liquid-liquid phase transition in P, where DP established the main features of the liquid phase diagram [98]. In particular, DP phase diagrams of Ga and liquid P clearly indicate how DPs accurately describe the PES of different phases and transitions between them.

The applications in elemental bulk systems are among the first applications of DPs in materials science. Their successes for different crystal structures and various classes of properties (mechanical, catalytic, irradiation properties, phase transformation, thermal conductivity, etc) encouraged their extension to multi-element bulk systems and increasingly complex phenomena.

3.2   Multi-element bulk systems

Al-Mg was the first alloy system for which an accurate DP was developed [83]. This DP was used to describe the 28 crystalline Al-Mg alloys structures in the Materials Project (MP) database [202]; include accurate prediction of formation energies, equilibrium volumes, elastic constants, vacancy and interstitial formation energies, and unrelaxed surface energies. Wang [128] applied this Al-Mg DP and validated its reliability for crystal structure prediction by using DP+CALYPSO. This work interfaces the DeePMD-kit package and crystal structure prediction software (e.g. CALYPSO [203], USPEX [80], and Pychemia (web link 13 in section S2 of the SM)). Andolina et al [129] developed a Al-Mg DP based on the original DP to investigate anisotropic surface segregation. Based on the Al-Mg DP, a Al-Cu-Mg ternary DP was developed for the entire compositional space [100]. 2.73 billion alloy configurations were explored in the DP-GEN process. The resulting DP yields more accurate results for energetic, mechanical, and defect properties of 58 crystalline structures as compared with MEAM potentials [10]. The multi-component DP approach can readily be applied to high-entropy alloys for which adequate empirical potentials are difficult to obtain.

Bourgeois et al [130] built an Al-Cu DP to simulate the aggregation of vacancies around embedded ${\theta }^{\mathrm{\prime }}$ precipitates and investigated the nucleation of this strengthening phase onto a template structure. An Al-Mg-Zn DP was developed [131] and applied to confirm the co-segregation of Mg and Zn atoms at a precipitate and matrix interface. DPs have also been found to be powerful and promising for the prediction of the structure and dynamics of metallic liquids, glasses, and quasi-crystal [133-140]. The DP developed for Pd-Si accurately represented the structure of liquid and crystal structures, melting points, and glass-forming ability at compositions near Pd₃Si and Pd₉Si₂ (more accurately than existing EAM potential) [133]. Accurate liquid structure and dynamic properties were also obtained with DPs for Nb₅Si₃ [134], Zr₇₇Rh₂₃ [135], and Bi₂Te₃ systems [136].

Tang et al [137-140] performed DP MD simulations of a series Al-based alloys; we focus now on Al-Cr quasicrystals [139] as an example and a demonstration of how DPs can be used together with experimental studies. Dendritic growth of metastable quasicrystals were observed in the Al₁₃Cr₂ approximant phase (formed from Al₉₀Cr₁₀ thin film) by pulsed laser deposition [139] which is structurally similar to quasicrystal of the Al₁₃Cr₂ matrix. The Al-Cr DP was used to simulate the quenching of the Al₉₀Cr₁₀ alloy from 2200 to 700 K at 10¹¹ K s^-1. There are three types of 13-atom icosahedra in the approximant Al₁₃Cr₂ phase and one icosahedral Al-Cr quasicrystal motif. These 4 icosahedral motifs are similar despite slightly different Cr-Al bond lengths. All 4 types of 13-atom icosahedral motifs were Cr-centred. The icosahedral motif appears in both the quasicrystal and approximant structures, which results in the survival of the 13-atom icosahedron after laser irradiation. This was observed in both simulation and experiment; the success of the DP was attributed to the excellent reproduction of the liquid structure in the DP-based simulations.

A DP was developed to simulate liquid and glassy silica which proved to have satisfactory accuracy based upon a relatively small training dataset [143]. Other DPs were developed to calculate transport properties of silicate in the mantle [144-146]. DPs were also employed in large-scale calculations of thermodynamic, transport, and structural properties in different molten salts [149-157].

Another exciting class of DP applications was for Li and/or Na-based battery materials [105, 158-162]; here we focus on the example of Li₁₀GeP₂S₁₂-type superionic conductors [105]. The DP-GEN scheme was used to generate DPs for three solid-state electrolyte materials (Li₁₀GeP₂S₁₂, Li₁₀SiP₂S₁₂, and Li₁₀SnP₂S₁₂) and applied to diffusion over a wide temperature range with 1000 atoms. The predicted diffusion coefficients slightly overestimated the experimental values but were within the experimental uncertainty. These DP-based simulations provided a starting point for large size scale and long time scale MD investigations of solid-state electrolyte materials.

Additional DPs were developed for a wide-range of other multi-element bulk systems, including metal oxide [163-165], metal sulfide [166, 167], thermoelectric SnSe materials [168], metal borides [169, 171], and metal carbide [170] systems.

3.3   Aqueous systems

Since the original DP paper [48], water has been widely studied with DPs. Ko et al [172] applied DPs to perform extensive sampling of thermal and nuclear quantum fluctuations on an accurate PES. In particular, a DP was used to investigate isotopic effects on structural properties of liquid water (H₂O and D₂O). To understand the experimentally observed isotope effect in the x-ray absorption spectra between liquid H₂O and D₂O, DP-based, path-integral MD simulations were performed [174]. A more comprehensive study by the same group [177] examined isotope effects on radial distribution functions, O-O-O triplet angular distributions, density and found that the DP-based simulations were consistent with the experimental observations. Recently, Calio et al [183] performed DP-based simulations to interpret experimental observations on the hydrated excess proton in water.

Sommers et al [173] trained a DP to predict the polarisability of liquid water with ab initio accuracy in order to calculate the Raman spectra in long time scale. Gartner et al [175] trained a DP for water to examine the questions of the existence of a liquid-liquid transition in water. Andreani et al [176] combined neutron scattering experiments and DP MD to investigate hydrogen dynamics in supercritical water. The vibrational densities of states observed in DP MD showed coupling between intramolecular vibrational and intermolecular librational and rotational motions. Piaggi et al [178] used a DP to investigate ice nucleation in water, hexagonal ice, and cubic ice and obtained quantitative agreement between DP and experiment (better than the state-of-the-art semi-empirical potentials). A more complete description of the phase equilibrium between different phases of H₂O was achieved through constructing the DP phase diagram from low temperature and pressure to 2400 K and 50 GPa [95]. Although both the DP and empirical TIP4P/2005 potential [204] work well at low and intermediate pressures, DP outperforms TIP4P/2005 at higher pressure compared to experimental results, especially at the phase boundaries between ice VIII, VII, and VI (figure 11). Unlike in experiments, the TIP4P/2005 predicts a first-order transition from ice VII to a plastic phase (no such phase transformation occurs using DP). This is a milestone for DP considering the importance of water, the vast range of temperatures and pressures, and the high accuracy required for free energy representation.

Figure  11.  Phase diagram of water. (a₁) DP model (red solid lines) and experiment (gray solid lines) for $T<$ 420 K. (a₂) Phase diagram at high T and P. (b) Phase diagram of TIP4P/2005 [204] water. Reprinted figure with permission from [95], Copyright (2021) by the American Physical Society.

DownLoad: Full-Size Img PowerPoint

Tisi et al [179] calculated the thermal conductivity of water using both DFT (SCAN) and DP. Interestingly, both approaches yield the same conductivity which was 50% higher than the experimental value. Zhang et al [180] applied DP to help improve the exchange functional in DFT from SCAN to SCAN0 for water. Similarly, Torres et al [181] evaluated the errors of DFT-based simulations on structural and dynamical properties due to time- and size-scale limitations by using DP MD. In these two examples, the high efficiency and accuracy of DP provided rapid screening of different properties to feedback into DFT exchange-correlation functional optimisation. Shi et al [182] extended DP to produce accurate molecular multipole moments in the bulk and near interfaces consistent with AIMD simulations. These moments were used to compute the electrostatic potential at the centre of a molecular-sized hydrophobic cavity in water.

The DP approach has also been used in a wide range of aqueous systems. Xu et al [184] developed a DP to perform MD study of zinc ions in liquid water. The experimentally observed zinc-water radial distribution function, as well as the x-ray absorption near edge structure spectrum, was well-reproduced by the DP MD simulation. Recently, the Limmer group [185, 186] applied DP to study liquid-vapour interfaces. They found that the DP yielded accurate interfacial properties by incorporating explicit models of the slowly varying long-ranged interactions and training neural networks only on the short-ranged components [185]. In addition, they trained a DP for solvated N₂O₅ and bulk ambient water and applied DP MD and importance sampling to study the uptake of N₂O₅ into an aqueous aerosol [186]. In contrast to the previous understanding that the uptake process occurs within the bulk of an aerosol, interfacial processes dominate the uptake process due to facile hydrolysis at the liquid-vapour interface and competitive re-evaporation. This work not only brings new insights to a long-standing questions, but also extends the application of DPs to the liquid-vapour interface. Other examples of DP applications in aqueous systems include TiO₂-water interfaces [187], the ice Ih/XI transition [188], and dynamical states of high-pressure ice VII [189].

3.4   Other systems

Here we briefly list several applications of DPs to other classes of systems including molecular systems, clusters, surfaces, and low-dimensional systems. Jiang et al [190] developed DPs for sulfuric acid-sulfuric acid, dimethylamine-dimethylamine, and sulfuric aciddimethylamine organic molecular systems to investigate the atmospheric aerosol nucleation process. Zeng et al [191] trained a DP based on a dataset for the pyrolysis of n-dodecane and performed a reactive DP MD simulation to reveal the detailed pyrolysis mechanism, in good agreement with experiment. Chen et al [192] used a DP to accurately represent the ground- and excited-state PES of CN₂NH. This DP accurately reproduced excited-state topological structures, photo-isomerisation paths, and conical band structure intersections. Yang et al [99] first generated training datasets through active learning with enhanced sampling and then developed a DP to study the urea decomposition process in water. Wang et al [194] presented a data-driven coarse-grained simulation of polymers in solution and validated the accuracy of this method with DPs to construct a coarse-grained potential. Pan et al [195] extended the DP-approach to incorporate external electrostatic potentials in a molecular system; the resultant DP was accurate for energies and forces of representative configurations along the Menshutkin and chorismate mutase reactions pathways.

A study of metal and alloy clusters and surfaces demonstrated a conflict between Al bulk and cluster energies [196]. This indicates that the compromise between properties are, on occasion, necessary and suggest that DPs should be developed for target properties (this is the specialisation discussed above). Andolina et al applied a DP to study the nucleation and growth of seeded core-shell Ag and Au clusters [119] and predict the anisotropic surface segregation for Al-Mg alloys [129]. Wang et al [103] successfully applied DPs in Ag-Au, Au {111} surface reconstruction and segregation of Au on the Ag-Au nanoalloy surfaces.

Achar et al [197] proposed a DP for graphane and showed that it outperforms empirical interatomic potentials for phonon density of states, thermodynamic properties, velocity autocorrelation function, and stress-strain curve up to the yield point. Wu et al [198] developed a DP for In₂Se₃ monolayers and used it to predict a range of thermodynamic properties of In₂Se₃ polymorphs and lattice dynamics of ferroelectric In₂Se₃ with ab initio accuracy. Chen et al [199] applied a DP to simulate the synthesis of amorphous CoFeB during a rapid cooling process. The applications of DPs in low-dimensional systems are in their early stage, but existing evidence suggests that the DP method is promising for the simulation of low-dimensional materials.

4.   Accuracy and efficiency of DPs in practice

The summary of recent applications of the DP method in different systems suggests the wide-applicability and high accuracy of DPs. Nonetheless, it is appropriate to return to the competition between DP accuracy vs the efficiency of DP-based simulations. On the accuracy issue, we reflect on what we have learned from the applications performed to-date and our experience, to address: (1) How good is DP? (2) What have we learned? (3) When can we rely on DPs and when can we not? On the efficiency issue, we summarise how fast DPs are by comparison with other approaches from pair potentials to DFT.

4.1   Accuracy

4.1.1   How good are DPs?.

Broadly speaking, DPs are more accurate than other types of empirical interatomic potentials; this was the main reason behind the development of the DP method and its application to a wide range of systems. The improvement in accuracy often leads to qualitatively’ new results. We return to two examples where DPs are qualitatively’ better. The first example is metal Ti [101]. Previous experiments and DFT calculations [113] confirmed that the screw $⟨\mathbf{a}⟩$ dislocation is more stable on the Pyramidal I than on the Prism plane. DFT calculations, using pseudopotentials with different numbers of valence electrons [205], showed that the dislocation core energy of the screw $⟨\mathbf{a}⟩$ dislocation on the Pyramidal I plane is 18.4 meV b^-1 lower than that on the Prism plane. This energy ordering is a prerequisite for the phenomena of screw $⟨\mathbf{a}⟩$ dislocation locking’ and unlocking’ observed in experiment at 150 K [113]. While the screw $⟨\mathbf{a}⟩$ dislocation is more stable on the Pyramidal I plane, it has a high energy barrier for glide; hence the screw is locked’ on this plane. However, the screw $⟨\mathbf{a}⟩$ has a much lower glide barrier on the higher energy Prism plane; hence on this plane, the screw unlocks’ and glides easily. Unfortunately, all extant empirical Ti potentials predict the incorrect energy ordering; i.e. the $⟨\mathbf{a}⟩$ screw is more stable on the Prism plane (including the generally excellent MEAM empirical potential [109]), in qualitative’ disagreement with experiment. The powerful representability and flexibility of the DP method allows the Ti DP to reproduce this important feature of defects in Ti and enables new insights through large-scale MD simulations.

Another example of qualitative’ improvements made possible by a DP is the calculation of the phase diagram of water [95]. TIP4P/2005 is one of the most accurate empirical water models available today for phase diagram prediction [95]. At high temperature and pressure, TIP4P/2005 predicts a first-order transition from ice VII to a plastic phase, in which the BCC oxygen sublattice coexists with freely rotating molecules. This prediction has not been experimentally confirmed. The DP for water predicts a super-ionic ice VII in that region, in agreement with recent experimental observations [206]. DP also better reproduces the phase boundaries between ice VI, VII, and VIII, at high pressures better than TIP4P/2005. These are two of many examples of where the improved quantitative predictions of DP enable qualitatively’ correct phenomena not accessible through other empirical interatomic potentials.

The accuracy of DP is routinely compared with DFT results especially when benchmarks from widely accepted empirical interatomic potentials do not exist. In most examples (and our own experience), the agreement between energy and force (RMSE) values obtained using a DP and DFT is commonly smaller than 10 meVatom^-1 and 100 meV Å^-1, respectively. In some systems, such as the water phase diagram, the error is typically 1 meVmolecule^-1. This accuracy should meet the requirements for most applications. In particular, almost half of the examples listed in table 2 focus on non-crystalline systems (liquids, amorphous systems; including liquid metallic alloys, and metallic glasses). DP-based predictions are especially accurate for liquid properties (DFT references) on liquid structures, diffusion coefficients, thermal conductivities.

While the flexibility of DPs is desirable for describing complex potential energy surfaces, the ultimate accuracy of DPs are often limited by the accuracy of the DFT training set. Since DFT is itself an approximation, perfect agreement with nature is not expected. Improvements are possible with improved exchange and correlation DFT functionals; climbing the Jacob’s ladder from LDA to GGA to meta-GGA to hybrid-GGA to fully non-local approaches [75]. In practice, the ultimate limit of DP accuracy can often be associated with the DP training sets, the accuracies of which are often associated with the choice of the DFT functionals. This choice is often dictated by the associated computational cost which rises rapidly on climbing the rungs of the Jacob’s ladder.

4.1.2   What have we learned?.

The two examples, above, demonstrate that DPs can lead to very high accuracy results; comparable to the underlying DFT approaches and higher than the vast majority of empirical interatomic potentials. The examples, above, also show us that DPs can be employed to add new understanding and insight in situations which were previously inaccessible to other computational approaches and experimental observation. The DP approach is also useful in the development of new ab initio methods and pseudopotentials [180, 181]. Finally, the examples showed how DPs can be readily specialised to describe phenomena for which general purpose DPs do not suffice.

At early stages of DP development and for simple applications (e.g. simulation of liquid state structure) where the demands on accuracy are not too high, sampling efficiency is not so critical. For such simple applications, an initial DP based on a small set of AIMD trajectory training sets, often leads to efficient and accurate results. However, sampling efficiency is critical when a general-purpose DP is required or target properties are subtle (e.g. where energy differences between phases are very small or some defect properties). AIMD may not adequately/efficiently sample atomic configurations that represent those associated with properties of interest. This may be addressed by starting with the general purposed DP developed through DP-GEN and tweaking it through identification of property-appropriate and incorporation in the specialisation step, as described above. The specialisation process varies with systems and properties of interest and may be viewed as the art of tweaking DPs based upon physical understanding.

4.1.3   When can we rely on a DP?.

After performing our normal suite of property testing on a DP (e.g. see [101]), experience shows that such DPs yield reliable results in atomistic simulationsespecially compared with DFT calculations that are of insufficient spatial scales or empirical interatomic potentials that are of insufficient accuracy. Critical issues for all ML potentials are representability and transferability. Representability implies the ability of the functional form to accurately reproduce the target properties. Transferability is the ability of a potential to describe the properties which were not included in the training process.

The DP approach usually performs well from a representability perspective; DPs are usually able to provide fits that adequately represent all of the training datasets. In some cases DP failed to distinguish similar configurations (e.g. the configurations along the transition path of the screw dislocation in BCC W); the representability of the DP can be improved by using more expressive descriptors, such as a three-body embedding descriptor [102]. On the other hand, transferability can be non-trivial and subtle for DPs. Transferability can be classified as in-distribution and out-of-distribution transferability. In-distribution or out-of-distribution transferability is the ability of model to interpolate within or extrapolate out of the sampled configuration distribution, respectively. The in-distribution transferability of DPs trained with the DP-GEN scheme is generally quite good, providing reliable and accurate predictions of configurations similar to those in the sampled distribution. However, ML potentials will fail in the out-of-distribution transferability where the explored configuration is far’ from those sampled configurations used in training [34]. A simple example of an out-of-distribution issue is a DP trained using only liquid datasets; such a DP normally shows poor transferability with respect to crystal datasets because there is little overlap between the liquid and solid configuration distributions. For defect properties, we do not know a priori whether defect configurations represent an in- or out-of the distribution with respect to those sampled by DP-GEN. Model deviation (equation (31)) serves as a good indicator of the transferability of a DP (without the need for additional DFT calculations). In the case of transferability failure, the DP can be specialised by adding configurations to the training set that more closely represent the configuration of interest; this is converting out-of-distribution transferability failure to in-distribution transferability agreement.

In the cases when a user is not confident whether a DP is transferable, the DP can be used in conjunction with model deviation after validating against relevant DFT or experimental benchmarks. DPs are, like other empirical potentials subject to the adages that interatomic potentials will work only before they fail’ or will work until they do not’. From this perspective, the replacement of DFT by ML potentials including DP are not completely reliable; DFT will remain the method of choice where very high accuracy property prediction is necessary. The combination of DFT and DP provides a practical strategy for the needs of the materials science community, providing the implicit trade-offs between accuracy and computational efficiency.

4.2   Efficiency

Using the previous DP application examples, we see that (1) on CPUs, compressed DPs are faster than DFT by a factor of over 10⁶ and slower than empirical interatomic potentials such as EAM (MEAM) by 100 (10) times; (2) on GPUs, DP compression model is slower than potentials like EAM by, of order 10 times. Of course, the actual efficiency is application dependent (especially for comparison with DFT). Compressed DPs can be faster than the original DPs by a factor of over 10 and consume an order of magnitude less memory. Additional optimisation of the neural networks at the heart of DPs is possible by optimisation of different operators on the computational graph and through hardware changes [107]. Both DPs and empirical potentials show linear scaling with the number of atoms on both CPU and GPU machines. This linear scaling is the enabler of large-scale atomistic simulations. The lower speed of DPs, compared with empirical potentials, is reasonable considering the vast number of parameters in DPs (often at the order of 10⁵). Because DPs are and will continue to be slower than empirical potentials, simpler, empirical potentials will continue to play an important role in materials science. The competition’ between MLs and simpler, empirical potentials drives the continuous improvement of each. The ML potential community continues to focus on improving the potential efficiency (computational speed), while empirical potential development continues to develop new formalism increasing accuracy. The concept of ML potentials guided by the physics in the empirical potentials’ is also an exciting area; e.g. see the recent work of Mishin et al [31, 207].

5.   Conclusions and outlook

With increasing need for atomistic simulations with higher accuracy, larger length scales, longer time scales and computational efficiency, ML-based interatomic potentials are rapidly gaining acceptance in the broad materials science community. This is especially true in areas where the phenomena of interest are subtle and those in which the material system is complex. In this review, we examined the DP approach (for ML potentials); summarising the basic theory, how to develop DPs and apply DP-development software and database, how to make DPs more efficient in applications, how to specialise DPs for subtle application, reviewed several DP applications, and discussed DP accuracy and efficiency. After several years of evolution, the DP method is now relatively mature, yet continuing to improve in both accuracy and efficiency within an open-source community framework. We envision the DP method to continue developing in the coming years and the continued expansion of the data base of useable DPs.

Continued development will likely proceed along several avenues. The first is the development of new and more intelligent descriptors for better predictability. We see from the W example [102] that the Peierls barrier (the barrier for dislocation glide) can only be accurately reproduced by expanding the DP descriptor to include three-body embedding. We suspect that such examples will continue to arise as DPs are expanded to include a broader set of applications in different materials. Another issue is related to magnetism; how can magnetic moment degrees of freedom be incorporated into DPs? Empirical EAM and MEAM potentials deal with this issue through a set of assumptions and approximations. More intelligent descriptors also improve the ease of DP training. For example, questions arise for the current hybrid descriptors, with two-body and three-body embeddings about how much to weight three-body embeddings in the hybrid descriptors. Current strategies in this area tend to be based largely on empirical experience; hence, there are opportunities to transform this into a machine-driven process.

An important second area of development will be improvement of the automation of DP training and specialisation. In DeePMD-kit and DP-GEN software, different settings do, on occasion, influence the performance of trained DPs. Although we presented our experience on choosing these settings in section 2.3, ideally this experience should be replaced in future generations of the DeePMD-kit and DP-GEN software. The automatic selection of the trust levels in DP-GEN is already in the testing stage. The specialisation step should also be more automated to reduce user intervention in determining: (1) what types of specialisation datasets are needed? (2) how many specialisation datasets are needed to combine with DP-GEN datasets? (3) when to include specialisation datasets and how to modify DeePMD-kit and DP-GEN settings? More automated training and specialisation schemes would accelerate the development of new DPs for more systems and applications.

A third area of DP future development is further optimisation of the computational speed of DP. Currently, many empirical potentials are faster than DP (at least a factor of 10) which leads potential users to prefer empirical potentials when accuracy demands are not high. In principle, DPs cannot be faster than empirical potentials considering the vast number of parameters involved, but decreasing the computation efficiency loss in using a DP can change the speedaccuracy tradeoff and enable more effective material simulations. As seen above, the greatly improved accuracy of DPs over most empirical potentials opens the door to applications where empirical potentials are simply qualitatively’ incorrect.

A fourth area is the combination of DP and transformer [208]. Transformers were introduced in 2017 and are now commonly used in the fields of natural language processing (NLP) and computer vision. Large pre-trained models based on the transformer architecture (e.g. BERT [209] and GPT-3 [210]) are changing the research paradigm of NLP. Traditionally, models for a specific task are trained from scratch. The Transformer approach is based upon fine-tuning pre-trained models such that new records are being established for specific tasks [211]. For DP, pre-trained models are trained using a large database based upon different elements in the periodic table; learning the embedding mapping from the element type to the abstract space. For example, when we have training datasets for binary systems and want to train a DP model for a ternary system, the traditional step would be using DPGEN to generate datasets for ternary system from scratch. However, with the introduction of pre-trained binary system DP models, the number of additional datasets required to train a ternary system DP model is greatly decreased. From this perspective, pre-trained models can accelerate the training process and decrease the computational time required for generating the requisite training datasets.

Finally, the open-source DP Library (a database for DP including training datasets, training schemes, DP, and testing results) should be greatly expanded to include most of the periodic table, figure 6 and alloys. This will be an on-going effort requiring contributions from the entire user group. Another task is to make the DP Library easier to use. Both of these will be enabled by improvements in the openness of scientific computing community (e.g. appropriate acknowledgment of contributions).

The development of the DP approach (and related ML potential approaches) represents an important milestone for the field of atomistic simulations of materials that rests on advances in ML technology and descriptors of atomic environments. DPs routinely provide high (near DFT) accuracy with reasonable computational efficiency, as compared with empirical potentials. The accuracy and efficiency of DPs open the door to qualitatively’ new applications of atomistic simulation.