Global minimum beryllium hydride sheet with novel negative Poisson's ratio: first-principles calculations

Abstract

As one of the most prominent metal-hydrides, beryllium hydride has received much attention over the past several decades, since 1978, and is considered as an important hydrogen storage material. By reducing the dimensionality from 3 to 2, the beryllium hydride monolayer is isoelectronic with graphene; thus the existence of its two-dimensional (2D) form is theoretically feasible and experimentally expected. However, little is known about its 2D form. In this work, by a global minimum search with the particle swarm optimization method via density functional theory computations, we predicted two new stable structures for the beryllium hydride sheets, named α–BeH₂ and β–BeH₂ monolayers. Both structures have more favorable thermodynamic stability than the recently reported planar square form (Nanoscale, 2017, 9, 8740), due to the forming of multicenter delocalized Be–H bonds. Utilizing the recently developed SSAdNDP method, we revealed that three-center-two-electron (3c–2e) delocalized Be–H bonds are formed in the α–BeH₂ monolayer, while for the β–BeH₂ monolayer, novel four-center-two-electron (4c–2e) delocalized bonds are observed in the 2D system for the first time. These unique multicenter chemical bonds endow both α– and β–BeH₂ with high structural stabilities, which are further confirmed by the absence of imaginary modes in their phonon spectra, the favorable formation energies comparable to bulk and cluster beryllium hydride, and the high mechanical strength. These results indicate the potential for experimental synthesis. Furthermore, both α– and β–BeH₂ are wide-bandgap semiconductors, in which the α–BeH₂ has unusual mechanical properties with a negative Poisson's ratio of −0.19. If synthesized, it would attract interest both in experiment and theory, and be a new member of the 2D family isoelectronic with graphene.

---

1. Introduction

The prediction of new materials is always exciting,^1–3 especially when it has a reduced dimensionality, e.g. graphene^4–6 and graphene-like 2D materials, Cu₂X (Si, Ge),^7,8 which were named “graphene mimics of the future.” As the lightest element in the group of alkaline earth metals, beryllium is a magic element that not only can improve many physical properties when added as an alloying element to metal elements, but also can form both ionic and covalent bonds with many non-metal elements, such as hydrogen, carbon, and nitrogen. Among them, beryllium hydride (BeH₂) has been well investigated over the last 40 years as an effective hydrogen storage material.^9–11 Dating back to 1978, the crystalline BeH₂ was first prepared by G. J. Brendel et. al. under high-pressure compaction-fusion of amorphous BeH₂ catalyzed by lithium.¹² While suffering from a toxicity issue, bulk BeH₂ has the highest hydrogen weight percentage of 18.2 wt% and moderate hydrogen desorption temperature of about 250 °C.^13,14 Recently, its non-bulk forms, BeH₂ nanoparticles, and line-style [BeH₂]_n nanoclusters, have been successfully synthesized in experiment for hydrogen storage due to the large surface to volume ratios at the nanoscale level,^15–19 which raised further interest in its low dimensions such as two-dimensional (2D) form.

The concept of “2D beryllium hydride monolayer” was first proposed by M. Seel as early as in 1989,²⁰ thirteen years before the realization of the first atomic-thickness material graphene in 2004. Under the inspiration of the geometric structure of bulk LiBeH₃, this BeH₂ monolayer is supposed to be a very flat sheet, as shown in Fig. 1(c). After then, a debate on the stability and electronic property (metallic or insulating) for this planar BeH₂ monolayer had been developed.^21,22 Very recently, this planar square form of BeH₂ monolayer has again attracted interests from researchers due to its topological nontriviality.²³ In which, it has Fermi circles centered at Γ (0,0) and K (1/2,1/2) points, respectively, and the corresponding Fermi velocities are even higher than that in graphene. Then, a natural research question comes: is this planar form the global minimum structure for BeH₂ monolayer?

Fig. 1 Top (upper) and side (bottom) views of the geometric structure of (a) α–BeH₂, (b) β–BeH₂ and (c) γ–BeH₂ monolayers. Big blue balls represent Be atoms, and small white balls represent H atoms.

