General-purpose neural network interatomic potential for the $\alpha$-iron and hydrogen binary system: Toward atomic-scale understanding of hydrogen embrittlement

To understand the physics of hydrogen embrittlement at the atomic scale, a general-purpose neural network interatomic potential (NNIP) for the $\alpha$-iron and hydrogen binary system is presented. It is trained using an extensive reference database produced by density functional theory (DFT) calculations. The NNIP can properly describe the interactions of hydrogen with various defects in $\alpha$-iron, such as vacancies, surfaces, grain boundaries, and dislocations; in addition to the basic properties of $\alpha$-iron itself, the NNIP also handles the defect properties in $\alpha$-iron, hydrogen behavior in $\alpha$-iron, and hydrogen-hydrogen interactions in $\alpha$-iron and in vacuum, including the hydrogen molecule formation and dissociation at the $\alpha$-iron surface. These are superb challenges for the existing empirical interatomic potentials, like the embedded-atom method based potentials, for the $\alpha$-iron and hydrogen binary system. In this study, the NNIP was applied to several key phenomena necessary for understanding hydrogen embrittlement, such as hydrogen charging and discharging to $\alpha$-iron, hydrogen transportation in defective $\alpha$-iron, hydrogen trapping and desorption at the defects, and hydrogen-assisted cracking at the grain boundary. Unlike the existing interatomic potentials, the NNIP simulations quantitatively described the atomistic details of hydrogen behavior in the defective $\alpha$-iron system with DFT accuracy.

Figure 1

Comparison of DFT and NNIP (a) energies and (b) forces of the structures in the training and testing data sets. The line with a slope of 1 corresponds to a perfect training.

Figure 2

NNIP performance for $\alpha$-iron. (a) Phonon dispersion curve, (b) $\gamma$ surface for the $\left(112\right)$ crystallographic orientation, (c) misorientation-energy relationship for the symmetric tilt GB with $\langle 110\rangle$ tilt axis, and (d) 2D Peierls potentials. The axes along [101] and $\left[11\overline{2}\right]$ directions in (d) are scaled by $\sqrt{3}{a}_0/6$; $a_0$ is the lattice constant of $\alpha$-iron. The available reported DFT [16, 50] and experimental [51] results for phonon dispersion curves and DFT results for symmetric tilt GB formation energy [52, 53, 54] are also shown.

Figure 3

(a) Migration barriers of the hydrogen-vacancy complex. DFT and EAM results extracted from Ref. [10] and Ref. [81] are also plotted. Locations of vacancy and H atom indicated by ${\mathrm{V}}_m+{\mathrm{T}}_n(m=0$–1 and $n=0$–5), where ${\mathrm{V}}_m$ and ${\mathrm{T}}_n$ indicate the location of vacancy and H atom, respectively, and the schematic diagram is presented in the inset of (a). (b) Hydrogen dissociation on the (001) surface. Positions of hydrogen atoms in several interesting locations are depicted in the inset of (b).

Figure 4

(a) Atomic structure of $\mathrm{\Sigma }5$ GB showing several sites of interest at and around the GB plane. Red dotted line indicates the predesigned fracture path, and grain-1 and grain-2 represent two separated grains after fracturing. (b) Solution energies of H atom at various sites ($E_{\mathrm{s}}^{\mathrm{GB}}$) presented in (a). (c) Diffusion barriers of the H atom along the $\left[1\overline{3}0\right]$ direction in the GB plane. DFT results extracted from Ref. [96] are plotted in (b).

Figure 5

(a) Positions of H solution sites around the easy core configuration of screw dislocation. The $E_{\mathrm{k}}$ site is equivalent to $E_1$ and is used to calculate the influence of the H atom on kink pair nucleation energy (see text for details). The energy basins are shaped like curvy caterpillars. The side view of (a) is displayed in (b). (c) Solution energies of the H atom at various sites marked in (a) and T site in bulk, presented with DFT results from Ref. [98]. (d) Diffusion barriers of H atom in the paths marked in (a); the solution energy of the H atom at the T site is set as the benchmark, which is shown by the dotted line. (e) $\gamma$ surface of pure iron with and without H atom. (f) Kink nucleation energy of dislocation with and without H atom (at $E_{\mathrm{k}}$).

Figure 6

(a) Snapshots of the system at the simulation time of 0 ps. Several kite units are designed in the $\left[1\overline{3}0\right]$ direction and their sequence number is marked next to it individually. The large view of kite unit and H-atom diffusion path is also presented. (b) Snapshot at 150 ps. Trajectories of H atoms of interest, i.e., atoms ${\mathrm{H}}_{\mathrm{O}}$ (olive line) and ${\mathrm{H}}_{\mathrm{R}}$ (red line), at the simulation time interval of 0–40 ps, ${\mathrm{H}}_{\mathrm{B}}$ (black line) at 100–130 ps are presented. (c) H volume density (${\rho }_{\mathrm{H}}$) evolution with time in each kite unit marked in (a). (d) The screw dislocation model used in the simulation. The fixed and free Fe atoms are indicated by dark red and dark gray balls, respectively, and the three columns of blue balls denote the dislocation core. (e) Three typical trajectories of H atoms, presented in a MD simulation for 130 ps and viewed from two directions. (f) Interior trajectories of H atoms in the H-charging simulation; all iron atoms are set to be invisible. The thick green line in (d)–(f) indicates the dislocation line obtained by the DXA tool implemented in the ovito code [107].

Figure 7

(a) Radial distribution function (RDF) of volume densities of the number of H atoms around vacancy, GB, and screw dislocation at various temperatures. (b) Volume density of the number of H atoms around defects within 2 Å [dotted lines, replotting of the first column data of (a)] and the hydrogen desorption rates at various temperatures (solid lines).

Figure 8

Tensile tests for clean and H-charging $\mathrm{\Sigma }9$ twists. (a) Clean $\mathrm{\Sigma }9$ twist model. (b) Left, middle, and right panels present the tensile stress, dislocation volume density (${\rho }_{\mathrm{Dis}}$), and H atomic concentration (${\rho }_{\mathrm{H}}^a$) in the models versus the tensile strain. (c) Volume density of H (${\rho }_{\mathrm{H}}^V$) in the H-charging twist at various tensile strains. (d), (e) Snapshots of clean and H-charging twists at various tensile strains (Fe atoms with bcc lattice are set as invisible), and the corresponding strains are printed under each panel. Blue and red dots denote Fe and H atoms, respectively. The green surface and thick lines in (d) and (e) are the regions that differ from bcc lattice [113] and dislocations, respectively, which are analyzed by the DXA tool implemented in the ovito code [107].