Modeling and stabilities of Mg/MgH₂ interfaces: A first-principles investigation

Magnesium (Mg) is a promising candidate^1,2 for application of mobile hydrogen storage material due to the high gravimetric hydrogen content³ (7.6% in MgH₂), high volumetric density³ (6.7 × 10¹⁹ H atoms/L), and high energy density⁴ (9MJ/kg) of its hydride MgH₂. To make Mg applicable, there are still two technical obstacles^1,2 to be solved: (1) Slow hydrogenation/dehydrogenation kinetics; (2) A too high working temperature of 350 - 400°C. In order to improve the material hydrogen sorption performance, extensive theoretical and experimental researches including strain engineering,^5,6 material surface/bulk ratio increasing,^7–9 transition metal catalyzing^10–12 and mixing^13,14 with other compounds have been done to mitigate the thermodynamic and kinetic drawbacks of Mg.

Hydrogen sorption in Mg concerns to the Mg/MgH₂ or other related interfaces. Recently, some experimental and theoretical efforts have been made to study the Mg/MgH₂ interface structure and improve MgH₂ dehydrogenation properties. Epitaxial growth of MgH₂(110) on Mg(0001) is experimentally observed.¹⁵ The significant increase of equilibrium hydrogen pressure in Pd/Mg systems was first believed to be caused by the clamping effect¹⁶ but later found to be induced by the interfacial alloy hydride.¹⁷ Mooij et al.⁹ have found that the nanoconfined MgH₂ is more easily destabilized compared to its bulk counterpart because of an interface energy effect. Moreover, decreased dehydrogenation enthalpies could be also observed in Mg/Mg₂Ni interface system¹⁸ or on MgH₂ clusters⁷ with smaller size. Noticeably, the MgH₂ dehydrogenation mechanism was found¹⁹ not to be simply a reverse of the adsorption one. Although the barrier for Mg nucleation is small¹⁹ in dehydrogenation process, the interfacial void formations^19,20 accompany the Mg nucleation thus leading to a difficult determination of Mg/MgH₂ interface structures. Theoretically, it is found²¹ that the dehydrogenation facilitation of (100) MgH₂/TiH₂ multilayers attributes to structural deformations and interfacial energy variations. Moreover, by means of Car-Parrinello molecular dynamics method (CPMD), it is found that the atomic hydrogen diffusivity at MgH₂(110)/Mg(010) interface steeply increase²² at 700 K, and is facilitated²³ when Fe is doped. Though these experimental and theoretical studies provide interesting and illuminating knowledge of Mg/MgH₂ interface properties, there is still a lack of systematic study of relative stabilities of various Mg/MgH₂ interfaces, i.e. various combinations of different Mg and MgH₂ surfaces. Additionally, the theoretical studies have not fully taken into account the influences of relative position,²² interlayer spacing^22,24,25 and lattice mismatch^24,26–29 at the first-principles level, which may be determining the Mg/MgH₂ interface stabilities.

In this paper, we firstly investigate the matching between Mg and MgH₂ surfaces by considering their surface energies, lattice constants and cell shapes. Then we consider the interfacial elastic effect by determining the mutual lattice constant which leads to minimized energy cost on elastic effect. Furthermore, we indicate that the most stable interfacial atomic structures would be found according to their lowest interfacial energies and the interfacial geometric characters would be discussed in detail. Lastly, the interfacial properties and strain effects on defect formations are also provided. The methods are specified in Section II, the results and discussions are in Section III, while conclusions are in Section IV.

The incoherent interface between two phases, namely magnesium and magnesium hydride, is chosen to be modeled in this paper. The two phases are modeled to be thick enough so that they could be seen as the atom reservoirs, so as to stimulate the hydrogen sorption process. During the establishment of interface model, we do not consider the diffusion of H atoms at the interface of Mg/MgH₂, though it does not reflect the practical scenario as the hydrogenation/dehydrogenation often occurs above 600 K for Mg-H system. The trend of the H diffusion at the interface will be discussed from the formation energies of interstitial and vacancy H defects around the interface.

We start with the investigation of surface energies. Afterwards, some combinations of two phases could be found. Of these the chosen are determined according to the proposed rules, where the relative positions of two phases and the biaxial strains along each axis are considered. Finally, the comparison between the stabilities of interface models is made to discuss what the factors that dominate the interface stabilities are.