To answer this question, we first look into the bonding behaviors for Be and H elements themselves. On one hand, Be element is electron-deficient, which has an electron configuration of 1s²2s². Hydrogenation is an effective strategy to compensate the electron-deficiency of beryllium. The stoichiometric BeH₂ is isoelectronic with graphene,²⁴ thus, it is expected to be very stable in energy. The stability of graphene is represented by the delocalized six-center-two-electron (6c–2e) pairs.²⁵ On the other hand, Be element can not only form the traditionally chemical bonds but also can participle in the forming of the multicenter delocalized bonds, e.g. Be₅C₂ monolayer,²⁶ Be₂C monolayer.²⁷ In a recent literature, a tug-of-war between classical and multicenter bonding is demonstrated in Be containing H–(Be)_n–H species.²⁸ By sharing electrons between Be and H atoms, multicenter bonding would exist and then reduce the total energy significantly. In this way, low dimensionality and unique structural topology may lead to exceptional properties of materials, e.g. negative Poisson's ratio.

2. Computational details

Our density functional theory (DFT) calculations were performed by the Vienna ab initio simulation package (VASP).²⁹ The ion–electron interaction was treated by the projector-augmented plane wave (PAW) approach, and the electron exchange-correlation functional was described by the generalized gradient approximation (GGA) as expressed by Perdew, Burke, and Ernzerhof (PBE).³⁰ Since the standard PBE functional is unable to describe weak interactions correctly, we adopted a Grimme vdW corrected PBE approach by adding a semi-empirical dispersion potential to the conventional Kohn–Sham DFT energy.³¹ This approach introduces damped atom-pairwise dispersion correction of the form C₆R⁻⁶ in the DFT formalism. The accuracy of PBE has been well validated in recent works.³² A cutoff energy of 600 eV for the plane-wave basis set was adopted in all computations with the energy precision of 10⁻⁶ eV and force precision of 10⁻³ eV Å⁻¹. A 9 × 9 × 1 Γ-centered Monkhorst−Pack k-points grid is good enough for numerical convergence. A vacuum space 15 Å was adopted in the z-direction to prevent the interaction effect from adjacent cells.

The CALYPSO code³³ using a particle swarm optimization (PSO) method was employed to search for the lowest-energy BeH₂ monolayers. This PSO algorithm has been proved as an efficient structure prediction method in many works.³⁴ The required structure relaxations were performed by PBE functional, as implemented in VASP code. In our PSO calculations, the population size is set to 48, and the number of generation was maintained at 30. Besides the linear polymeric form of BeH₂, two low-energy stable BeH2 monolayers were obtained. To assess the kinetic stability of BeH₂ monolayers, the Phonopy code³⁵ was used to perform a phonon dispersion analysis with density functional perturbation theory (DFPT),³⁶ as implemented in VASP. In the VASP-DFPT lattice dynamics calculation, a 4 × 4 × 1 q-point grid is employed to compute the interatomic force. To evaluate the thermal stability, the ab initio molecular dynamics (AIMD) was carried out with a Nosé−Hoover algorithm.³⁷ To be accurate and reliable, a 5 × 5 × 1 supercell of Be₂₅H₅₀ is used with a total simulation time of 20 ps.

3. Results and Discussion

Following the “bottom-up” analysis strategy, we first study the simplest beryllium hydride, Be₂H₄ molecule. With two Be, two terminal and two bridging H atoms, the ground state of the Be₂H₄ molecule has a D₂h structural symmetry. A total of four electrons are shared between two Be atoms and two bridging H atoms, resulting in two three-center-two-electron (3c–2e) bonds, as well known as “banana” bonds. Inspired by the work of the borane monolayer constructed from the bonding nature of diborane molecules,³⁸ we designed a BeH₂ chain from the bonding nature of the Be₂H₄ molecule. It is constructed by both the “banana” bonds and the covalent Be–Be bonds. This BeH₂ chain is easy to be expanded to a ring type or a link type configurations, while not in a periodic 2D form. To solve this issue, a comprehensive first-principles-based particle swarm optimization (PSO) search is used to predict the global minimum structures for the 2D BeH₂. The optimized geometric structures for two most stable configurations are shown in Fig. 1, namely α– and β–BeH₂. For comparison, the planar form of the BeH₂ monolayer is also given here, namely γ–BeH₂.²⁰

The symmetric groups for α–, β– and γ–BeH₂ are P–4M2, P–3M1 and P4/MMM, with the lattice constants of 2.38, 2.30 and 2.95 Å, respectively. According to our careful computations, α–BeH₂ has the lowest total energy, which is 0.009 eV per unit cell lower than that of β–BeH₂. Since the room temperature (300 K) corresponds to thermal energy of K_BT ≈ 0.026 eV, both α– and β–BeH₂ may exist as stable structures. In comparison, γ–BeH₂ is 1.3 eV higher in energy than those of α– and β–BeH₂. As zero-point energy (ZPE) could be relevant, we then included the ZPE to the structural energies, and found that the β–BeH₂ is slightly more stable than α–BeH₂ by 0.16 eV in energy. Due to the complexity of the real experimental environment, we believed that α– and β–BeH₂ are equally important and thus discussed them together in this work.

