General-purpose neural network interatomic potential for the α -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 α -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 α -iron, such as vacancies, surfaces, grain boundaries, and dislocations; in addition to the basic properties of α -iron itself, the NNIP also handles the defect properties in α -iron, hydrogen behavior in α -iron, and hydrogen-hydrogen interactions in α -iron and in vacuum, including the hydrogen molecule formation and dissociation at the α -iron surface. These are superb challenges for the existing empirical interatomic potentials, like the embedded-atom method based potentials, for the α -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 α -iron, hydrogen transportation in defective α -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 α -iron system with DFT accuracy.


I. INTRODUCTION
The hydrogen-induced mechanical performance degradation of metals and all phenomena associated with its deleterious effects are known as hydrogen embrittlement (HE) [1,2], which is a prominent problem in current societal and technological requirements [3]. Iron and its alloys are the most widely used structural materials in various industries due to their excellent mechanical properties and low cost. The development of iron alloys with HE resistance requires a deep understanding of the interaction between hydrogen and the defects in the iron alloy's matrix [4]. Experimental examination of the role played by hydrogen in metals is extremely challenging owing to its high diffusivity, low solubility, and small size [5]. As an alternative solution, computer simulations at the atomic level can provide insightful information. Electronic structure calculations, such as density functional theory (DFT), supply the most accurate potential energy surfaces (PESs). However, the high computational cost limits their application to large-scale simulations despite the remarkable progress made in computational power and algorithms. Unlike DFT calculations, the molecular dynamics (MD) method with empirical interatomic potentials (EIPs) can simulate large systems at long time scales, but they strongly rely on the quality of EIPs; the accuracy and transferability of EIPs are always questionable [6].
Several EIPs, including the embedded atom method (EAM) and modified EAM (MEAM) potentials, for H atoms in α-iron have been proposed [7][8][9][10] and have made remarkable contributions to the understanding of HE. Among these EIPs, two EAM potentials for the Fe-H system [10] based on iron potentials, proposed by Mendelev et al. [11] and Ackland et al. [12], can produce the corrected atomic structure of screw dislocation. Unfortunately, they cause strong attractive interaction among H atoms, causing H atoms to cluster in the matrix at moderate temperatures. Although these potentials have been improved by introducing an additional repulsive term for H-H interactions [13] or by refitting the parameters of the Fe-H and H-H interactions [14], the problem remains in scenarios with high hydrogen concentration. Another severe drawback of such potentials is their poor description of the energetics of screw dislocation [15,16]. More information on the existing Fe-H EIPs is provided in the latest review paper on this topic [17].
Consequently, an interatomic potential for the α-iron and hydrogen binary system that combines the advantages of the accuracy of DFT calculations and computational cost of EIPs needs to be constructed. Advances in machine learning methods allow the attainment of this goal. In principle, machine-learning-based interatomic potentials (MLBIPs) can approximate any function with arbitrary accuracy [16]. ML-BIPs, such as neural network interatomic potentials (NNIPs) [18] and Gaussian approximation potentials (GAPs) [19], have been proposed and successfully applied to many systems [20][21][22][23]. Recently, Davidson et al. [24] proposed a purposemade H-Fe Gaussian approximation potential that enables rapid sampling of H atoms in a vacancy with DFT-level accuracy, and reported novel results, while the potential does not contain the Fe-Fe interaction. Although this is a wise strategy for constructing a purpose-made MLBIP that solves a specific target problem, a target involving multiple defects requires a general-purpose MLBIP.
In this study, the NNIP framework was adopted and a general-purpose NNIP for the α-iron and hydrogen binary system was constructed. The constructed NNIP can well reproduce the interactions of hydrogen with various defects in α-iron, such as vacancies, surfaces, grain boundaries, and dislocations; in addition to the basic properties of α-iron itself, the NNIP can also handle the defect properties, hydrogen behavior, and hydrogen-hydrogen interactions in α-iron and in vacuum, including the hydrogen molecule formation and dissociation at the α-iron surface. These properties are essential for understanding the physics of HE, which cannot be fully described by the existing interatomic potentials. Eventually, the constructed general-purpose NNIP was applied to several key phenomena underlying HE in α-iron with point defect, dislocation, surface, and grain boundaries. These phenomena include hydrogen charging and discharging, hydrogen transportation, hydrogen trapping and desorption at defects, and hydrogen-assisted cracking at the grain boundary. The simulations uncovered the atomistic details of these phenomena in DFT accuracy.

A. Reference database and neural network architecture
To learn and reproduce is the essence of MLBIPs [25,26]. The NNIP quality relies on the reference database, and NNIP transferability depends to a certain extent on the similarity of structures presented in actual simulations against those in the reference database. A broad spectrum of structures but limited to important structures that directly related the target properties were considered. In total, 21 928 configurations, equivalent to 1.9 × 10 6 atomic environments, were prepared in the reference database (see Sec. SI of the Supplemental Material [27] for details). Both the total energy and three force components of each atom were included for each structure [28]. The structures in the reference database were randomly distributed into training data set (90%) and testing data set (10%), which were used to fit the adjustable parameters and check the quality of the potential for unlearned structures. The neural network (NN) architecture has two hidden layers, each of which possesses 15 neurons (see Fig. S2 in the Supplemental Material [27]). The activation function is set as a logistic form in both hidden layers. Sixteen radial and 48 angular atom-centered symmetry functions (ACSFs) [29] were adopted for each element (see Sec. II B for details of the ACSF settings); 1231 adjustable parameters (i.e., 1200 weights and 31 biases) were used for each element. The entire NNIP training task was completed using the neural network potential package (n2p2) [30]. The constructed NNIP can be applied in the LAMMPS code [31] via the interface implemented in the n2p2 package [32], and both the reference database and potential files are shared online [33].