Our calculations are performed using Vienna ab initio Simulation Package.^30–32 We use the exchange-correlation with generalized-gradient approximation (GGA)³³ and projector-augmented wave (PAW) potentials.³⁴ In all the calculations, adequate tests with k-points meshes³⁵ have been done to show excellent converged results and the plane-wave cutoff energy was set to be 310 eV. Monkhorst-Pack³⁵ k-point samplings of 6 × 2 × 1, 5 × 2 × 1 and 5 × 2 × 1 are used in the Mg(⁠|$10\bar 10$|$101¯0$⁠)/MgH₂(210), Mg(0001)/MgH₂(101) and Mg(⁠|$10\bar 10$|$101¯0$⁠)/MgH₂(101) interfaces, respectively, for the Brillouin-zone integration. The internal optimization is converged with the force on each atom to be less than 0.01 eV/ Å. The energy convergence with respect to the k-points and cutoff energy has also been tested to be in the precision of ∼0.01 eV per cell.

Surface energies are seen as the energy cost of creating a surface by dividing a crystal into two parts. Here, the surface energies of magnesium (Mg) and magnesium hydride (MgH₂) are calculated by the following formula

where E^tot(n) and E_bulk represent the total energies of the slab model and the bulk with n formula units, A is the area of the surface slab and the factor 2 denotes two surfaces being created. Of note, the Fiorentini approach³⁶ is adopted to get a reasonable convergence of surface energies.³⁷ 

The interfacial energy E_i at the grain boundary between two phases is defined as the formation energy of the interface system in a particular context, where the chemical potentials at which two phases are at equilibrium could be determined. The calculated formula of interfacial energy proposed by Arras et al.³⁸ is given by

where E_tot is the total energy of interface model, N_Mg and N_H are the total number of atoms of magnesium and hydrogen respectively in the interface system, |$\mu _{{\rm Mg}}^0$|$μ Mg 0$ and |$\mu _{\rm H}^0$|$μH0$ are the chemical potentials of magnesium and hydrogen when two phases, namely Mg and MgH₂, are at equilibrium. In our model in which the two phases represent the atom reservoirs,³⁹ |$\mu _{{\rm Mg}}^0$|$μ Mg 0$ and |$\mu _{\rm H}^0$|$μH0$ are calculated by

and

where ɛ is the total energy per formula unit and x is the atomic concentration of an element (subscript) in a phase (superscript). Specifically, |$x_{{\rm Mg}}^{{\rm Mg}}$|$x Mg Mg$⁠, |$x_{\rm H}^{{\rm Mg}}$|$xH Mg$⁠, |$x_{{\rm Mg}}^{{\rm MgH}_{\rm 2} }$|$x Mg MgH 2$ and |$x_{\rm H}^{{\rm MgH}_{\rm 2} }$|$xH MgH 2$ are equal to 1, 0, 1/3 and 2/3, respectively.

The H defect formation energy E_f is defined as the formation energy of the H vacancy in the MgH₂ phase, or the H interstitial in Mg phase of the interface model. The calculated formula of the defect formation energy is given by

where E_tot and E_perfect are the total energies of interface models with and without H defect, N_defect indicates the number of H defects, and the signs + and - represent H vacancy and interstitial, respectively. Here, in order to stimulate the H-rich context, μ_H is equal to 1/2E_H2.

The surface energies strongly reflect the possibilities of formation of a particular surface. It should be noted that our first-principles calculation of surface energies is conducted at 0K, where the bond cutting and the surface reconstruction are the two major factors that dictate the surface energies. Previous papers^37,40,41 showed that both of these two effects played an important role on surface energies. Moreover, the contribution of vibration entropy, another factor in determining the surface energies, goes up when the temperature rises so that the two major factors may not be the only determining ones and the surfaces with high surface energies at 0 K may appear.⁴² As a matter of fact, the hydrogen desorption in magnesium hydride requires a temperature much higher than the room temperature. Therefore, the surfaces with high surface energies at 0 K also have to be considered in dictating the combination of two surfaces of two phases.

The periodic interface model which consists of two phases requires that the two phases share the same shape and close lattice constants. The perfectness of matching leads incoherent interfaces to coherent or semi-coherent interfaces.

