Interatomic Fe–H potential for irradiation and embrittlement simulations

Interatomic Fe–H potential for irradiation and embrittlement simulations

Pekko Kuopanportti Erin Hayward Chu-Chun Fu Antti Kuronen Kai Nordlund School of Physics and Astronomy, Monash University, Victoria 3800, Australia Department of Physics, University of Helsinki, P.O. Box 43, FI-00014 Helsinki, Finland CEA, DEN, Service de Recherches de Métallurgie Physique, F-91191 Gif-sur-Yvette, France
January, 2016

The behavior of hydrogen in iron and iron alloys is of interest in many fields of physics and materials science. To enable large-scale molecular dynamics simulations of systems with Fe–H interactions, we develop, based on density-functional theory calculations, an interatomic Fe–H potential in the Tersoff–Brenner formalism. The obtained analytical potential is suitable for simulations of H in bulk Fe as well as for modeling small FeH molecules, and it can be directly combined with our previously constructed potential for the stainless steel Fe–Cr–C system. This will allow simulations of, e.g., hydrocarbon molecule chemistry on steel surfaces. In the current work, we apply the potential to simulating hydrogen-induced embrittlement in monocrystalline bulk Fe and in an Fe bicrystal with a grain boundary. In both cases, hydrogen is found to soften the material.

Interatomic potential, Molecular dynamics, Tensile testing, Hydrogen-induced embrittlement
61.43.Bn, 75.50.Bb, 61.72.Mm, 81.40.Jj
journal: and published in Comput. Mater. Sci.\biboptions


1 Introduction

Hydrogen, although not soluble in iron in equilibrium, can be introduced into it by irradiation, nuclear decay, or chemical processes. Hydrogen is well known to cause embrittlement in iron and steel Oriani (1978); Kimura and Kimura (1986); Bernstein (1970); Zhong et al. (2000); Nagumo et al. (2001); Tateyama and Ohno (2003); Geng et al. (2005); Tian et al. (2011), which is a serious issue in, e.g., the automotive and nuclear industries. In the former, the high mechanical resistance desired from the body steels must often be traded off against their increased susceptibility to hydrogen embrittlement Oriani et al. (1985); Lovicu et al. (2012); Koyama et al. (2012a, b), while the nuclear processes in the latter will, on long time scales, induce hydrogen buildup in the reactor steels Alter and Weber (1965); Hirth and Johnson (1976); Scott (1994); Greenwood (1994). Moreover, the recent changes in the design of the ITER fusion reactor are to render some of its steel components directly exposed to the fusion plasma Merola (), making it important to study how the energetic H isotopes escaping from the plasma interact with Fe.

Atomic-level molecular dynamics (MD) simulations have proven to be a good tool for examining irradiation effects Averback and Diaz de la Rubia (1998); Nordlund and Djurabekova (2014), mechanical properties of materials Keblinski et al. (1997); Schiøtz et al. (1999); Bringa et al. (2005), and plasma–wall interactions Salonen et al. (1999); Nordlund et al. (2014). The key physical input for MD is the interatomic potential. Since steels by definition contain Fe and C Callister, Jr. (1993), simulations of H effects in steels require, at a minimum, a potential that can describe all interactions in the ternary Fe–C–H system.

In this work, we develop a potential for Fe–H interactions in the same reactive Tersoff–Brenner formalism Tersoff (1988); Brenner (1990); Albe et al. (2002a) we used previously to construct a potential for the stainless steel Fe–Cr–C system Henriksson et al. (2013). The potential is fitted to a database of properties of FeH molecules and H in bulk Fe, obtained from literature and our own density functional theory (DFT) calculations. By using already available C–H parameters Brenner (1990); Juslin et al. (2005), the potentials can be directly combined to model the entire ternary Fe–C–H system. The potential allows simulating H in bulk Fe as well as ion irradiation and chemical reactivity of H at Fe and Fe–C surfaces. We demonstrate its use in simulating hydrogen-induced softening in bulk Fe and in Fe grain boundaries.