B. Calculation settings
ACSFs settings. Two types of ACSFs [29], namely, radial symmetry functions and angular symmetry functions, were adopted to distinguish the local atomic environments within the cutoff radius. The radial symmetry function is defined as where R i j is the atomic distance between an atom j and the central atom i. The angular symmetry function is given by where θ i jk is the angle enclosed by the vectors of R i j and R ik of two neighboring atoms j and k, respectively. Both types of ACSF have a common function f c called the cutoff function, which is defined as follows: Here, R c is the cutoff radius. The parameters η, R s , ξ , and λ in Eqs. (1) and (2) can be determined by the strategy reported in Ref. [34]. The cutoff radii are 6.5 and 3.6 Å for Fe and H atoms, respectively. DFT calculation settings. The reference database was produced by spin-polarized DFT calculations using the VASP code [35]. The electron-ion interaction was described using the projector augmented wave [36] method. The exchangecorrelation energy was treated with the GGA-PBE approach [37]. The Brillouin zone (BZ) was sampled using the Monkhorst-Pack scheme, and the BZ integration employed the Methfessel-Paxton method for relaxations. The smallest allowed spacing between k-points was 0.03 Å −1 , and the energy smearing was 0.1 eV. A cutoff energy of 360 eV was employed to truncate the plane-wave expansion of the wave functions. Geometry optimizations were continued until the force in the system was less than 10 −2 eV/Å. The stopping criterion for the electronic structure optimization was set to 10 −5 eV. To enlarge the coverage area of the PES, the reference database was expanded by ab initio MD (AIMD) sampling. AIMD was performed in the canonical ensemble with a Nosé-Hoover thermostat implemented in the VASP code. All systems were simulated within 100-1000 K (up to 2000 K in liquid configurations only) with a time step of 1-2 fs. The transition state between the initial and final atomic configurations was calculated using the climbing image nudge elastic band (CI-NEB) method. The zero-point vibration energy was not included, but if needed, can be incorporated (for example) by the path-integral method [38] when using the NNIP in a simulation.

A. Overall NNIP accuracy
Owing to the possibility of trapping in local minima in the NN learning process, many NNIP trial trainings should be performed with different initial states, and many candidates should be obtained. Among these candidates we should choose the best NN that well describes not only the training data set but also the testing data set. The NNIP quality is estimated by the root-mean-squared error (RMSE) of energies (E ) and forces (F ) primarily used to guide the optimization of adjustable parameters: where N s and N a denote the number of structures in the training and testing data sets and the number of atoms in individual structures, respectively. The NNIP with the highest quality among all candidates has a low RMSE and good reproducibility and transferability of the important properties; this was adopted and addressed in detail in the following text. As shown in Fig. 1, the training and testing data sets have RMSE(E )s of 2.98 and 3.01 meV/atom, and RMSE(F )s of 69.68 and 68.16 meV/Å, respectively. These RMSEs are sufficiently small because they approach the usual accuracy of DFT analysis, and the small difference between the RMSEs of training and testing data sets denotes that no serious overfitting occurred [39]. Just in case, an additional data set including 3240 configurations, produced by AIMD simulations at higher temperature for the different systems, was prepared and further confirmed the no-overfitting issue (see Sec. SIII of the Supplemental Material [27]). All points for energy and force are very close to the line with a slope of 1, implying that the trained data can be well reproduced. Several concerned properties are further confirmed in the following sections. These properties following the data source categories are fully included, partially included, and not included, and also of the test results are tabulated in Table S3 in the Supplemental Material [27].