In the α–BeH₂, as similar as in γ–BeH₂, the Be atoms are tetra-coordinated with four H atoms around it. However, not like all the atoms in the same plane in γ–BeH₂, the configurations of α–BeH₂ could be seen as one Be atomic layer sandwiched between two H atomic layers. In this way, α–BeH₂ can be considered as a distorted γ–BeH₂. The H–Be–H bond angle is increased from 90° in γ–BeH₂ to 109.4° in α–BeH₂, meanwhile, the bond length of Be–H is slightly decreased from 1.48 to 1.46 Å, which is still longer than that of 1.43 Å in bulk BeH₂.¹⁶ As a result, the structural distortion of the α–BeH₂ results in a much higher stability. In the β–BeH₂, each Be atom is hexacoordinate with six neighboring H atoms, and each H atom is bonded to three Be atoms on both sides of the plane in an alternating manner. Interestingly, the bond length of Be–H in the β–BeH₂ is computed to be 1.59 Å, which is much longer than that in bulk BeH₂ (1.43 Å) at the same theoretical levels. To the best of our knowledge, this is the longest Be–H bond length which has been reported to date.

So, why the structures of α–BeH₂ and β–BeH₂ are more stable than that of γ–BeH₂? Even the former structures have much longer Be–H bond lengths? To understand this, we evaluated the cohesive energies to access the experimental feasibility of these two predicted BeH₂ monolayers, defined as E_c = (n_BeE_Be + n_HE_H − E_BeH2)/(n_Be + n_H). In which, E_BeH2 is the energy of α– or β–BeH₂. E_Be and E_H are the energies of single Be and H atoms. In this definition, a larger E_c represents a higher stability. According to our computations, the E_c values of α–BeH₂ and β–BeH₂ are similar as much as 2.76 eV per atom, which are 0.14 eV per atom smaller than that of bulk BeH₂ (2.90 eV per atom), higher than that of the experimentally synthesized nanoclusters Be_nH_2n (n = 3–9) (2.54–2.71 eV per atom),¹⁶ and much higher than that of 2.34 eV per atom in the γ–BeH₂. For comparison with other popular 2D materials, using the same theoretical method, the cohesive energies of graphene, Cu₂Si monolayer, and germanene are 7.85, 3.46 and 3.26 eV per atom, respectively. Thus, both α– and β–BeH₂ have high structural stabilities.

To access the kinetic stabilities of α– and β–BeH₂, we then calculated their vibrational normal modes. As shown in Fig. 2, there are no appreciable imaginary modes in the phonon spectrums in the α– and β–BeH₂, indicating high kinetic stabilities of these two nanomaterials. Remarkably, the highest frequencies of α–BeH₂ and β–BeH₂ are ∼1785 cm⁻¹ (= 53.54 THz) and ∼1470 cm⁻¹ (= 44.09 THz), respectively, which are comparable to the highest frequency of 1600 cm⁻¹ in graphene,³⁹ showing the robust bonds in both α– and β–BeH₂.

Fig. 2 Phonon spectrums of (a) α–BeH₂ and (b) β–BeH₂ monolayers with no appreciable imaginary modes.

Then, the thermodynamic stabilities of α– and β–BeH₂ have been confirmed through the AIMD simulations. Three different temperatures of 298, 500 and 550 K with a simulation time of 10 ps were performed. Snapshots of the geometry structures of both α– and β–BeH₂ at the end of each AIMD simulation are presented in Fig. 3. Our computations demonstrated that the structures of α– and β–BeH₂ can be well kept throughout 20 ps AIMD simulation up to 500 K, but collapses at 550 K, indicating a melting point in the range of 500–550 K. This melting point could be approximately considered as the hydrogen desorption temperature, which is comparable to that of bulk beryllium hydrides at 523 K.⁴⁰ After full atomic relaxation, the final structures of α– and β–BeH₂ at 298 and 500 K can easily recover their initial configurations. These results suggest that both α– and β–BeH₂ have good thermodynamic stabilities. Furthermore, since 500–550 K is moderate hydrogen desorption temperature, this result renders BeH₂ monolayers as promising hydrogen storage materials.