Fe–H potentials have been already devised using the embedded-atom method Lee and Jang (2007); Ramasubramaniam et al. (2009). However, its associated functional form cannot realistically describe the C–H bonding chemistry Brenner (1990). Since our aim is to obtain a potential for the entire Fe–C–H system, we choose to develop the Fe–H potential in the Tersoff–Brenner formalism which allows combining the Fe–H part with both Fe and C interactions, similar to what was done earlier for the W–C–H system Juslin et al. (2005).

The remainder of this article is organized as follows. In Section 2, we summarize the Tersoff–Brenner potential formalism and describe our fitting procedure. Section 3 presents the obtained Fe–H potential and evaluates its performance against experimental and ab initio data. In Section 4, we employ the potential in tensile-test simulations of hydrogen-containing iron. We discuss the implications and limitations of the study in Section 5. Finally, Section 6 concludes the article with a brief summary.

2 Potential formalism and fitting procedure

The reactive Tersoff–Brenner formalism Tersoff (1988); Brenner (1990); Albe et al. (2002a) used in this work originates from the concept of bond order proposed by Pauling Pauling (1960), and it has been shown Albe et al. (2002b) to resemble both the tight-binding scheme Cleri and Rosato (1993) and the embedded-atom method Daw and Baskes (1984); Brenner (1989). Since the formalism has been described extensively elsewhere Albe et al. (2002b, c); Nord et al. (2003), we will give here only a brief overview.

The total cohesive energy of the system is written as a sum of individual bond energies:


where is the distance between atoms and , is a cutoff function for the pair interaction, is a repulsive and an attractive pair potential, and is a bond-order term that describes three-body interactions and angularity. The pair potentials are of the Morse-like form


where and are the bond energy and length of the dimer molecule, respectively. The parameter is related to the ground-state vibrational frequency and the reduced mass of the dimer according to


The bond-order term is given by




Here is the angle between the vectors and , and the angular function is defined as


where , , , and are adjustable parameters. The range of the interaction is restricted to the next-neighbor sphere by the cutoff function


where and determine the locus and width of the cutoff interval.

If the analytical potential is used for modeling nonequilibrium phenomena involving short-distance interactions, such as high-energy particle irradiation processes or melting, the short-range part of the potential must be adjusted to include a strong repulsive core that follows, i.a., from the Coulomb repulsion between the positively charged nuclei. To this end, the potential is modified in the manner already used for other Tersoff-like many-body potentials Albe et al. (2002b); Nordlund et al. (1996): The total potential is constructed by joining the universal Ziegler–Biersack–Littmark potential  Ziegler et al. (1985) with the equilibirium potential using


where is the potential implied by Eq. (1) and is the Fermi function . The values of the parameters and are chosen manually such that the potential is essentially unmodified at equilibrium and longer bonding distances and that a smooth fit at short separations with no spurious minima is obtained for all realistic coordination numbers.

In order to devise a well-performing Fe–H potential in the Tersoff–Brenner formalism, we use the following fitting procedure: The parameter sets for the H–H and Fe–Fe interactions are taken unchanged from Refs. Brenner (1990) and Müller et al. (2007), respectively, so that only the parameter set for the Fe–H interactions is fitted. From the outset, we fix the parameters pertaining to the properties of the dimer FeH—i.e., , , and —according to their experimentally observed values. To avoid unwanted side effects, we set the three-index parameters to zero. The values of the remaining seven parameters (, , , , , , and ) are then fitted to a structural database comprising the molecules FeH and FeH, the stoichiometric FeH with the rock-salt crystal structure, and the energies of the lowest-lying hydrogen point defects in body-centered cubic (bcc) iron. The fitting is formulated as a nonlinear least-squares minimization problem, which we solve using the trust-region-reflective algorithm Moré and Sorensen (1983); Byrd et al. (1988); Branch et al. (1999) implemented in MATLAB mat ().

3 Obtained potential