B. NNIP performance for pure α-iron
The NNIP of the Fe-H binary system must be reliable also for α-iron. The basic properties of iron, such as lattice constant [40][41][42][43], elastic constants [44], vacancy cluster formation energies [45][46][47], self-interstitial formation energies [48,49], low index surface energies [50], Bain path, and energy-volume relations [16] are summarized in Table S4 and Fig. S4 in the Supplemental Material [27]. All of the above-mentioned properties are in good agreement with the DFT results, and most of them are included in the reference database (see Sec. SIV of the Supplemental Material [27]).
Phonon dispersions of α-iron under an equilibrium lattice condition with the supercell size of 3 × 3 × 3 were computed using the PHONOPY code [55] and displayed in Fig. 2(a). The results produced by NNIP show an excellent agreement with DFT [16,50] and with the experimental [51] results. The training data set contains the energy and force information of the AIMD snapshots at finite temperatures, which is sufficient to obtain the phonon and vibration properties of iron.
γ surfaces, or generalized stacking fault (GSF) energies [56], of the (110) and (112) planes (the dominant slip planes in bcc metals [57]) were considered. The γ surface was calculated by in-plane shifting between upper and lower halves in the supercell. In-plane atomic motions were fixed while those in the normal direction were allowed. The results are displayed in Fig. 2(b) and Fig. S4(f) in the Supplemental Material [27]. We confirmed that the NNIP can reproduce the reported γ surface [50]. Note that as mentioned above, the γ surface energies are included in the reference database.
Grain boundaries (GBs) are another important defect for the polycrystalline plasticity and fracture. Therefore, it is strongly expected that the NNIP can describe various GB structures and energetics. If the NNIP can well describe the 113606-3  [101] and [112] directions in (d) are scaled by √ 3a 0 /6; a 0 is the lattice constant of α-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. GB structure and energetics, it should also well describe edge dislocations because low-angle GBs are constituted by an array structure of edge dislocations. We mapped the misorientation versus the formation energy of symmetric GBs with tilt axes of 110 , 001 , and 111 . The result for the 110 tilt axis is plotted in Fig. 2(c). As marked by red circles in Fig. 2(c), 10 GBs with misorientation in the range of 20 • to 148 • are included in the reference database. Compared to the DFT results, NNIP accurately reproduces the formation energies of the trained GBs, and, more importantly, it correctly predicts the formation energies of other GBs [52][53][54], which are not included in the reference database. The maps for GBs with tilt axes of 001 and 111 also show good agreement with DFT results (see Figs. S4(d) and S4(e) in the Supplemental Material [27]).
Screw dislocation motion is essential for the plasticity of bcc metals [58], and the core structure of screw dislocation and motion of kink pairs is believed to be attributed to that behavior [59][60][61]. To examine the core structure of screw dislocation, we used a supercell with dimensions of 137.4, 2.45, and 40.3 Å in the [112], [111], and [110] directions, respectively, comprising a total of 1200 Fe atoms. Periodic boundaries are taken along the [112] and [111] directions. A 1 2 [111] screw dislocation with easy core configuration is embedded at the center of the simulation cell. To study the dislocation motion, a large model, that is 40 times larger than the 1200-Fe-atom model in the [111] direction, was constructed (see Fig. S5(a) in the Supplemental Material [27]). This periodic array of dislocation (PAD) configuration is widely used in atomic simulations [62,63]. Our NNIP can correctly predict the atomic structures of easy core, hard core, and split core (differential displacement maps for the three core configurations are shown in Figs. S5(b)-S5(d) in the Supplemental Material [27]) and can also satisfactorily predict their energetics. The two-dimensional (2D) Peierls potential [ Fig. 2(d)] was obtained using numerical interpolation from the corresponding energy profiles along the pathways among the above-mentioned atomic configurations; in total, 28 configurations were adopted. The 2D Peierls potential vividly matches that reported in Refs. [15,64], although only the data on the straight paths from one easy core state to a neighboring easy core state and from hard core to split core states are included in the reference database using a small 150-Fe-atom model with various screw dislocation core configurations at center (see Fig. S1(a) in the Supplemental Material [27]). Clearly, the Peierls barrier between the two neighboring easy core configurations has only one hump, which has important consequences on the kink pair formation mechanism [65] and most of the EAM potentials cannot reproduce the onehump shape [50]. Compared to the energy of the easy core configuration, the energies of the configurations of the hard core, split core, and the middle point, respectively, are 47.4, 82.3, and 38.2 meV/b, which reasonably agree with DFT results of 39.3, 108, and 37.9 meV/b [15], and 33.2, 87.9, and 34.9 meV/b [64], and 57.9, 110.3, and 49.2 meV/b [66], where b is the Burgers vector length, i.e., √ 3a 0 /2 (a 0 is the lattice constant of α-iron). Screw dislocation gliding and kink pair nucleation under shear stress were further tested (see Sec. IV of the Supplemental Material [27]).

C. Hydrogen solubility and diffusivity in bulk α-iron
The solution energy and diffusion barrier of the H atom should be reproduced to successfully demonstrate the hydrogen charging and diffusion phenomena. In α-iron, H atoms preferentially settle at two interstitial sites: the tetrahedral (T) and octahedral (O) sites. The solution energy of the H atom in bulk, E s , is calculated as follows: where E Fe+H and E Fe are the energies of Fe bulk with and without H atom at the interstitial site, respectively, and E H 2 is the energy of a gaseous H 2 molecule. The E s of a H atom at the T site obtained by NNIP is 0.236 eV, which reasonably agrees with the DFT results of 0.22 [67], 0.234 [68], and 0.21 eV [69]. The E s for the H atom at the O site is determined as 0.389 eV using NNIP, which matches the reported DFT result of 0.35 eV [70]. The diffusion energy barrier, E , for H-atom diffusion between neighboring T sites predicted by NNIP is 0.108 eV, which is in line with our DFT result of 0.092 eV as well as other DFT and experimental results of 0.088 [71], 0.096 [72], and 0.035-0.14 eV [53]. The distance from the saddle point in the T-T diffusion path to the nearest O site is 0.368 Å, which is consistent with a previously reported DFT result of 0.407 Å [71]. All of the above data are included in the reference database. The diffusivity of H atoms in α-iron bulk [73][74][75] and the volumetric engineering strain dependent hydrogen solution energy and diffusion barrier were further tested and provided in the Supplemental Material (Secs. SV and SVI) [27].