The matching of shape takes priority over the matching of lattice constants. Specifically, the cells of rectangle shape (NB: the hexagonal shape can be transformed into rectangle shape) could always be matched easily after some unit repeating operations, which is in pursuit of existence of a coincidence site lattices^43,44 (CSL). Comparatively, the cells of parallelogram with unequal neighboring sides and a non-right angle are difficult to be matched with the others. For instance, it is found that the Mg(⁠|$10\bar 11$|$101¯1$⁠) surface cell is in its difficulty to be matched due to its shape of parallelogram where |$a_{{\rm Mg}(10\bar 11)} = 6.11\,{\rm {\AA}}$|$a Mg (101¯1)=6.11Å$⁠, |$b_{{\rm Mg}} (10\bar 11) = 6.38\,{\rm {\AA}}$|$b Mg (101¯1)=6.38Å$ and the included angle γ = 74.87°. The cell of this shape, unlike the rhombus, is also impossible to be transformed into the periodic cell of rectangle shape due to its specified atomic arrangement.

The unit repeating operations follow two geometrical criteria:^45,46 (1) The lattice misfit strain between two phases should be as small as possible. (2) The size of repeated cells, manifested by the period of the CSL, is as small as possible. The first criterion is easily understandable. The second one is based on the observation^43,45 that the grain boundaries between two identical materials are often found to have a CSL with small unit cell volume. The information about the size of our repeated cells is specified in Table II.

In summary, the rectangle cell (original or transformed) is always preferred in the process of choosing the surface of each phase. In practice, we choose the Mg(0001), Mg(⁠|$10\bar 10$|$101¯0$⁠), MgH₂(101) and MgH₂(210) surfaces whose information is specifically listed in Table I, which will be further discussed in following section. Fig. 1 provides an intuitive illustration of interface modeling and interface pairs in which we are interested.

For an interface system with lattice mismatch, both of the two phases would be possibly under elastically tensile or compressive effect. Therefore, the elastic cost subject to the mutual lattice constants is required when the interface is created. Here, the mutual constants are defined as the lattice constants of the interface model where two phases coexist with minimized elastic energy cost. The minimum of this elastic energy cost is expected to be found in our interface models. Basing on the assumption that both of two phases would compromise energetically to reach the mutual constant, we propose that the energy cost on elastic effect is calculated by

where E_EEC represents the total elastic energy cost of two phases, ΔE_Phase.1 and ΔE_Phase.2 denote the elastic energy costs of two phases with respect to the most energetically favored structure.

The search of the mutual lattice constants is required to determine the mutual constants. In practice, the calculations of ΔE_Phase.1 and ΔE_Phase.2 are calculated respectively in bulk forms of two phases by fitting the E-c curve at given a and b. Thus, the mutual constants could be determined by the minimized elastic energy cost. It is in line with widely adopted theoretical models, although it is relatively simple to the real interface where the lattice constants of both phases will be restored beyond corresponding critical thickness.⁴⁷ 

The stability of interface, after the determination of the mutual lattice constants, is highly influenced by the bonding nature in interfacial region, namely the interfacial structure. Naturally, the factors that lead to a change of interfacial structure should be focused on. One may easily come up with several factors including the relative position of two phases, the concentration gradient between two phases and the effect of defect or dopants in interfacial region.

The authentic hydrogen sorption process always exhibits migration of grain boundaries, which is affected not only by the thermodynamic stability but also by the kinetic one. From the point of view of theoretical calculation, it is not easy to describe all these factors simultaneously. However, when it comes to the discussion of the interfacial bonding, we simply assume that the relative position of two phases should be the primary factor in affecting the interfacial structure in our calculations since the variation of the relative position would lead to a dramatic change of interfacial bonding, whose effect is similar as the sliding shear acted on the system.

The crystal structure of magnesium is hexagonal close packing and the calculated Mg bulk lattice constants are a_Mg = b_Mg = 3.190Å with the c/a value of 1.623, which are in good agreement with the experimental values^48,49 and previous theoretical results.^50–52 

Magnesium hydride, whose crystal structure is simple tetragonal, is a highly stable hydride with strong Mg-H bond, which is a mixture of ionic and covalent bonding.⁵³ The calculated bulk lattice constants are a_MgH2 = b_MgH2 = 4.452 Å, c_MgH2 = 2.992 Å with the c/a value of 0.672, which agree well with the previous experimental^54–56 and theoretical^57–59 results.

Here, surface energies of two phases are calculated and compared to find the relative stabilities of surfaces. It is known that the low-index surfaces are often more stable than the high-index ones since the number of breaking bones⁶⁰ that are required to create a low-index surface is relatively less. However, it has been reported³⁷ that the low-index Mg(⁠|$10\bar 10$|$101¯0$⁠) has a high surface energy thanks to its massive basal bond cuttings. MgH₂(210) is of the highest surface energy in Table I, while the other MgH₂ surfaces share low surface energies.