Fig. 3 Top (upper) and side (bottom) views of geometric configuration of α–BeH₂ and β–BeH₂ monolayers after 20 ps AIMD simulations at (a and d) 298 K, (b and e) 500 K and (c and f) 550 K, respectively. Big blue balls represent Be atoms, and small white balls represent H atoms.

To understand the high structural stabilities for α– and β–BeH₂, a deep insight into the chemical bonding patterns were analyzed by the recently developed SSAdNDP method.^41–43 Since the electron configurations of Be and H atoms are 1s²2s² and 1s¹, respectively, there are four valence electrons in one unit cell of BeH₂. According to our computational results, multicenter delocalized bonds are obtained in both of these two configurations. In the α–BeH₂ (see Fig. 4(a)), the bonding within Be–H–Be ring is proved to be a 3c–2e bond, with an occupation number of 1.94 |e|, well known as “banana” bonds. And there are two “banana” bonds on both sides of the Be atomic layer, accounting for total 4 valence electrons. This novel “banana” bond formed in the α–BeH₂ can also be found in the [BeH₂]_n nanoclusters.¹⁶ While in the β–BeH₂ (see Fig. 4(b)), the bonding nature is quite unique: the bonds distributed between the Be atomic layer and two H atomic layers are turned out to be two 4c–2e bonds. The occupation number of 4c–2e bonds is 1.96 |e|. To the best of our knowledge, this is the first report of 4c–2e bonds in the 2D beryllium hydride systems, quick after its first finding in beryllium hydride cluster of Be₄H₂ in 2008.²⁸ It will explain the high stability for the longest Be–H bond (1.59 Å). In summary, the above results reveal that the forming of the “banana” bonds in α–BeH₂ and 4c–2e bonds in the β–BeH₂ are responsible for their high structural stabilities. As a building block, these multicenter bonds endow the α– and β–BeH₂ with a 2D network structure. Furthermore, the Hirshfeld charge population analysis showed that the net charges on Be in the α– and β–BeH₂ are only 0.15 and 0.10 e, respectively, indicating that the Be–H chemical bonds are mainly of covalent nature.

Fig. 4 Schematic of SSAdNDP chemical bonding pattern for unit cells of (a) α–BeH₂ and (b) β–BeH₂ monolayers. Big blue balls represent Be atoms, and small white balls represent H atoms. “ON” represents the occupation number. The iso-surface value is 0.10 e Å⁻³.

Since the α– and β–BeH₂ are fully hydrogenated, both configurations are expected to have a wide bandgap. For example, hydrogenation of the graphene, silicene, germanene would significantly induce these zero bandgap semiconductors to wide bandgap semiconducting for the graphane (5.40 eV), silicane (4.07 eV) and germanane (2.37 eV).⁴⁴ We then calculated the band structures to examine the electronic properties of the α– and β–BeH₂. Unlike the metallic nature in the γ–BeH₂,²³ as shown in Fig. 5, both the α– and β–BeH₂ have very large band gaps of 4.68 and 4.88 eV, respectively. In the band structure of the α–BeH₂, the valence band maximum (VBM) is located at the M point and the conduction band minimum (CBM) is located along the Γ–X line, showing an indirect bandgap semiconducting feature. While for the β–BeH₂, it is a direct bandgap semiconductor due to the VBM and CBM locating at the same symmetric point along the Γ–M line. It needs to be pointed out that the large band gaps of the α– and β–BeH₂ restrict their applications in microelectronic and optoelectronic devices, thus, further band modulations are needed, e.g. doping, edge decorations.

Fig. 5 Electronic band structures of (a) α–BeH₂ and (b) β–BeH₂ monolayers at the PBE level.

