Deep learning inter-atomic potential model for accurate irradiation damage simulations

TL;DR

This paper introduces a hybrid deep learning potential model that accurately simulates irradiation damage in materials, capturing physics at atomic scales and outperforming traditional potentials in aluminum simulations.

Contribution

A novel hybrid scheme combining ZBL potential with a deep learning model for improved irradiation damage simulations in materials.

Findings

Better prediction of defect formation energies

More accurate simulation of collision cascades

Enhanced modeling of residual point defects

Abstract

We propose a hybrid scheme that interpolates smoothly the Ziegler-Biersack-Littmark (ZBL) screened nuclear repulsion potential with a newly developed deep learning potential energy model. The resulting DP-ZBL model can not only provide overall good performance on the predictions of near-equilibrium material properties but also capture the right physics when atoms are extremely close to each other, an event that frequently happens in computational simulations of irradiation damage events. We applied this scheme to the simulation of the irradiation damage processes in the face-centered-cubic aluminium system, and found better descriptions in terms of the defect formation energy, evolution of collision cascades, displacement threshold energy, and residual point defects, than the widely-adopted ZBL modified embedded atom method potentials and its variants. Our work provides a reliable and feasible scheme to accurately simulate the irradiation damage processes and opens up new opportunities to solve the predicament of lacking accurate potentials for enormous newly-discovered materials in the irradiation effect field.

Tables

Table 1. Table 1: Equilibrium properties of Al: atomization energy E am subscript 𝐸 am E_{\rm am} , equilibrium lattice constant a 0 subscript 𝑎 0 a_{\rm 0} , vacancy formation energy E vf subscript 𝐸 vf E_{\rm vf} , interstitial formation energy E if subscript 𝐸 if E_{\rm if} for octahedral interstitial (oh) and tetrahedral interstitial (th), independent elastic constant C 11 subscript 𝐶 11 C_{\rm 11} , C 12 subscript 𝐶 12 C_{\rm 12} , and C 44 subscript 𝐶 44 C_{\rm 44} , Bulk modulus B V subscript 𝐵 𝑉 B_{V} (Voigt), shear modulus G V subscript 𝐺 V G_{\rm V} (Voigt),stacking fault energy γ sf subscript 𝛾 sf \gamma_{\rm sf} , twin stacking fault energy γ tsf subscript 𝛾 tsf \gamma_{\rm tsf} , melt point T m subscript 𝑇 m T_{\rm m} , enthalpy of fusion Δ ​ H f Δ subscript 𝐻 𝑓 \Delta H_{f} and diffusion coefficient D 𝐷 D at T 𝑇 T = 1000 K.

Al111The results above, unless specified with a reference, are computed by the authors.	EXP.	DFT	DP-ZBL	DPZhang et al. (2019)	MEAM-ZBL	EAM-ZBL222The interstitial and stacking fault configurations were unstable upon relaxation with the EAM-ZBL potential, so their formation energies are not reported here.
$E_{\mathrm{am}}$ [eV/atom]	-3.49J. D. Cox (1990)	-3.75	-3.74	-3.65	-3.36	-3.39
$a_0$ [Å]	4.04Cooper (1962)	4.04	4.04	4.04	4.05	4.01
$E_{\mathrm{vf}}$ [eV]	0.66Triftshäuser (1975); Fluss et al. (1978)	0.67Hood, Kent, and Reboredo (2012)	0.73	0.79	0.67	1.14
$E_{\mathrm{if}}\left(oh\right)$ [eV]	-	2.91Hood, Kent, and Reboredo (2012)	2.57	2.45	3.12	-
$E_{\mathrm{if}}\left(th\right)$[eV]	-	3.23Hood, Kent, and Reboredo (2012)	3.23	3.12	3.83	-
$C_{11}$ [GPa]	114.3Kamm and Alers (1964)	111.2	112.8	120.9	113.5	106.9
$C_{12}$ [GPa]	61.9Kamm and Alers (1964)	61.4	57.6	59.6	61.6	81.5
$C_{44}$ [GPa]	31.6Kamm and Alers (1964)	36.8	41.2	40.4	45.4	44.2
$B_V$ [GPa]	79.4Kamm and Alers (1964)	78.0	76.0	80.1	78.9	90.0
$G_V$ [GPa]	29.4Kamm and Alers (1964)	32.1	35.8	36.5	37.6	31.6
${\gamma }_{\mathrm{sf}}\left[J/m^2\right]$	0.11-0.21Kannan and Thomas (1966); Dobson, Goodhew, and Smallman (1967); Tartour and Washburn (1968); Mills and Stadelmann (1989)	0.142Zhao et al. (2016)	0.070	0.132	0.184	-
${\gamma }_{\mathrm{tsf}}\left[J/m^2\right]$	-	0.135Zhao et al. (2016)	0.075	0.130	0.184	-
$T_m$ [K]	935Ross, Yang, and Boehler (2004)	950($\pm 50$)Bouchet et al. (2009)	885	918	950	1050
$\Delta H_f$ [KJ/mol]	10.7($\pm$ 0.2)McDonald (1967)	-	9.3	10.2	11.5	8.8
$D\left[{10}^{-9}{m}^2/s\right]$	7.2-7.9Meyer (2015)	-	6.8	7.1	4.9	6.8