In a theoretical calculation, the adopted model should be periodic. Consequently, the shapes and lattice constants of two phases are identical, which are widely accepted in previous literatures^{24,26,38,61–63} concerning interface. It means that the shapes and the lattice constants should be taken into account in combining the two surfaces of two phases.

Considering these factors mentioned in the above two paragraphs, i.e. the surface energies, the shape and lattice constants of a surface-slab cell, we choose the combinations, as listed in Table I, of which two surfaces of two phases are identical in shape and close in lattice constants (repeated along the axes if necessary). In practice, the low-energy Mg(0001) and MgH₂(101) are chosen while the high-energy Mg(⁠|$10\bar 10$|$101¯0$⁠) and MgH₂(210) are selected to exemplify our approaches.

Formation of interfaces always exhibits energy compromise of two phases, which can be reflected by elastic effect. We define here, for an ideal condition, the mutual constant as a constant which leads to minimized energy cost on elastic effect. Therefore, the search of mutual constants for interface models could be interpreted as searching the minimum E_EEC among a wide range of lattice constants.

In practice, the search of mutual lattice constant for interface could be divided into two steps: (1) The elastic costs on each phase are estimated. (2) Take into account the whole contribution from two phases caused by elastic effect. It should be pointed out here that the interface-induced strain effect on both phases may lead to different strains along axes a and b. Therefore the consideration of nonequivalent strain on axes a and b makes sense. Since the total energies of Mg and MgH₂ per formula unit are −1.56eV and −8.91eV respectively, the contributions from two phases, caused by elastic effect, also differ. It is clear that the thickness of the two phases significantly affects the minimization of the elastic energy in our model. Practically, the critical thickness, including the strain release effect of the layers within the critical thickness, should be included in the minimization of the elastic energy. Considering the difference in contribution of two phases, the energies of ΔE_Mg and |$\Delta E_{{\rm MgH}_{\rm 2} }$|$ΔE MgH 2$ are calculated with respect to the most stable structure within the searching range. Moreover, the number of Mg layers is nearly twice as that of MgH₂ layers so that the contributions from two phases could be competing. Here, we assume that our choice of thicknesses, made by the two above mentioned strategies, is capable of investigating the elastic effect, i.e. determining the mutual constants. Similar ratio between numbers of two phases, i.e. Ge and GeO₂, is used in previous literature.⁶⁴ 

Fig. 2 shows the case of Mg(⁠|$10\bar 10$|$101¯0$⁠)/MgH₂(210) interface. We apply various structures, which are distinguishable from their lattice constants, to two phases in their bulk form cell. The elastic effects are found to be obviously influencing the phase stabilities. The relative stabilities of varying lattice constants could be reflected by the energy difference of the respective structures with respect to the most energetically favored one (as shown in Figs. 2(a) and 2(b)). Considering the Possion effect, we fixed a and b but varied c, and then fitted the E-c curve to determine the energy for each structure. Afterwards, the whole elastic contribution from two phases are estimated by simply calculating E_EEC = ΔE_Phase.1 + ΔE_Phase.2 [see Fig. 2(c)]. The most energetically favored structure could be found among many structures that we considered. In this Mg(⁠|$10\bar 10$|$101¯0$⁠)/MgH₂(210) interface case, the mutual constants are found to be a = 3 Å and b = 10.2 Å. Here, we have made a cell repeating operation for Mg(⁠|$10\bar 10$|$101¯0$⁠) in order to match the interface (see Table II).

Similarly, the mutual constants searching proceeds on the cases of Mg(0001)/MgH₂(101) and Mg(⁠|$10\bar 10$|$101¯0$⁠)/MgH₂(101) interfaces as shown in Fig. 3. The details of three interface combinations are specifically listed in Table II. With the obtained lattice constants, the interfaces could be modeled.

The determination of mutual constants mainly focuses on elastic effect but does not involve the interfacial bonding, which primarily affects the interface stabilities. Recent papers have shown that interface properties highly depend on the lattice matching⁶⁵ or interface chemistry.⁶⁶ At Mg/MgH₂ interfaces, the Mg-H bonding is the most important bonding that leads to interface stability.

Authentically, there are interfaces consisted of two phases with distinctive symmetries in many cases, especially for the cases with high-index surfaces. For example, correlation between interface trap densities and interface anisotropy are found⁶⁷ on the interfaces consisted of high-index surfaces. Another case⁶⁸ is that the interface of Ag/high-index substrate enhances resonance by enhancing the scattering cross-sections. From a point of view of theoretical calculation, it is rather difficult to find the most energetically-favored structure for interfaces where two phases are with different symmetries. Therefore, the full consideration of relative position is essential.