The optimized parameter values for the analytical Fe–H potential are given in Table 1. We also show the parameter values used for the H–H and Fe–Fe potentials; it should be noted, however, that during the fitting process, H–H and Fe–Fe interactions play a role only in the evaluation of the hydrogen point-defect energies.

H–H Fe–Fe Fe–H
(eV) 4.7509 1.5 1.630
(Å) 0.7414 2.29 1.589
) 1.9436 1.4 1.875
2.3432 2.0693 4.000
12.33 0.01158 0.01332
0.0 1.2899 424.5
1.0 0.3413 7.282
(Å) 1.40 3.15 2.497
(Å) 0.30 0.2 0.1996
) 15.0 2.9 16.0
(Å) 0.35 0.95 1.0
Table 1: Tersoff–Brenner parameters [Eqs. (1)–(9)] for the Fe–H system. The H–H potential is taken from Ref. Brenner (1990) and the Fe–Fe potential from Ref. Müller et al. (2007); the Fe–H potential is derived in this work. The parameters are zero in all cases. In Section 4, we use  Å for the Fe–Fe interactions instead of  Å.

Table 2 presents the fitting database together with the results from the analytical potential. The potential exactly reproduces the experimentally observed dimer properties that are also in good agreement with all the ab initio calculations. For the linear trimer FeH, the experimentally measured bond lengths in Refs. Körsgen et al. (1999) and Körsgen et al. (1996) are, respectively, 1.1% and 2.1% greater than the analytical prediction, while differing from each other by roughly the same relative amount. The analytical potential yields a bond length for the trigonal planar FeH molecule that falls between the values given by the different ab initio calculations. The lattice constant , bulk modulus , and its pressure derivative for the rock-salt FeH are also in line with the DFT results.

Expt. Ab initio calculations AP
FeH CI/ECPSodupe et al. (1990) MRCP4Tanaka et al. (2001) MRSDCITanaka et al. (2001)
(Å) 1.589 1.578 1.596 1.582 1.589
() 1774 1701 1735 1778 1774
FeH linear CASSCFSiegbahn et al. (1984) CISiegbahn et al. (1984) B3LYPWang and Andrews (2008)
(Å) 1.648 1.746 1.689 1.645 1.630
FeH planar UHFBalabanov and Boggs (2000) MBPT4Balabanov and Boggs (2000) CCSDTBalabanov and Boggs (2000)
(Å) 1.667 1.603 1.609 1.619
FeH rock salt USPEXBazhanova et al. (2012) MBPPElsässer et al. (1998) FLAPWElsässer et al. (1998)
(Å) 1.828 1.839 1.833 1.839
(GPa) 270.8 216 200 238.9
4.25 3.7 3.7 4.749

Balfour et al. Balfour et al. (1983)
Schultz and Armentrout Schultz and Armentrout (1991)
Stevens et al. Stevens et al. (1983)
Körsgen et al. Körsgen et al. (1999)

Table 2: Properties of the Fe–H molecular and rock-salt phases as obtained from experiments, ab initio calculations, and the analytical potential (AP) derived in this work. The notation is as follows: , bond length; , wave number for the ground-state vibrational frequency; , cohesive energy per atom; , lattice constant; , bulk modulus; , pressure derivative of the bulk modulus. For the abbreviations of the ab initio methods, see Table 1 in the appendix.

Table 3 lists the formation energies of hydrogen point defects as obtained using DFT and our analytical potential. The formation energies are defined as


where is the total cohesive energy of the defect-containing cell with iron and hydrogen atoms; and are the atomic cohesive energies of bcc iron and the H molecule, respectively. For the analytical potentials of Table 1, and . A summary of the DFT methods is provided in the Appendix.