Table 2. Table 2: The average displacement threshold energy ( E d subscript 𝐸 d E_{\rm d} ) of fcc Al.

fcc Al	RecommendedNorgett, Robinson, and Torrens (1975b)	DP-ZBL	MEAM-ZBL	EAM-ZBL
$E_d\left(\mathrm{eV}\right)$	25.0	26.54	22.67	16.73

Equations

$$E^{\text{ZBL-DP}}=\sum _i{E}_i^{\text{ZBL-DP}}$$

$$E_i^{\text{ZBL-DP}}=E_{s(i)}^{\text{ZBL-DP}}({\mathbf{R}}_i,\{{\mathbf{R}}_j|j\in {\mathcal{N}}_{R_c}(i)\})$$

$$E_i^{\text{ZBL-DP}}=w_i{E}_i^{\text{ZBL}}+(1-w_i)E_i^{\text{DP}},$$

$$w_i=\{\begin{array}{rlrl}& 1& & {\sigma }_i<R_a,\\ & -6u_i^5+15u_i^4-10u_i^3+1& & R_a\le {\sigma }_i<R_b,\\ & 0& & {\sigma }_i\ge R_b,\end{array}$$

$$u_i=\frac{{\sigma }_i-R_a}{R_b-R_a},$$

$${\sigma }_i=\frac{{\sum }_{j\in {\mathcal{N}}_{R_c}(i)}{R}_{ij}{e}^{-R_{ij}/\alpha }}{{\sum }_{j\in {\mathcal{N}}_{R_c}(i)}{e}^{-R_{ij}/\alpha }},$$

Full text

††preprint: AIP/123-QED

With the rapid growth of computing science and computer performance, computational simulations, including molecular dynamics (MD) Jorgensen, Maxwell, and Tirado-Rives (1996); Dumpala et al. (2015); Hooda et al. (2016) and density functional theory (DFT) method Car and Parrinello (1985); Kohn and Sham (1965), are becoming increasingly important to evaluate the properties of materials. However, the accuracy of empirically constructed atomic potential models for MD simulations are often in question, while the quantum mechanics approaches, such as DFT, are limited by the time and size scale of the simulated systems. Therefore, a solution that combines the advantages of both methods is needed.

Recently, machine learning (ML) methods have been used to solve this dilemma Zhang et al. (2018); Yao et al. (2017); Han et al. (2018); Chmiela et al. (2017); Behler and Parrinello (2007); Bartók et al. (2010); Schütt et al. (2017). Several studies have demonstrated that the ML-based potential energy surface can reach the accuracy of DFT, with the cost comparable to classical empirical potentials Bartók et al. (2010); Schütt et al. (2017); Bartók, Kondor, and Csányi (2013); Zhang, Wang, and E (2018a). Nevertheless, challenges have remained for ML-based methods to describe very short-distance interactions, e.g. those in the irradiation damage processes. In these processes, the distance between atoms can be very short, and the interactions can hardly be treated as quasi-static, wherein conventional DFT approaches may fail, so only the Ziegler-Biersack-Littmark (ZBL) screened nuclear repulsion potential Ziegler and Biersack (1985) has been validated for a good description of the corresponding interactions. In other words, in this case, energies and forces from DFT calculations may no longer be accurate training data for ML-based potentials. Moreover, the magnitude of energies and forces is much larger than that in systems near equilibrium, which may pose additional difficulties for the training of ML-based potentials. Therefore, it is necessary to develop a new scheme that is applicable for irradiation damage simulations while still remains the accurate predictions of material properties for both near-equilibrium state and short-distance interaction.

To solve this problem, we interpolate the ZBL potential into a deep learning model, so that short-distance collisions between atoms can be accurately described. In our previous studies, we have developed the Deep Potential (DP) scheme, an end-to-end symmetry preserving machine learning-based inter-atomic potential energy model, which can efficiently represent the properties of a wide variety of systems with the accuracy of ab-initio quantum mechanics models Han et al. (2018); Zhang et al. (2018). This ZBL-modified deep learning scheme (DP-ZBL), which can be seen as an improved and specialised version of the original DP model, makes it possible to accurately simulate the irradiation bombardment damage for materials. Here in this letter, we use face-centered-cubic (fcc) aluminium as the reference material, for which many irradiation experiments and collision simulations results have been reported Jelinek et al. (2012); Gupta et al. (2016); Qiu et al. (2017); Trung et al. (2018); Mendelev et al. (2008); Botu and Ramprasad (2015); Botu et al. (2017), to validate the feasibility and reliability of this method.

In the DP-ZBL model, we assume that the system under consideration is composed of $N$ atoms with coordinates denoted by $\{\bm{R}_{i},\dots,\bm{R}_{N}\}$. Similar to the original DP model Zhang et al. ; Zhang et al. (2018), the DP-ZBL potential assumes the system energy is decomposed into atomic contributions, i.e.,

[TABLE]

with $i$ being the indexes of the atoms. The atomic contribution of atom $i$ is fully determined by the coordinates of atom $i$ and its near neighbors,

[TABLE]

where $s(i)$ denotes the chemical species of atom $i$, and $\mathcal{N}_{R_{c}}(i)$ denotes the set of near neighbors within cut-off radius $R_{c}$, i.e. $\mathcal{N}_{R_{c}}(i)=\{j|R_{ij}=|\bm{R}_{ij}|\leq R_{c}\}$. The atom contribution of DP-ZBL is the interpolation of the ZBL screened nuclear repulsion potential $E^{\textrm{ZBL}}_{i}$ and the standard deep potential $E^{\textrm{DP}}_{i}$

[TABLE]

where $w_{i}$ is the scale of ZBL potential that smoothly changes from 1 to 0 as the distance between atom $i$ and its nearest neighbor goes from 0 to a threshold value. To be more specific, the scale $w_{i}$ is defined as

[TABLE]

with $u_{i}$ being the short-hand notation defined by

[TABLE]

and $[R_{a},R_{b})$ denoting the range in which the ZBL potential and the deep potential are interpolated. It is noted that the switch function $-6u_{i}^{5}+15u_{i}^{4}-10u_{i}^{3}+1$ is continuous at 0 and 1 up to the second order derivative. The symbol $\sigma_{i}$ denotes the smooth-minimal distance of atom $i$’s near neighbors, which is defined by

[TABLE]

with $\alpha$ being a tunable scale of the distances between atoms. In the current work, we fix the scale to $\alpha$ = 0.1 Å.

As the schematic diagram is shown in FIG. 1, a ZBL-modified layer is added to better describe the strong repulsion at short inter-atomic distances in the DP-ZBL model, through the smooth switch function $w_{i}$ in Eq. 4. Then the DP-ZBL model neural network was trained with the same dataset generated by the deep potential generator (DP-GEN) in Ref. [Zhang et al., 2019], a scheme employing the idea of active learning Settles (2012) and reinforced dynamics Zhang, Wang, and E (2018b). Note that this dataset contains a vast range of configurations explored and labelled during the active-learning process, which ensures that the DP-ZBL model can be trained with enough possible configurations with high accuracy. Apart from the interpolation with the ZBL potential, the cut-off radius adopted by the current work is 6 Å, and the total training steps are 640, 000. These differences in training will not lead to a significant difference in the accuracy. We have also tested several switching ranges ($R_{a}$, $R_{b}$) to generate the DP-ZBL potentials. The best one (1.2Å, 2.0Å) was selected for collision cascades simulations that are presented below.

In order to evaluate the behavior of the DP-ZBL potential on irradiation effects, two classical potentials have also been employed for comparison, including the ZBL joined embedded atom method (EAM) potential (EAM-ZBL)Jacobsen, Norskov, and Puska (1987) and the state-of-the-art modified EAM (MEAM) potential Pascuet and Fernández (2015) with self-implemented ZBL (MEAM-ZBL), which are widely used in the previous irradiation simulations and give satisfying results Popok et al. (2009); Urbassek (1997); Lu, Li, and Brenner (2015). In this work, we used the DeePMD-kitWang et al. (2018a) for training the DP-ZBL potential, LAMMPSPlimpton (1995) for molecular dynamic simulations, VASPKresse and Hafner (1993, 1994); Kresse and Joubert (1999) for ab-initio calculations, and OVITOStukowski (2009) for the defect identification.

First, we have calculated some material properties using the three potentials and the DFT (see Methods in the supplementary materials), as summarized in TABLE 1. It is no wonder that the MEAM-ZBL potential provides nearly the same vacancy formation energies ($E_{\rm vf}$) as the results of experiments because these basic solid state properties have been used to tune the parameters of the MEAM potentials. Besides that, the DP-ZBL potential gives reliable results in all these considered properties. These results demonstrate that the DP-ZBL potential can still provide accurate predictions about the material properties near the equilibrium state with the accuracy comparable to the DFT calculations, which also implies that the smoothly joined ZBL potential which dominates the interatomic interactions below 1.2 Å would not influence the accuracy of the original DP model.

Next, we did collision cascade simulations by using these three potentials. Collision cascades are the feature phenomena in irradiation effects Wallace et al. (2018); Wang et al. (2018b); Zhou et al. (2018). When the energetic particles including protons, neutrons, electrons, and ions inject the target material, it will transfer energy to the target atoms. If the transferred energy is higher than the displacement threshold energy ($E_{\rm d}$) of target atoms, they will displace from the original lattice sites. If these primary knock-on atoms (PKAs) still have enough energy, they can knock out other target atoms subsequently and so on. Thus a large number of atoms are displaced from their original lattice sites, which is called the collision cascade. However, as the cascade begins to thermally equilibrate with its surrounding environment, most of the displaced atoms regain position in the perfect lattice structureMeldrum et al. (1998); Nordlund et al. (1999), as illustrated in FIG. 2.

It can be observed in FIG. 2 that all the three potentials exhibit a similar trend of displaced atoms during the evolution. The number of displaced atoms increased sharply within 1 ps and reached a peak at 0.3 $\sim$ 0.4 ps. Then it decreased monotonically because of the recombination process, and only a few defects remained. The evolution of displaced atoms generated by PKA at other energies was also illustrated in the Supplementary Materials (SM). It can be concluded from FIG. S5 that the peak value of displaced atoms increased with the increasing PKA energy, but the peak of MEAM-ZBL model is significantly higher than other models when the PKA energy is larger than 2 keV. Though this transient process can hardly be examined by experiment or other models, which means we cannot give a reliable estimation which one is more accurate, we can still conclude that the collision cascade evolutions provided by DP-ZBL potential do not significantly deviate from the results obtained other existing methods.

In fact, the number of residual point defects is even more important than the peak value during evolution, for a broad range of fundamental science and applied engineering applications. To quantify the numbers of point defects caused by a single PKA, Norgett $et~{}al.$ have proposed the Norgett-Robinson-Torrens (NRT) model, based on the binary collision approximation method, to evaluate the bombardment damage Robinson and Torrens (1974); Norgett, Robinson, and Torrens (1975a). However, it has been recognized for several decades that the NRT model overestimates near 3 times the number of stable defects in pure metals after energetic cascades Averback, Benedek, and Merkle (1978); Zinkle and Singh (1993); Konings (2011). Therefore, we calculated the residual point defects by the NRT model and used one third of it as a benchmark to evaluate our cascade simulation results, and a brief introduction of NRT model was also introduced in the SM.

As shown in FIG. 3, the residual displaced atoms calculated by the three potentials all exhibit similar near-linear trends with the value of initial PKA energy, which is consistent with the NRT model according to Eq. S2. Note that the slope of the fitting line is inversely proportional to the displacement energy $E_{\rm d}$, which is usually defined as the minimum amount of energy transferred to a lattice atom to make it displace from original stable site Was (2016). Other widely used models, such as the Kinchin-Pease (KP) model Kinchin and Pease (1955) and the athermal recombination corrected DPA (arc-DPa) model Nordlund et al. (2018), typically take $E_{\rm d}$ as a key parameter to quantify the amount of displacement damage generated by energetic particles inject in materials. So we also made a comparison of average $E_{\rm d}$ value in all the possible directions calculated by these three models, as shown in TABLE 2 and FIG. S4.

According to our calculations, $E_{\rm d}$’s provided by the DP-ZBL and MEAM-ZBL potentials were quite close to the recommended value (25 eV for fcc Al) Norgett, Robinson, and Torrens (1975b); Was (2016), while the result of EAM-ZBL potential has a relatively large deviation. Furthermore, the DP-ZBL potential can also give the best prediction of residual displaced atoms in the three examined potentials, if we took the 1/3 NRT results as a reference. Therefore, the simulation results of DP-ZBL potential are consistent with most of the existing theoretical models in the field of low-energy ion irradiation damage.

Meanwhile, the displacement energy of three potentials are in the order of $E_{\rm d}$(DP-ZBL) > $E_{\rm d}$(MEAM-ZBL) > $E_{\rm d}$(EAM-ZBL), so the corresponding residual point defects values should be $N_{\rm d}$(EAM-ZBL) > $N_{\rm d}$(MEAM-ZBL) > $N_{\rm d}$(DP-ZBL), which is completely in accordance with the results in FIG. 3. Besides, it is worth noting that the DP-ZBL model usually produces smaller $N_{\rm d}$ than the classical EAM-ZBL and MEAM-ZBL potentials, which may be caused by the fact that the DFT data sets used to train the DP-ZBL model have considered the energy difference of configurations far from the equilibrium state, while the traditional EAM-ZBL and MEAM-ZBL potentials are simply constructed by fitting the material properties near the equilibrium state. From this point of view, the DP-ZBL model should provide a better description for the irradiation damage events than traditional analytical potentials.

In conclusion, we have proposed the DP-ZBL scheme by smoothly interpolating the accurate repulsive pair potential (ZBL) into the DP model. The DP-ZBL potential generated in this way can not only give accurate results regarding the material properties near equilibrium states but also be sufficient to describe the atomic collision cascades during the irradiation damage processes. This method can minimize the impact of subjective factors on potentials during their establishments, and provide higher agreement with experimental or DFT results compared with other widely used classical potentials. Moreover, due to the applicability of the DP model to a wide range of materials, see, e.g. Ref.[Zhang et al., 2018], this newly proposed method may be used in the irradiation effect studies of new advanced materials, such as high entropy alloys (HEAs), layered transition metal ternary nitrides and carbides ($\rm{M_{n+1}AX_{n}}$ phases), and two-dimensional (2D) materials, for whom suitable classical potentials are still lacking. We hope that with the development and improvement of the DP-ZBL potential database and corresponding computational algorithm, the predicament of accurate potentials lacking in the irradiation effect field could be better solved.