We propose that the relative position of two phases could be seen as the vital reason causing change of interfacial bonding in theoretical calculations. For theoretical modeling, considerable interfacial structures could be obtained by adjusting the relative position of two phases. Therefore, it is possible to compare the stabilities of relaxed interfaces with various structures.

We have calculated the effect of relative position in Fig. 4. It could be obviously observed the varying of relative position could cause a dramatic change in interfacial energies (the largest difference is 73.2 meV/Ų (1.171 J/m²), up to 22.73% of the maximum). In every structure (at particular x and y), an energy dependence on d could be seen. Here, d is defined as the distance between two phases of initial configuration along z axis, as shown in the inset of Fig. 4(c). This is easy to understand since the increasing of d lengthens the distance of two phases and thus greatly alters the Mg-H bond lengths at interface. It would lead to the dramatic change of Mg-H bond energies and thus influence the interfacial stabilities. We find there are local energy minimums in every structure within the proposed d range. It means that there is a preference for d, which dovetails with strengthening of interfacial bonds, in interface models. Logically, it could be interpreted that the optimization of d is always required for interface models.

Moreover, it is found the varying of relative position on x-y plane (at fixed d) also induce a change of interfacial energies although not as dramatic as that induced by d varying. Specifically, the varying of x-y relative position would cause a change of local structure. But this effect would be obviously offset by relaxations of interfacial atoms and thus just induce a slight effect on the change of interfacial energies. Comparably, the varying of d would be more effective on tuning interfacial bonding. Higher effectiveness is mainly attributed to the fact that all interfacial bonds have to respond as d varies. Compared with x-y relative position varying, it is more global, which relaxations of interfacial atoms are not easy to offset.

Similar findings could be observed at interfaces of Mg(0001)/MgH₂(101) and Mg(⁠|$10\bar 10$|$101¯0$⁠)/MgH₂(101) interfaces (see Fig. 5). In details, every structure experiences quadratic curves. Thus, the relative position of two phases, in theoretical interface models and calculations, should be always taken into account, especially for the influence of d varying on interfaces where no obvious symmetries are adopted.

Notably, it is also observed that the interfacial energies of three interface pairs in Fig. 4 and Fig. 5 are quantitatively distinct from each other. This draws our attention to the descriptions of interfacial bonding. Table III lists the figures of interfacial bonding of various structures. Observations show obviously that the most stable interfacial structures exhibit a higher percentage of Mg-H bonds in specified range and a lower one of Mg-Mg bonds. Specifically, the interfacial atom relaxations allow Mg-H bonds to be in their most energetically favored length. Thus it indicates that the Mg-H bonding leads to higher interfacial stability. Comparatively, the Mg-Mg bonding is much less contributive to interfacial stability. Therefore we argue that the Mg-H bonding is the primary cause leading to higher interfacial stabilities. The most energetically favored interfacial bonding structure could be obtained by varying the relative position, i.e. the search of relative position of two phases.

Interestingly, we found that the interfacial stabilities are influenced by surface energies of two phases. The Mg(⁠|$10\bar 10$|$101¯0$⁠)/MgH₂(210), the interfacial pair of surfaces with high surface energy, are expected to have low interfacial energy due to great number of dangling bonds that each phase has. However, the results show that Mg(⁠|$10\bar 10$|$101¯0$⁠)/MgH₂(210) are owning a high interfacial energy around 4 meV/Ų (0.064 J/m²) while the Mg(0001)/MgH₂(101) pair consisted of two low energy surfaces has a much lower one of nearly 2.1 meV/Ų (0.034 J/m²). Meanwhile, the interfacial energy of Mg(⁠|$10\bar 10$|$101¯0$⁠)/MgH₂(101) pair, which is with one high energy and one low energy surfaces, is about 2.8 meV/Ų (0.045 J/m²). This difference of interfacial energies could be explained by the fact that the partial saturation of interfacial dangling bonds is not capable enough to offset the instabilities of high energy surfaces. Experimentally, it is observed⁴⁴ that the grain boundaries comprised of low energy surfaces have relatively low energies. Our finding is in line with the experimental ones, and indicates that consideration of surface energies should be taken into when choosing interface combinations.