The mechanical properties are an important parameter for the potential application of 2D materials. Then, we examined the mechanical properties of both α– and β–BeH₂ by computing the elastic constants. To be mechanical stable, in theory, a 2D structure should meet the following criteria: C₁₁C₂₂–C₁₂² > 0, C₆₆ > 0, in which C_ij are the elastic constants. According to our calculations, the elastic constants of α–BeH₂ are C₁₁ = C₂₂ = 62.24 GPa nm, C₁₂ = C₂₁ = −11.84 GPa nm, and C₆₆ = 8.21 GPa nm, which satisfy well with the mechanical stability criteria. The in-plane Young's modules for the α–BeH₂ can be computed from the elastic constants by Y = (C₁₁C₂₂–C₁₂C₂₁)/C₂₂, which is as much as 59.98 N m⁻¹. For the β–BeH₂, the elastic constants are calculated to be C₁₁ = C₂₂ = 84.97 N m⁻¹, C₁₂ = C₂₁ = 13.98 N m⁻¹, C₆₆ = 35.50 N m⁻¹, and the corresponding in-plane Young's modulus is 82.67 N m⁻¹. The in-plane stiffness of α– and β–BeH₂ are both lower than that of graphene (355.35 N m⁻¹). However, they are comparable to that of Cu₂Si monolayer (Y = 93 N m⁻¹), 2D silicene monolayer (Y = 71.20 N m⁻¹),⁴⁵ and higher than that of germanene (42 N m⁻¹) computed at the same theoretical level. This result indicates the strong mechanical stabilities for both α– and β–BeH₂ monolayers.

Remarkably, we noticed that α–BeH₂ has a negative C₁₂, which results in negative Poisson's ratios of −0.19 (C₁₂/C₂₂ = C₁₂/C₁₁). Since Poisson's ratio is defined as the ratio of a lateral contraction to the longitudinal extension during the stretching of a material, the negative Poisson's ratio of α–BeH₂ indicates a lateral extension in response to stretching, which is generally believed to be rare in 2D materials. For validation, a uniaxial stretching of 5% is applied for the α–BeH₂ in the x-direction. According to expectation, the equilibrium lattice constant of the α–BeH₂ in the y-direction is elongated by ∼0.08%. Similarly, a uniaxial compressing of 5% in the x-direction shorters the equilibrium lattice constant by ∼0.1% in the other direction. This unusual mechanical response confirms α–BeH₂ has a negative Poisson's ratios, which render its potential applications of mechanics.

4. Conclusions

In summary, we extended the beryllium hydride system from bulk and linear form to 2D form by designing two new BeH₂ monolayers. Our computations demonstrated that both BeH₂ monolayers are structurally stable, as indicated by the absence of imaginary modes in phonon spectra, moderate melting points, favorable formation energies and high mechanical strengths. According to the SSAdNDP analysis, the “banana” bonds in α–BeH₂ and 4c–2e bonds in the β–BeH₂ are formed and are responsible for their high structural stabilities. Remarkably, the α–BeH₂ monolayer has a rather intriguing mechanical property featured with a negative Poisson's ratios. Since the BeH₂ monolayers are isoelectronic to graphene, the existence of them could be experimentally confirmed in the future. They would have a potential in the hydrogen storage due to the high hydrogen weight percentage of 18.2 wt% and moderate hydrogen desorption temperature 500–550 K.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Supported by National Natural Science Foundation of China (No. 21522305, 61605087 and 61704083), Natural Science Foundation of Jiangsu Province (No. BK20160881, 16KJB140010), and China Postdoctoral Science Foundation (No. 20161001). The authors gratefully acknowledge the computing time granted by the VSR commission on the supercomputer JURECA at Forschungszentrum Jülich and Shanghai Supercomputer Center.

Notes and references