D. Hydrogen-vacancy interaction
Hydrogen trapping by defects is an essential phenomenon in hydrogen embrittlement of structural metallic metals [76] because it leads to an extremely high localized H concentration [77], then leads to the defects growth, such as void or bubble formation [67] and crack initiation and propagation. First, we systematically confirmed the interaction between monovacancy and various number of H atoms, which demonstrated fairly good agreement with the reference database and DFT results in the literature [78][79][80][81]. The details of the trapping configurations and trapping energies can be found in the Supplemental Material (Sec. SVII) [27].
Next, we paid more attention to the diffusion of the Hvacancy complex because it is still challenging for empirical potentials [81]. The energy profile of the H-vacancy complex in the diffusion pathway outlined in Ref. [81] is predicted using the CI-NEB method with NNIP, as shown in Fig. 3(a). The diffusion barrier for the H-vacancy complex (V 0 + T 1 → V 1 + T 4 ) is 0.837 eV, which is higher than that of a H-free monovacancy of 0.73 eV. The corresponding DFT results are 0.79 and 0.69 eV, respectively [81]. The EAM potential can reasonably predict the diffusion barrier, 0.78 eV; however, the details of the energy profile differ from the DFT result [81], especially the double hump presented at the final vacancy jump (V 0 + T 4 → V 1 + T 4 ). In addition, the barrier that the H atom has to overcome to jump from vacancy to a neighboring T site (V 1 + T 3 → V 1 + T 5 ) is 0.591 eV; this is very close to the DFT results of 0.554 [81], 0.583 [72], and 0.650 eV [10]. The barrier for the H atom to jump into a vacancy from a neighboring T site (V 1 + T 5 → V 1 + T 3 ) is 0.069 eV; this is also in agreement with the DFT results of 0.078 [81], 0.09 [72], and 0.033 eV [10]. Note that these energy barriers are included in the reference database.
More importantly, it is difficult for the EAM potentials to predict the diffusion barrier between trapping sites in a vacancy [10], such as V 0 + T 0 → V 0 + T 1 [inset of Fig. 3(a)]. By virtue of the flexibility of NNIP, the above-mentioned energy barrier predicted by NNIP is 0.245 eV, which is in close agreement with the DFT result of 0.230 eV, whereas the EAM result is 0.061 eV [10]. Interestingly, the energy barrier of 0.245 eV is significantly larger than that recorded in bulk of 0.108 eV, even when the free volume in vacancy is larger than that in bulk. This could be explained by the fact that the H atom at T 0 (abbreviated as H@T 0 ) bonded with its neighboring Fe 1 atom through a covalent bond of Fe 3d-H 1s [78].
To accomplish a jump from T 0 to T 1 , the Fe 1 -H@T 0 covalent bond must be broken, for which a large energy barrier needs to be conquered before forming the Fe 2 -H@T 1 bond. This suggests that NNIP can accurately describe the bond-breaking and bond-forming processes.