DFT calculations
Defect Expt. Hirth (1980) PAWJiang and Carter (2004) USPPMiwa and Fukumoto (2002) This w. AP
Unrelaxed T interst. 0.29 0.484 0.515
Relaxed T interst. 0.296 0.20 0.30 0.234 0.240
Unrelaxed O interst. 0.76 0.822 1.186
Relaxed O interst. 0.33 0.259 0.256
Unrelaxed substit. 2.855 4.027
Relaxed substit. 2.526 3.145
Table 3: Experimental and theoretical formation energies in units of eV for the most relevant hydrogen point defects in bcc iron: the tetrahedral (T) and octahedral (O) interstitials and the substitutional defect. Theoretical values are given for both unrelaxed and relaxed atomic configurations. The penultimate column lists our DFT results, and the last column shows the values calculated from our analytical potential (AP). For the DFT methods, see the appendix.

Due to the low solubility and high mobility of hydrogen in iron, as well as the high probability of trapping at defect sites at low temperatures, little direct evidence for the site occupancy exists. Indirect evidence indicates that H resides in the tetrahedral (T) site of bcc Fe, with an experimental value of  eV per atom for the dissolution energy of H in Fe Hirth (1980). According to the DFT results in Table 3, the T site is more stable, both for the unrelaxed and relaxed structures. The DFT calculations also indicate that the octahedral (O) site occupancy gains significantly more stabilization from lattice distortion than the T site does. This is because the O site undergoes a greater structural distortion than the T site, which can be understood heuristically by considering the sizes of the two sites: Using the lattice constant  Å, the radii of the T and O sites are 0.36 Å and 0.19 Å, respectively. The hydrogen atom has a covalent radius of 0.37 Å, so it fits better in the T site and causes smaller lattice distortions than in the O site. The same argument also explains why the energy of the substitutional defect decreases only slightly when relaxed. The analytical potential qualitatively reproduces this behavior and yields very good quantitative agreement for all three relaxed defect energies.

Regarding the diffusion of hydrogen in bcc iron, Jiang and Carter Jiang and Carter (2004) used DFT to obtain the Arrhenius equation for the diffusion coefficient, , where and the zero-point-energy-corrected activation energy corresponds to direct hopping between two neighboring T sites. Our DFT calculations yield for the same transition path Hayward and Fu (2013). On the other hand, since H is easily trapped by impurities in Fe, the diffusion coefficients of H in Fe from laboratory measurements show a large scatter: Hayashi and Shu Hayashi and Shu (2000) compile experimental values of in the range from  eV to  eV. Using our analytical potential and the nudged elastic band method Mills et al. (1995), we get for the nearest-neighbor path. This value and the two DFT results all fall within the experimental range.

4 Effect of hydrogen on tensile testing of iron

As a demonstration of possible applications of the derived Fe–H potential (Table 1), we employ it in MD simulations Nordlund () to investigate the effect of hydrogen impurity atoms on the stress–strain response of crystalline iron subjected to uniaxial tensile stress. We consider two types of computational cells of iron atoms, one consisting of a regular bcc lattice of unit cells and one containing a bcc bicrystal with a grain-boundary plane at its center; the latter is illustrated Stukowski (2010) in Fig. 1. Both cells have  Å, and periodic boundary conditions are imposed in all three directions. The axis and angle of rotation for the grain boundary are chosen as and , so that due to the periodic boundary conditions, the structure corresponds to a stack of symmetric tilt boundaries with a separation distance of  Å and a grain-boundary energy of , where is the total energy of the computational cell and  Å is its cross-sectional area perpendicular to the direction. In crystallographic notation Lejček (2010), the grain boundary can be described as .

Figure 1: (a) Perspective and (b) side views of the computational cell used for the grain-boundary system with an atomic hydrogen concentration of 2.0%. Iron atoms () are shown in gray (light) and hydrogen atoms () in blue (dark). The dimensions of the cell are , and periodic boundary conditions are imposed in all three directions. The cell contains symmetric tilt boundaries in its top/bottom and middle sections. (b) Neighboring grains are tilted about the axis by the angle with respect to each other; the vectors and correspond to equivalent lattice directions in the adjacent grains. The vector lies along the grain boundary.