1. M. L. Cohen, Science, 1993, 261, 307–308.
2. S. Serra, C. Cavazzoni, G. L. Chiarotti, S. Scandolo and E. Tosatti, Science, 1999, 284, 788–790.
3. J. O. Sofo, A. S. Chaudhari and G. D. Barber, Phys. Rev. B, 2007, 75, 153401.
4. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Science, 2004, 306, 666–669.
5. J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. J. Booth and S. Roth, Nature, 2007, 446, 60–63.
6. A. K. Geim and K. S. Novoselov, Nat. Mat., 2007, 6, 183–191.
7. L. M. Yang, I. A. Popov, A. I. Boldyrev, T. Heine, T. Frauenheim and E. Ganz, Phys. Chem. Chem. Phys., 2015, 17, 17545–17551.
8. L. Yang, V. Bačić, I. A. Popov, A. I. Boldyrev, T. Heine, T. Frauenheim and E. Ganz, J. Am. Chem. Soc., 2015, 137, 2757–2762.
9. H. Wu, ChemPhysChem, 2008, 9, 2157–2162.
10. E. N. Koukaras, A. D. Zdetsis and M. M. Sigalas, J. Am. Chem. Soc., 2012, 134, 15914–15922.
11. N. V. Belkova, L. M. Epstein, O. A. Filippov and E. S. Shubina, Chem. Rev., 2016, 116, 8545–8587.
12. G. J. Brendel, E. M. Marlett and L. M. Niebylski, Inorg. Chem., 1978, 17, 3589–3592.
13. M. G. Ganchenkova, V. A. Borodin and R. M. Nieminen, Phys. Rev. B, 2009, 79, 134101.
14. A. Allouche, M. Oberkofler, M. Reinelt and C. Linsmeier, J. Phys. Chem. C, 2010, 114, 3588–3598.
15. A. D. Zdetsis, M. M. Sigalas and E. N. Koukaras, Phys. Chem. Chem. Phys., 2014, 16, 14172–14182.
16. E. N. Koukaras, A. P. Sgouros and M. M. Sigalas, J. Am. Chem. Soc., 2016, 138, 3218–3227.
17. A. Abdurahman, J. Phys. Chem. A, 2003, 107, 11547–11552.
18. A. Abdurahman, A. Shukla and M. Dolg, J. Chem. Phys., 2000, 112, 4801–4805.
19. C. B. Lingam, K. R. Babu, S. P. Tewari and G. Vaitheeswaran, Comput. Theor. Chem., 2011, 963, 371–377.
20. M. Seel, A. B. Kunz and S. Hill, Phys. Rev. B, 1989, 39, 7949–7954.
21. J. Z. Wu, S. B. Trickey and J. C. Boettger, Phys. Rev. B, 1990, 42, 1663–1667.
22. M. Seel, Phys. Rev. B, 1991, 43, 9532–9537.
23. B. Yang, X. M. Zhang and M. W. Zhao, Nanoscale, 2017, 9, 8740–8746.
24. Y. L. Jiao, F. X. Ma, J. Bell, A. Bilic and A. J. Du, Angew. Chem., Int. Ed., 2016, 55, 10292–10295.
25. I. A. Popov, K. V. Bozhenko and A. I. Boldyrev, Nano Res., 2012, 5, 117–123.
26. Y. Wang, F. Li, Y. Li and Z. Chen, Nat. Commun., 2016, 7, 11488.
27. Y. F. Li, Y. L. Liao and Z. F. Chen, Angew. Chem., Int. Ed., 2014, 53, 7248–7252.
28. K. A. Lundell and A. I. Boldyrev, Chem. Phys. Lett., 2018, 699, 85–87.
29. G. Kresse and J. Furthmuller, Phys. Rev. B, 1996, 54, 11169–11186.
30. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
31. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
32. Y. F. Li, F. Y. Li and Z. F. Chen, J. Am. Chem. Soc., 2012, 134, 11269–11275.
33. Y. Wang, J. Lv, L. Zhu and Y. Ma, Phys. Rev. B, 2010, 82, 094116.
34. T. Gu, W. Luo and H. J. Xiang, Wires Comput. Mol. Sci., 2017, 7, e1295.
35. A. Togo, F. Oba and I. Tanaka, Phys. Rev. B, 2008, 78, 134106.
36. X. Gonze and C. Lee, Phys. Rev. B, 1997, 55, 10355–10368.
37. G. J. Martyna, M. L. Klein and M. Tuckerman, J. Chem. Phys., 1992, 97, 2635–2643.
38. T. A. Abtew, B.-c. Shih, P. Dev, V. H. Crespi and P. Zhang, Phys. Rev. B, 2011, 83, 094108.
39. C. Si, W. Duan, Z. Liu and F. Liu, Phys. Rev. Lett., 2012, 109, 226802.
40. M. T. Kelly, in Fuel Cells and Hydrogen Storage, ed. A. Bocarsly and P. D. M. Mingos, Springer Berlin Heidelberg, Berlin, Heidelberg, 2011, pp. 169–201.
41. T. R. Galeev, B. D. Dunnington, J. R. Schmidt and A. I. Boldyrev, Phys. Chem. Chem. Phys., 2013, 15, 5022–5029.
42. D. Y. Zubarev and A. I. Boldyrev, Phys. Chem. Chem. Phys., 2008, 10, 5207–5217.
43. B. D. Dunnington and J. R. Schmidt, J. Chem. Theory Comput., 2012, 8, 1902–1911.
44. W. Wei, Y. Dai, B. B. Huang and T. Jacob, Phys. Chem. Chem. Phys., 2013, 15, 8789–8794.
45. Q. Peng, X. D. Wen and S. De, RSC Adv., 2013, 3, 13772–13781.

---