E. Hydrogen-surface reaction
The interactions of hydrogen with surfaces are always engaged in the topic of HE, such as the hydrogen charge and uncharge process, nanoscale cavity migration in a H 2 -gas environment [82]. The adsorption of H atom or H 2 molecule, and H 2 molecule formation and dissociation processes in iron surfaces has been extensively studied using DFT [83][84][85][86]. However, the description of these process remains difficult for the empirical potentials. The adsorption energies of a H atom on several low-index surfaces were reproduced with NNIP and show good agreement with the DFT results; these results are tabulated in Table S6 in the Supplemental Material [27]. Subsequently, the H 2 molecule dissociation process on the (001) surface is further demonstrated via CI-NEB calculations. The potential energy profile along the formationdissociation path and several typical atomic configurations are shown in Fig. 3(b).
To compute the H 2 formation-dissociation energetics and pathway, we used a supercell with a size of 3 [ Fig. 3(b), A], both DFT and NNIP predicted that the configuration of those H atoms at hollow sites with a separation of d H−H = √ 2a 0 (a 0 is the lattice constant) is energetically preferred. Thus, this configuration is set as the initial state. The final configuration is obtained by relaxing the configuration with the two H atoms at 6.0 Å above the surface [see Fig. 3 [87] and the experimental result of 0.187 eV [88]. The transition from state C to state E has the energy barrier of 0.268 eV. The energy barrier for H 2 molecule transition from state E to state C (D) is about 0.1 eV, indicating the repulsive interaction between the molecule and substrate. A movie for the H 2 molecule formation process on the (001) surface at 800 K can be found in the Supplemental Material and is named Movie S2 [27].
Note that most of the above hydrogen-surface reaction data are not included in the reference database, while structures from AIMD for low-index surfaces (see Table S2 in the Supplemental Material [27]) with various number of H atoms are included in the reference database.
First, the solution energy of hydrogen to interstitial sites [see Fig. 4(a)] in 5 GB (E GB s ) was calculated by Eq. (6). NNIP predicted that site 1 (short for s1, the same below), which is the equivalent of s4, is the most stable solution site having an E GB s of −0.180 eV; this well agrees with the DFT result of −0.187 eV [96,97]. The solution energies for other sites are presented in Fig. 4(b) and are found to align well with the results of our DFT calculations, although the solution energies were not fully included in the reference database.
The energy profile of H-atom diffusion toward the GB and in the GB was studied via CI-NEB calculations, which was not included in the reference database. The following three diffusion paths were considered: (a) between two adjacent most stable solution sites in the [001] direction, i.e., s1 ↔ s1 [001] ; (b) from the most stable solution site to a T site in bulk, i.e., s4 ↔ s8; and (c) along the entire GB, i.e., s1 ↔ s2 ↔ s3 ↔ s4.
Each of the above-mentioned three paths has been reported in previous studies [91,96], with the exception of the part of s2 ↔ s3 ↔ s4 in path (c). The unavailable part in path (c) was newly calculated in this study. The diffusion energy barrier for path (a) with NNIP is 0.237 eV, which well agrees with the reported DFT calculations of 0.250 eV [91,96]. For the second path, two humps in the diffusion path yield energy barriers of 0.451 and 0.568 eV, respectively, which also agree with the previously reported DFT results of 0.40 and 0.57 eV [96]. The diffusion pathway of path (c) is shown in Fig. 4(c). NNIP can properly describe the entirety of the migration process. The large energy barrier between s2 and s3 of 0.674 eV denotes that it is extremely difficult for the H atom to migrate across the entire GB. The energy barrier from s4 to s6 with NNIP is 0.165 eV, i.e., significantly lower than that of s2 to s3. The trapped hydrogen atom prefers to stay in the region enclosed by the sites of s4 and s6 and their equivalent sites in the [001] and [310] directions. This qualitatively meets the conclusion in the previous study [91].

G. Hydrogen-screw dislocation interaction
The presence of hydrogen around screw dislocation is believed to change its motion because the trapping, such as a segregation isotherm, stabilizes the screw dislocation by reducing total free energy. Therefore, the solution energies of the H atom at various sites around the screw dislocation (using the large PAD model described in Sec. III B and presented in Fig. S5(a) in the Supplemental Material [27]) should be reproduced with NNIP. Using NNIP we calculated them and compared with the available DFT results, as shown in Figs. 5(a)-5(c). Note that DFT data of the hydrogen-screw dislocation interaction using the small 208-Fe-atom model with an easy core configuration at the center (see Fig. S2(d) in the Supplemental Material [27]) are included in the reference database. The positions of Fe atoms in the outer layers are fixed to sustain the core structure of dislocation. The H atom is placed at various sites around the dislocation to represent the different binding sites in Fig. 5(a). The H-atom solution energy in the screw dislocation calculated with NNIP is shown in Fig. 5(c). The energetically favorable sites can be found as E 1 and E 2 . The solution energy for the H atom at which the site of C m (m = 1-3), is 0.13 eV; nonetheless, it is still energetically preferred over bulk sites. All solution energies are within small error ranges compared with the DFT calculations in Ref. [98] [see Fig. 5(c)].
The diffusion of interstitial solute along the screw dislocation is also of interest [99] because the trapped H atoms may mainly diffuse along the screw dislocation (this diffusion mechanism is called pipe diffusion). Thus, the diffusion barriers around the screw dislocation determine the diffusion pathway and diffusivity along the screw dislocation. As expected in Ref. [98], the region enclosed by E 1 and E 2 is an energy basin that the H atom can move with a tiny barrier of 0.048 eV. The barrier for an H atom to enter the core region of dislocation from the energy basin is 0.15 eV, whereas that to diffuse along the dislocation inside the core is only 0.057 eV, i.e., much lower than the path from the E 1 site to the neighboring E 1 site along the Burgers vector, which is 0.638 eV. The actual diffusion path at finite temperature can be observed in MD simulation for the random walking H atoms. This is shown in the following application section, whereas other considered diffusion barriers are shown in Fig. 5(d). Data on H-atom diffusion in a model with screw dislocation are not included in the reference database.
The effect of H on the core shape and mobility of screw dislocation can be indirectly described using the γ surface. In the presence of an H atom at the T site in the sliding interlayer, the maximum of the γ surface along the 111 direction is 5 meV/Å 2 lower than that of pure iron, as shown in Fig. 5(e), which is consistent with the DFT results [100]. The calculated kink pair nucleation energy at zero external stress for dislocation with one H atom at the site of E k [see Fig. 5(a)] located in dislocation gliding direction [112] is displayed in Fig. 5(f). We find that the nucleation energy is lowered because of the attractive interaction between screw dislocation core and solute H atom located in the gliding direction [inset of Fig. 5(f)], which is qualitatively consistent with the prediction using the LT+DFT model [98]. No data regarding the kink pair nucleation energy are included in the reference database. The reduction in the maximum of the γ surface along the 111 direction suggests that the presence of an H atom enhances the intrinsic mobility of dislocation; however, simultaneously the dragged effect by the presence of H atoms at the dislocation core should be considered, which may depress the mobility of dislocation and depends on the H concentration and temperature [101]. This requires further study.

IV. APPLICATIONS
Hydrogen-enhanced decohesion (HEDE) [102], hydrogenenhanced localized plasticity (HELP) [103], hydrogenenhanced strain-induced vacancies (HESIVs) [104], and other HE mechanisms have been proposed in the literature [105], and posterior interpretations of experimental observations have been made. Nevertheless, the physics of HE remains poorly understood. It has been suggested that all hydrogen-113606-7 induced embrittlement processes cannot be driven by any single mechanism [106]. Therefore, this section demonstrates several key phenomena concerning the HE, such as hydrogen charging and diffusion to defective α-iron, temperaturedependent hydrogen trapping and desorption to different defects in α-iron, and plasticity and hydrogen-enhanced GB fracture by MD simulations with the constructed NNIP.

A. Hydrogen charging and diffusion processes in GB and screw dislocation
The first step of H 2 gas permeation into metals is dissociation of H 2 molecules on the metal surface. H 2 molecules also dissociate and form on the surfaces of crack tips and nanocavities, which are all critical processes for HE. In this section, the charging process and diffusion behaviors of H atoms in tilt GB and screw dislocation of α-iron are simulated.
The hydrogen charging and diffusion behavior in a tilt GB were simulated using a bicrystal model with a 5 symmetric tilt GB [see Fig. 6(a)  (d) Perform MD runs and calculate the average for ρ H estimation. Figure 6(b) shows the snapshot of the bicrystal model simulated for 150 ps. The figure reveals that the H atoms prefer diffusing along the GB and the H concentration in the GB is, as expected, significantly higher than that in the bulk region. The time evolution of ρ H is shown in Fig. 6(c) and suggests that H atoms gradually diffuse from the two surfaces to the interior part; i.e., GB is a ditch for H transportation. The H diffusion path can be directly depicted using the H traveling trajectories, and three typical diffusion manners are shown in Fig. 6(b) Fig. 6(b). This trajectory indicates that H atoms tend not to diffuse directly through the GB unit but instead deviate to region C before jumping to region A in the next kite unit, as schematically illustrated in the subpanel of Fig. 6(a). With an increase in hydrogen concentration in the GB region, the GB serves as a source of hydrogen migration into the bulk, such as the migration path of H B presented in Fig. 6(b). Additionally, hydrogen segregation at and diffusion along the GB is a process strongly relying on the angle between grains [89]; one can directly observe the hydrogen diffusion path in any other GBs via this charging model with the NNIP.
The diffusion of H atoms in screw dislocation was further studied. The atomic model used in simulation is shown in Fig. 6(d). The radius of the atomic model in the (111) plane (r Dis ) and the length of the atomic model in the [111] direction  (L Dis ) are 31.2 Å and 40.1 Å respectively, with a total of 13 920 Fe atoms. The radius of the free region (r free ) is 26.2 Å. First, to see individual H-atom diffusion around the screw dislocation core, a single H atom was placed at different positions in the screw dislocation core. Thereafter, MD simulations were performed at 600 K on an NVT ensemble for 130 ps. As the results, the following three typical diffusion motions of H atoms were observed [ Fig. 6(e)]: (1) The first type of H-atom motion, indicated with purple trajectory, exhibited a spiral motion along the Burgers vector of dislocation and was frequently pulled out from the core region to the energetically preferred sites (energy basin).
(2) The second type of H-atom motion, indicated with the dark gray trajectory, exhibited the same diffusion path as the purple one, but it eventually escaped the constraints of the energy basin and performed a random walk in the bulk region.
(3) The third type of H-atom motion, indicated with the yellow trajectory, exhibited that such random walking atoms in the bulk region are susceptible to being recaptured by the dislocation.
To obtain a more general diffusion path of an H atom in a system featuring dislocation, a simulation similar to the above-mentioned H charging of the GB was conducted. In this simulation, the atomic model sizes of r Dis , L Dis , and r free are 42.5, 40.1, and 32.5 Å, respectively, with a total of 20 301 Fe atoms. A vacuum region with a thickness of 10 Å was set in the [111] direction to enable the injection of H 2 gas. The injected H 2 molecules were constrained within a cylindrical region in the [111] direction centered at the core location of the dislocation, taking the radius as the free iron region. H-atom diffusion was enhanced by gradually increasing the quantity of the H 2 molecules in that vacuum region. The trajectories of the diffused atoms are displayed in Fig. 6(f) and Fig. S11 in the Supplemental Material [27]. Although a diffusion path inside the dislocation core along the Burgers vector is formed [dark red trajectory in Fig. 6(f) and in Figs. S11(a)-S11(c)], more atoms prefer diffusing through the bulk region via random walking ([green trajectory in Fig. 6(f) and Figs. S11(d)-S11(p)], and H-atom diffusion is faster in bulk than that in dislocation [ Fig. 6(f) and Fig. S11]. After carefully checking the snapshots in the charging process, we further uncovered that the screw dislocation traps random-walking H atoms and those H atoms can be transported along the screw dislocation line through repeated trapping and detrapping from the screw dislocation core. However, hydrogen transportation efficiency along a screw dislocation was really poor because of the higher solution energy in the dislocation core relative to the basin region, frequently causing trapped H atoms to be pulled out from the core region and kept in the energy basin.

B. Hydrogen trapping and desorption at defects
Understanding the hydrogen distribution in defects not only helps to characterize defects in metals [108], but can also guide the artificial creation of defects designed to trap and reduce mobile H atoms in metals. However, the hydrogen concentration in and around defects (e.g., vacancies, GBs, and screw dislocation) in α-iron at various temperatures remains an open question.
Combining the grand canonical Monte Carlo (GCMC) method [109] with MD simulations (the MD/GCMC hybrid method), we studied the radial distribution function of hydrogen volume density (ρ V H ) around various defects. Herein, to adjust the hydrogen bulk concentration in α-iron to approximately 40 atomic ppm (appm), the H-atom chemical potential was set to −2.16, −2.25, −2.37, −2.48, and −2.59 eV at 200, 300, 400, 500, and 600 K, respectively. A 12 × 12 × 12 supercell with one vacancy at the center, 3455 iron atoms in total, was employed as the vacancy model, and the models used in Sec. IV A were adopted for the GB and screw dislocation (the one with 13 920 Fe atoms) cases without inserting the vacuum region. 100 GCMC trials were conducted in each 10 MD steps with a time step of 0.5 fs.
To analyze the radial distribution function (RDF) of the volume density of the number of H atoms, spherical shell regions co-centered with the vacancy, slab regions with equivalent location above and below the GB interface, and cylindrical shell regions coaxed with the screw dislocation are set respectively for the models of vacancy, GB, and screw dislocation. The thickness of each region (i.e., the mesh size of the RDF analysis) was 2 Å. The ρ V H in an individual region was determined as follows: where the n H and V shell/slab denote the number of H atoms and the volume of an individual shell or slab region, respectively. These values were averaged over 5 × 10 6 GCMC trials in 50 000 MD steps. As shown in Fig. 7(a), the trapped H atoms were distributed within a very small region (within 2 Å) around the defects at all considered temperatures. An exception appeared in the case of screw dislocation [inner panel of Fig. 7(a)]; after carefully checking snapshots, this wider distribution was contributed by the H atoms located at E 3 and E 4 [ Fig. 5(a)]. Naturally the number of trapped H atoms decreased with increasing temperature in all defects. Snapshots of typical hydrogen distributions at various temperatures can be found in the Supplemental Material (Fig. S12) [27]. The temperaturedependent desorption rate, such as the quantity reduction of the trapped H atoms with respect to temperature, −∂ρ V H /∂T , is plotted in Fig. 7(b).
For the vacancy case, the ρ V H was 0.15 H/Å 3 at 200 K, equivalent to 5.01 H per vacancy. The ρ V H gradually decreased from 0.14 H/Å 3 to 0.06 H/Å 3 as the temperature increased from 300 to 500 K (equivalent to a desorption of three H atoms per vacancy). The desorption peak presented around 400 K in Fig. 7(b) matches the peak obtained at 400-440 K in experiment [110].
For the GB case, at 200 K, H atoms were trapped in regions A and C [ Fig. 6(a)]. As shown in Sec. III F, the solution energy of an H atom at s5 (corresponding to region C) is positive at 0 K, but is still 0.2 eV lower than at the T site in bulk. The desorption of these H atoms gives rise to the peak at 300-400 K in Fig. 7(b), which aligns with the experimental desorption peak at 400 K ascribed to the desorption of H atoms released from GBs [111]. For a screw dislocation model with a length of 1b, there are three energy basins, each lying between two neighboring core atoms [ Fig. 5(a)]. At 200 K, each energy basin can accommodate at most 1 H atom. Importantly, the desorption peak should exist below 300 K at least [see Fig. 7(b)]. This argument might explain the experimental desorption peak observed at 220 K [111].
Before closing this section, we would like to emphasize that any trapped H atoms can convert to mobile atoms at higher temperatures. This hydrogen state transition is critical for the understanding and prevention of HE.

C. MD-GCMC tensile test for H-charging bicrystal model with twist GB
Applying MD with EIP, Wan et al. recently proposed a fresh model for HE. They suggested that GB activation by the GB-dislocation reaction is a key component in the HE process of metals [106]. Nevertheless, the quantitative behavior must be validated by DFT calculations to ensure the accuracy of the result [112]. Here, MD-GCMC hybrid tensile tests for a twist grain boundary with and without hydrogen charging were performed to create the connects among hydrogen charging, GB activation, GB-dislocation reaction, and crack nucleation and propagation at DFT accuracy.
Before performing tensile tests, the twist GB formation energies (γ twist ) of 9{110} and 5{100}, with 288 and 60 Fe atoms, respectively, were calculated. The γ twist of the 9{110} twist GB obtained with NNIP is 0.95 J/m 2 and shows perfect agreement with the DFT result of 0.98 J/m 2 . The γ twist of the 5{100} twist GB with NNIP is 2.21 J/m 2 , which well agrees with the reported DFT result of 2.12 J/m 2 [52] and our DFT result of 2.29 J/m 2 . The performance of NNIP for twist GB shows its good transferability because no twist data were included in the reference database.
The 9{110} twist GB (hereafter referred to as 9 twist or twist) was adopted for the tensile test. The dimension of the H-free clean 9 twist is 51.2, 48.2, and 96.9 Å in [114], [221], and [110] directions, respectively; periodic boundary conditions were employed in all directions, and the model contained a total of 20 736 Fe atoms. The twist model is illustrated in Fig. 8(a).
Two comparison tensile tests were performed for the clean and H-charging twists. A constant MD tensile strain rate of 1 × 10 9 s −1 in the [110] direction and the temperature of 300 K were adopted for both tensile tests. The time steps of 1.0 and 0.5 fs were used in the tensile tests for the clean and Hcharging twists, respectively. For the clean twist, normal MD simulation on an NPT ensemble with constant tensile engineering strain rate was employed. For the H-charging twist, H atoms were introduced via the GCMC method with a constant H chemical potential of −2.25 eV (bulk H concentration is 40 appm at 300 K) during MD simulation on NVT(MD) and μVT(GCMC) ensembles. Therefore, H atoms can be inserted in and extracted from the supercell, and then the Fe-H system may approach an equilibrium state. 500 GCMC trials were conducted every 0.05 ps of MD simulation; finally, approximately 340 H atoms were introduced into the twist model within 100 ps of MD simulation. Then, the MD-GCMC hybrid tensile test on NPT and μVT ensembles was conducted, while the lateral dimensions of the model were adjusted according to the zero-pressure conditions. Note that time cannot be real in this simulation because GCMC trial steps were inserted during the MD simulation. However, for simplicity, here we measured time according to the number of MD steps, thus assuming that H concentration is always nearly equilibrium. The clean twist model exhibits elastic deformation when the tensile engineering strain is lower than 14.0%, as shown in the left panel of Fig. 8(b). Its stress starts to decrease after reaching the maximum of 21 GPa at a strain of 14.5%. The release of the stress was introduced by the emission of dislocation, as shown in the middle panel of Fig. 8(b); the dislocation volume density (ρ Dis , dislocation length divided by the volume of the model) sharply increased from 0 to 2.8 × 10 −3 Å −2 at that strain. Several snapshots in the tensile test are shown in Fig. 8(d).
The H-charging twist model was also subject to elastic deformation within the tensile stain of 11%. The initiation of plastic deformation of the H-charging twist was presented around the strain of 12.0% with a stress of 17 GPa, which is  [113] and dislocations, respectively, which are analyzed by the DXA tool implemented in the OVITO code [107]. lower than that of the clean twist. As presented in the snapshot (0.0%) in Fig. 8(e), even before applying tensile strain, the H concentration in the GB region was much higher than that in the bulk, i.e., H segregation. The atomic concentration of H atoms (ρ a H ) in the model increases with the tensile strain, as shown in the right panel of Fig. 8(b). This phenomenon can be also observed from the volume density of the number of H atoms (ρ V H ) in the slabs with the thickness of 5 Å around the GB at various tensile strains, as shown in Fig. 8(c). The solubility of hydrogen in the twist model increases with the tensile strain (stress) in both GB and bulk zones. A crack embryo is nucleated at the strain of 12.2% that was introduced by the emission of dislocation with an excitation of the GB structure, as shown in Figs. 8(b) and 8(e). The crack embryo was gradually expanded until it crossed the entire interface at the strain of 12.5%; i.e., the H-charging twist model was completely fractured. Then, 70% of the H atoms in the model were uncharged during the fracture process, as shown in the right panel in Fig. 8(b).
Defects of vacancy (Vac) and self-interstitial atom (SIA), created by the gliding of dislocation (Dis), were left in the separated grains, and the H atoms were mainly distributed in these defects and the newly created surfaces, as shown in the snapshot (17.0%) in Fig. 8(e). There is an interesting phenomenon shown in the snapshot (14.0%) in Fig. 8(e): the hydrogen concentration in the upper grain is significantly higher than that in the lower grain. The defect density (or the stress field introduced by the defects) might have a strong influence on hydrogen desorption.
Based on the above-mentioned two tensile tests, we observed that the presence of H atoms can facilitate dislocation emission from GB. The interaction between dislocation and GB, including dislocation emission from GB and impingement on GB, plays an important role in the process of crack embryo nucleation [106]. The actual fracture should be a result of several aspects; in addition to the interaction of dislocation and GB, the vacancy diffusion or clustering and nanovoid growth should also be involved. However, the effects of vacancy and nanovoids were not observed owing to the high strain rate adopted in this study.
Note that in actual materials, many preexisting dislocations exist. Thus, these dislocations can move at a much lower stress level than the stress level observed in the above simulation starting from the preexisting dislocation-free condition. Some of the moving dislocations may impinge on the GB and excite the GB structure at a much lower stress level than that in our simulation, eventually leading to the nucleation of the crack embryo. However, the stress level of the nucleation of the crack embryo should be at a similar level as that in the above simulation because the nucleation of the crack embryo requires a bond breaking. Moreover, the local stress level at the tip of a propagating crack must be at the same level as that in our simulation; thus, similar phenomena can be observed at the tip of the crack.

V. CONCLUSION
A general-purpose NNIP for the α-iron and hydrogen binary system was constructed based on the DFT reference database. This NNIP has been comprehensively tested for both trained and untrained data. The presented NNIP is reliable for the description of pure α-iron and the interactions of hydrogen with bulk and various defects in α-iron, including vacancies, surfaces, grain boundaries, and screw dislocation. Importantly, by virtue of the excellent flexibility of NNIP, this NNIP can solve several challenges facing Fe-H EAM potentials, such as the bond-breaking and bond-forming processes of H atoms at the iron surfaces and the energetics of screw dislocation motion.
The constructed NNIP was applied to several systems at finite temperature with the consideration of the HE process. The hydrogen charging and diffusion processes for 5(310)[001] tilt GB and 1 2 [111] screw dislocation can be directly observed at the atomic level: the tilt GB is a ditch for hydrogen diffusion, while the screw dislocation is a low-efficiency pipe for hydrogen transportation. H atoms are trapped and distributed within a small range around the defects of vacancy, tilt GB, and screw dislocation. The desorption temperature of the trapped H atoms agrees with experimental observations. In comparison in MD-GCMC hybrid tensile tests for a twist GB with and without hydrogen charging, the presence of hydrogen can facilitate dislocation emission from the GB, eventually leading to nucleation of the crack embryo at the GB region. The growth of eth crack embryo results in intergranular fracture. This proved that the GB-dislocation reaction is an important piece of the HE process. These applications showed that the constructed NNIP has a potential contribution to the understanding of the physics in the HE process.