To introduce the impurities into system, we randomly place a variable number of hydrogen atoms into the computational cell, subject to the condition that the added atoms are at a minimum distance of 1.55 Å from the already existing atoms. For , we use the values 87, 174, 260, 346, 432, 519, 605, 691, 778, and 864, which corresponds to atomic hydrogen concentrations [defined as ] of 0.0%, 1.0%, 2.0%, 2.9%, 3.9%, 4.8%, 5.7%, 6.5%, 7.4%, 8.3%, and 9.1%, respectively. We let the stresses in the hydrogen-containing cell relax to zero by evolving the system for 20 ps at 300 K while applying the Berendsen pressure control Berendsen et al. (1984) in all directions.

Next, exertion of uniaxial tension in the direction is modeled in a quasistatic, stepwise manner: First, we increase the length of the simulation box by 0.02 Å, scaling the coordinates of all atoms by the ratio of the new and previous . Second, we evolve the system for 50 ps with fixed , while applying the Berendsen pressure control in the and directions and the Berendsen temperature control Berendsen et al. (1984) at 300 K, and extract the axial normal stress as a time average over the last 25 ps. These two steps are repeated 500 times, resulting in a maximum strain of 11–13%. We carry out the whole procedure for different values of and average the results for each over ten independent initial configurations of H atoms.

Since the tensile-test simulations are carried at 300 K, we have increased the value of the cutoff parameter for the Fe–Fe potential Müller et al. (2007) to 3.5 Å from the original 3.15 Å. Otherwise, the second-nearest neighbors of the bcc iron would—due to thermal vibrations—experience the onset of the cutoff function [Eq. (8)], resulting in an unphysical increase in the Young’s modulus of elasticity , where denotes the normal tensile strain. By extending to the middle of the second- and third-nearest-neighbor distances, this effect is avoided.

Figure 2(a) shows the average stress–strain curves for the regular monocrystalline (bulk) bcc iron, for six different atomic hydrogen concentrations. Figure 2(b) depicts the corresponding curves for the grain-boundary system, and Fig. 3 combines data from the two configurations.

Figure 2: Uniaxial nominal stress as a function of normal tensile strain for (a) monocrystalline bulk iron and (b) an iron bicrystal with a symmetric tilt boundary (Fig. 1). The legend shows the atomic hydrogen concentration for each curve. Uniaxial tension is exerted along the axis, perpendicular to the grain-boundary plane.
Figure 3: Comparison of the uniaxial stress–strain responses of bulk [see Fig. 2(a)] and grain-boundary [Fig. 2(b)] iron crystals for atomic hydrogen concentrations of 0% and 4.8%.

To quantify the stress–strain response in the linear, elastic regime, we determine the Young’s modulus at different values of the hydrogen concentration . This is done by performing a linear least-squares fit to each of the stress–strain curves in Fig. 2. In the case of the bulk system, we use the fitting intervals for , for , and for ; for the grain-boundary system, they are for and for . The reason for not always starting the fitting interval from zero strain is that due to the high mobility of hydrogen in iron, there was some reorganization of the hydrogen atoms during the first few steps of the stretching procedure, producing nonlinear stress-strain behavior. As can be seen from Fig. 2(a), the nonlinearity is particularly pronounced for bulk samples with large .

Figure 4(a) shows the obtained Young’s moduli at 300 K for  not (a). The error bars are calculated as the standard deviation of each set of ten simulations. Considering first pure iron at 0 K, its Young’s moduli can be determined directly from the elastic moduli predicted by the potential Müller et al. (2007); Nye (1985). This gives the values 115 GPa for the bulk bcc and 164 GPa for a bcc system in which the lattice is rotated by the same angle and in the same direction as in our grain-boundary system. When the temperature is increased, the Young’s moduli are expected to decrease. For pure iron at 300 K, our simulations yield 101 GPa for the bulk and 141 GPa for the grain-boundary system. The experimental value Adams et al. (2006) for the bulk system at 300 K is 132 GPa not (b).