The properties of Mg/MgH₂ interface, especially the hydrogen sorption property, are intriguing. As hydrogen sorption dovetails with the migration of Mg/MgH₂ interface, both of the hydrogen adsorption or desorption are considered. We here report comparative formation energies of hydrogen vacancies and interstitials, representing interfacial effect on hydrogen desorption and adsorption respectively.

Fig. 6(a) shows the interfacial effects on hydrogen vacancies. Clearly, it could be observed that hydrogen vacancy formations are highly facilitated in near-interface region. Overall, the formation of hydrogen vacancies at interface is much easier than that in MgH₂ bulk. It could be well understood as the formation of interface dovetails with some dangling (under-coordinated atoms) or floating (over-coordinated atoms) bonds. The existence of dangling or floating bonds means the interfacial H atoms are less bound by Mg atoms and weaker stabilities than those in MgH₂ bulk and thus the interface dramatically lowers the formation energies of hydrogen vacancies. Tensile strain monotonically increases the formation energies of |${\rm V}_{\rm H}^{{\rm 1st}}$|$VH1 st$ and |${\rm V}_{\rm H}^{{\rm 2nd}}$|$VH2 nd$ while the curves of far-interface |${\rm V}_{\rm H}^{{\rm 3rd}}$|$VH3 rd$⁠, |${\rm V}_{\rm H}^{{\rm 4th}}$|$VH4 th$⁠, |${\rm V}_{\rm H}^{{\rm 5th}}$|$VH5 th$⁠, |${\rm V}_{\rm H}^{{\rm 6th}}$|$VH6 th$ and |${\rm V}_{\rm H}^{{\rm 7th}}$|$VH7 th$ are slightly parabolic. When hydrogen vacancy is more distanced from interface, the Mg-H bonding approaches to that in MgH₂ bulk. Thus similar trends induced by interfacial effect on hydrogen vacancy formation could be clearly observed under various strain states. Biaxial strain effect affects the H vacancy formation energy as expected. As the interface lattice constants (tuned by strain) are close to those of MgH₂ bulk, the V_H formation energies go up due to increased stability of MgH₂ phase of interface thus parabolic curves appear. Interestingly, the MgH₂ bulk data is found to be monotonic, which is not parabolic. Actually, it is raised from difference between interface and bulk data. Specifically, the C axis of bulk is allowed to relax when strains are applied. It means all the atoms are able to relax in response to the strain. Two MgH₂ atom layers in the interface model are fixed in order to stimulate the stable MgH₂ substrate therefore the atom relaxations in far-interface regions are not enough so that similar parabolic curves could be observed. In authentic interface environments, the latter situation is more likely to happen since inadequate responses happen to a substrate that are ‘thick’ enough, which we usually deem to be strains. Based on the above analyses, we conclude that hydrogen vacancies are more energetically preferred in interface regions, especially the near-interface one, compared with the bulk region.

The interface also facilitates hydrogen interstitials. It is extremely easy for H_inte^1st to form, which is reflected by its negative value of formation energy. H interstitial formations are much more approachable to the bulk region situation as it is more distance from interface. The reason for this phenomenon is that Mg bulk context would facilitate hydrogen uptake thanks to the strong Mg-H bonding. The strong bonding strongly increases stability of system. When a tensile strain is applied, the space for H interstitials is enlarged and thus resulting in a monotonic trend of formation energies. It is observed that the lengths of Mg-H bond under strains of −6% and +6% are 1.92 Å and 1.96 Å for |${\rm H}_{{\rm inte}}^{{\rm 4th}}$|$H inte 4 th$⁠. The altering of Mg-H bonding rationalizes the monotonic trend of formation energies.

Three combinations of Mg/MgH₂ interface models are built to investigate the interface stabilities. The criterions including surface energies, interface matching, elastic effects and relative positions are discussed for Mg/MgH₂ interface modeling. It is found that the elastic effects and relative positions strongly affect the incoherent interface stabilities by affecting incoherent interfacial structure. Thus the searches of mutual constants and relative positions provide feasible solutions and make prediction of interfacial structure possible. Lastly, the defect formations are highly facilitated at Mg/MgH₂ interfaces.

This work is supported by MOST under project 2010CB631302, NSFC (Grant no. 11174082), the Fundamental Research Funds for the Central Universities (Grant No. 2013ZZ0082 and 2013ZM108), and KLGHEI (KLB11003). The computer times at National Supercomputing Center in Shenzhen (NSCCSZ) and ScGrid of the Supercomputing Center, Computer Network Information Center of CAS are gratefully acknowledged.