For both configurations, hydrogen is observed to induce softening of the material, i.e., to reduce . The effect is noticeably stronger for the bulk system: when the hydrogen concentration increases from zero to 9.1%, the Young’s modulus of the bulk system decreases by 55%, whereas for the grain-boundary system the decrease is only 22%. One possible explanation for this difference is that most of the hydrogen atoms in the grain-boundary system resided within  Å from one of the boundaries. Therefore, the concentration of hydrogen in the intact lattice was significantly lower than the nominal hydrogen concentration (e.g., at it was less than 4%), while in the bulk system, these two quantities were obviously equal and the H atoms were homogeneously distributed throughout the computational cell.

Let us next investigate the extreme plastic behavior in terms of the tensile strength , defined as the maximum stress reached by the stress-strain curve. The resulting values for the bulk and grain-boundary crystals are presented as a function of the atomic hydrogen concentration in Fig. 4(b). From there, we see that the addition of hydrogen decreases the tensile strength of both configurations. This can be understood by noting that the H atoms introduce disorder in the Fe lattice. As in the case of the Young’s modulus, the decrease is more substantial in the bulk than in the grain-boundary system. Without hydrogen, their tensile strengths are approximately equal ( GPa for the bulk and 7.27 GPa for the grain-boundary system). At a hydrogen concentration of 9.1%, however, the tensile strength of the bulk system has decreased by 54%, while in the grain-boundary system the decrease is only 29%. The reason for the weaker effect in the grain-boundary system is likely the same as mentioned above for the Young’s modulus.

Introduction of hydrogen into the system also modifies the shape of the stress–strain curves near the maximum stress. For pure bulk iron, there is a sudden drop in the stress at 9% strain [Fig. 2(a)]. Visual inspection Stukowski (2010) of the simulation system reveals that this is caused by a slip process that creates stacking-fault ribbons extending through the system in the direction. A similar but less distinct drop occurs for the grain-boundary system [Fig. 2(b)]. In this case, the stacking fault cannot extend through the whole system because the grain boundary interrupts the crystal structure. The presence of hydrogen smooths the abrupt drops by significantly disturbing the crystal lattice already before the stacking-fault ribbons appear.

Figure 4: (a) Young’s modulus and (b) the tensile strength as functions of the atomic hydrogen concentration for bulk bcc iron (triangles) and for the grain-boundary system (circles). The error bars are calculated as the standard deviation of ten simulations. The lines are guides to the eye.

5 Discussion

We have developed an analytical Tersoff–Brenner potential for interactions between hydrogen and iron atoms. It was fitted to a set of experimental and ab initio data on iron hydride molecules, rock-salt-structured crystalline FeH, and hydrogen point defects in iron. The obtained potential reproduces the experimentally measured bond energy, bond length, and ground-state vibrational frequency of the FeH dimer and describes with good accuracy the molecules FeH and FeH as well as the rock-salt FeH. The point-defect energies it predicts are also consistent with our own DFT calculations.

The constructed potential enables atomistic computer simulations of a wide range of materials problems involving iron and hydrogen. Since it can also model nonequilibrium phenomena such as sputtering and the formation of mixed materials, the potential is well-suited for MD studies of plasma–wall interactions in fusion reactors. In view of the recent design updates of the ITER reactor, which would result in direct exposure of steel to the fusion plasma Merola (), being able to incorporate both iron and hydrogen into these investigations is an important advancement. With the Fe–H part now available, the only Tersoff–Brenner potential missing from the quaternary Fe–Cr–C–H system is Cr–H, whose development we leave for future work.

In Section 4, we applied the potential to tensile-test simulations of iron in two different configurations, a bulk bcc monocrystal and a symmetric tilt boundary, using different concentrations of hydrogen impurity atoms. The simulations indicated that hydrogen softens iron; i.e., the Young’s modulus and the tensile strength decrease when the hydrogen concentration increases. The effect was much stronger in the bulk bcc monocrystal than in the tilt-boundary system. This was explained by noting that most of the hydrogen in the grain-boundary system was concentrated near the grain boundaries, thereby leaving the rest of the system depleted in hydrogen in comparison to the bulk system, where hydrogen was homogeneously distributed.

Our simulations demonstrate that the potential can be used to study hydrogen-induced embrittlement phenomena in iron and steel. We emphasize that the current simulation setup is constructed as a simple model system for potential testing. In likely experimental scenarios, most H would (due to its low solubility) be trapped in defects or in grain boundaries. Thus, the potential’s prediction of hydrogen-induced grain-boundary weakening is at least qualitatively consistent with the well-known effect of grain-boundary embrittlement by H in steels Bernstein (1970); Zhong et al. (2000); Nagumo et al. (2001); Tateyama and Ohno (2003); Geng et al. (2005); Tian et al. (2011). Future work could examine this more systematically for other grain boundaries and hydrogen distributions.

6 Summary

We constructed a DFT-based interatomic potential for the Fe–H system in the Tersoff–Brenner formalism. The potential can be directly combined with our previously developed potential for the steel Fe–C system, to enable simulations of hydrocarbon chemistry in steel bulk and at steel surfaces. We applied the new potential to investigating the effect of hydrogen on the mechanical properties of monocrystalline bulk Fe and an Fe bicrystal with a grain boundary. In both cases, hydrogen was found to soften the material, reducing the Young’s modulus as well as the tensile strength.


We thank T. Ahlgren, C. Björkas, F. Granberg, K. O. E. Henriksson, A. Lasa, and M. Nagel for insightful discussions. We are grateful for the computational resources granted by the CSC – IT Center for Science in Espoo, Finland. P. Kuopanportti acknowledges financial support from the Emil Aaltonen Foundation, the Finnish Cultural Foundation, the Magnus Ehrnrooth Foundation, and the Technology Industries of Finland Centennial Foundation. E. Hayward and C.-C. Fu thank the ANR-HSynThEx Project and GENCI (x2015096020). This project has been carried out within the framework of the EUROfusion Consortium and has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement number 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Appendix A Ab initio methods

Here we outline our DFT calculations of the hydrogen point-defect energies in Table 3. They were performed with the SIESTA code Soler et al. (2002), were spin-polarized within the collinear approximation, and used the generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof exchange-correlation functional Perdew et al. (1996). Core electrons were replaced by nonlocal norm-conserving pseudopotentials, and valence electrons were described by linear combinations of numerical pseudoatomic orbitals. We represented the charge density on a real-space grid with a spacing of  Å and employed a Methfessel–Paxton smearing Methfessel and Paxton (1989) of 0.3 eV. All calculations used a 128-atom supercell with a  -point grid. Zero-point energy corrections, calculated for hydrogen within the Einstein approximation, are included in our quoted values.

In Tables 2 and 3, we also employed a number of abbreviations for the ab initio methods used in other studies. The abbreviations are defined in Table 1.

Method Description
B3LYP Becke’s three-parameter hybrid functional Becke (1993) with
Lee–Yang–Parr correlation functional Lee et al. (1988)
CASSCF Complete active space self-consistent field
CCSDT Coupled-cluster method in singles and doubles approximation
with connected triple excitation terms
CI Configuration interaction
DFT Density functional theory
ECP Energy-adjusted quasi-relativistic effective core potential
FLAPW Full-potential linear-augmented-plane-wave method with
von Barth–Hedin formula von Barth and Hedin (1972)
GGA Generalized gradient approximation for DFT
MBPP Mixed-basis pseudopotentials with Perdew–Zunger
parametrization Perdew and Zunger (1981) of Ceperley–Alder data Ceperley and Alder (1980)
MRCP4 Fourth-order approximation of coupled pair approach
MRSDCI Single and double excitation multireference CI with Davidson
quadruple excitation correction
PAW Projector-augmented-wave method within GGA
MBPT4 Fourth-order many-body perturbation theory restricted to
single, double, and quadruple excitations
UHF Unrestricted Hartree–Fock method
USPEX Evolutionary crystal structure prediction method Oganov and Glass (2006) (GGA)
USPP Ultrasoft pseudopotentials within GGA
Table 1: Abbreviations used for the ab initio methods in Tables 2 and 3 (in alphabetical order).


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description