Density functional theory simulations of complex hydride and carbon-based hydrogen storage materials

Abstract

This critical review covers the mechanisms underlying density functional theory (DFT) simulations and their relevance in evaluating, developing and discovering new materials. It is intended to be of interest for both experimentalists and theorists in the expanding field of hydrogen storage. We focus on the most studied classes of materials, metal-hydride, -amide, and -borohydride mixtures, and bare and transition metal-doped carbon systems and the utility of DFT simulations for the pre-screening of thermally destabilised reaction paths (170 references).

---

Stephen Shevlin

Dr Stephen Shevlin received an MSci from the University of Birmingham and a PhD from University College London. He has worked extensively in the field of electronic structure methods, performing research on novel dielectric materials, metal nanoclusters, and new hydrogen storage materials. He has been based at IBM Zurich Switzerland, Oak Ridge National Laboratory, and Queen Mary, University of London. Currently he is working with Xiao Guo at the Department of Chemistry, University College London.

Z. Xiao Guo

Z. Xiao Guo is a Professor of Materials and Chemistry at University College London (UCL). He received a PhD from the University of Manchester, and then worked as a research scientist at the Universities of Strathclyde and Oxford. He was a Lecturer and then a Professor at Queen Mary, University of London before joining UCL. He has contributed over 130 journal publications and a similar number of conference presentations. He has received the Beilby Medal from the Society of Chemical Industry, the Royal Society of Chemistry, and the Institute of Minerals, Metals and Materials. He is on the editorial boards of several international journals.

---

1. Introduction

One of the major issues in the 21st century is the generation and utilization of non-polluting energy. At present, the consumption of energy is intimately linked to CO₂ emission, a major and increasingly prevalent artificial contributor to climate change. Several solutions can ameliorate this: including increases in efficiency, generation by renewable sources, and increasing usage of nuclear power. However fossil fuels with a high energy density will still be dominant in the near term. If CO₂ emission is to be significantly reduced there should be a switchover to a clean alternative, of which hydrogen (H₂) is the preferred candidate.^1,2Hydrogen however is difficult to store safely at low cost, light weight, and low volume. Thus storage in solid materials that have fast H₂ release and adsorption kinetics at 85 °C (a typical operating temperature) is desired.^1–4

Simulations can aid in the understanding and discovery of hydrogen storage materials. In this review, we will briefly outline the fundamentals underpinning density functional theory (DFT), a commonly used technique for modelling molecular and solid systems as it is faster than wavefunction-based ab initio methods, albeit at the cost of reduced accuracy. We will then discuss how this technique may be applied to calculate the following quantities of interest: hydrogen binding energies, enthalpies and free energies of formation, diffusion and reaction pathways, activation barriers, and transition states. Finally, we will review the literature on simulations of two important classes of hydrogen storage materials, metal hydride mixtures and doped carbon systems.

2. Theoretical background

2.1 Density functional theory (DFT)

The properties of a material are determined by the behaviour of the electrons that comprise the bonds of the system in question. In DFT the ground state of an interacting electron gas is mapped onto the ground state of a non-interacting electron gas which experiences an effective potential, in contrast to other techniques such as Hartree–Fock which solve for the full set of wavefunctions in the Schrödinger equation. This mapping in principle gives the exact ground-state properties (such as cohesive energy, lattice parameters and phonon spectra). The Hamiltonian that describes the electron gas density is
where T̂ and V̂_ee are the kinetic and electron–electron interaction operators, respectively and V_ext is the external potential. The two basic theorems of DFT are

1. There exists a functional F[n] of the electronic density n that represents the ground state energy E_GS.

2. The only density which gives an energy equal to the energy of the ground state is the electronic density of the ground state n_GS, with all other densities giving higher energies

∫drV_ext(r)n_GS(r) + F[n_GS] ≡E_GS

The functional E[n] thus satisfies a variational principle.

We can use a separation of the Hamiltonian first devised by Kohn and Sham,⁵ wherein we can write the energy as a sum of functionals that depend on electron density

E[n(r)] = T₀[n(r)] + ∫dr n(r)[V_ext(r) + ½Φ(r)] + E_XC[n(r)]

where T₀ is the kinetic energy of a system with density n in the absence of electron–electron interactions, Φ is the classical Coulomb potential and E_XC is the exchange and correlation energy. The exchange energy is the contribution to the energy that results from the Pauli principle, that no two electrons of the same spin energy can reside in the same state. The correlation energy is the contribution to the energy from the Coulomb repulsion that an electron experiences from all other electrons.

The exchange and correlation energies thus depend on the electronic density. As it is not known how to solve the n-body problem exactly we have no idea a priori how E_XC[n] behaves. The oldest and simplest way of solving the variational equations is the local-density approximation (LDA).⁶ This assumes that the exchange–correlation potential for an electron is equal to the exchange–correlation potential of an electron in a uniform gas of interacting electrons, and can be inserted into the above equation as

E^LDA_XC = ∫drn(r)ε_XC[n(r)]

where ε_XC[n(r)] is the exchange–correlation energy density of a homogeneous electron gas. In actuality, although the LDA is surprisingly accurate for many systems, it has been surpassed by newer and more sophisticated functionals. A useful scheme for classifying all of these functionals is to divide them into classes based on the types of variables that are involved, the so-called “Jacob’s ladder” of functionals.⁷ At the lowest level, the exchange–correlation functional depends solely on the electron density n, e.g. the LDA. The next level up are the generalised gradient approximations (GGA).^8–12 We note that in this review we discuss in more depth functionals that are nonempirical, that is functionals that are constructed by constraint fitting to known properties. In contrast, semiempirical functionals such as BLYP, B3LYP etc. are fitted to data sets provided by high quality wavefunction based methods.

Historically the GGAs originated in the presumption that as real systems are not homogeneous then improvements to the LDA can be made by treating the electron gas as slowly varying and expanding to second order, the gradient expansion approximation (GEA)¹³

U^GEA[n(r)] = ∫(a(n(r)) + b(n(r))|∇n(r)|² + ⋯)

Unfortunately this expansion does not lead to improvements, as in reality the electron gas is rapidly varying so the expansion breaks down at large distances due to large oscillations in the electron density.¹⁴ At shorter distances however the GEA is an improvement. Therefore in order to utilise the advantages of the GEA there needs to be a truncation at the point where unphysical oscillations in the exchange–correlation density manifest. This truncation is known as the GGA. Placing exchange and correlation into separate terms, for any GGA the exchange energy may be written as

E^GGA_X = ∫d³r ne^unif_x(n)F_X(s(r))

where F_X(s(r)) is the factor of exchange enhancement caused by the GGA and s(r) = |∇n(r)|/2k_Fn is a dimensionless density gradient parameter dependent on the Fermi wavenumber of the material k_F.¹² To recover the correct uniform gas limit, set F_X(0) = 1. For the correlation energy there is simply a renormalization of all radii such that r is now set to r_s. The full GGA correction to the exchange–correlation energy E_XC is

E^GGA_XC[n] = ∫d³r ne^unif_x(n)F_XC(r_s,s)

and is dependent on both n and ∇n. Different flavours of GGA are simply choices on how this truncation is performed. The Perdew–Wang (PW91) GGA is an analytic fit to this nonempirical GGA^10,13 however the parameterisation is complex. For example, the exchange enhancement factor in this functional is
The Perdew, Burke and Enzerhof (PBE) GGA¹² is a modification of this nonempirical functional that only satisfies conditions that are energetically relevant. For example, the exchange enhancement factor in the PBE is
leading to a much smaller and more transparent parameterisation.

The next rung in the ladder of increasing accuracy of density functionals are the meta-GGAs. As the LDA is a functional that depends only on n, and the GGAs are functionals that depend on n and ∇n, meta-GGAs are functionals that are often of the form¹⁵

E^MGGA_XC[n] = ∫d³r ε^MGGA_XC(n,∇n,∇²n,τ)

where τ are the Kohn–Sham orbital kinetic energy densities
which appear in the Taylor expansion of the exchange correlation hole and thus are implicit functionals of the density. Several meta-GGAs have been developed, but the only nonempirical functional is the TPSS functional of Tao, Perdew, Staroverov, and Scuseria.¹⁵

Higher rungs up the ladder of functionals would introduce more nonlocality. For example, it has been posited that the next rung up would be the so-called “hyper-GGA” functionals which would include exact exchange energy density.⁷ Currently there do not exist any nonempirical formulations of a hyper-GGA. “Hybrid” functionals, which admix a portion of exact Hartree–Fock exchange into the exchange–correlation functional,^16,17 can be regarded as analogous functionals to the hyper-GGA. These tend to be the functionals of choice for the chemistry community and often are the most accurate functionals for calculations of barrier heights. However these functionals are optimised semi-empirically on a specific set of small molecules, therefore transferability to other systems may be constrained. Indeed for large molecules B3LYP performs less well than alternative functionals.¹⁸ A summary of all of these different types of functionals is illustrated in Fig. 1.

Fig. 1 Schematic of relationship between commonly used exchange correlation functionals, increasing sophistication of treatment of exchange and correlation, and accuracy. Figure adapted from Fig. 1 of ref. 7.

2.2 Relevant properties for hydrogen storage

Density functional theory can accurately determine energetics and structures, and on modern day parallel computers is able to simulate fairly large (few hundred atom) systems in a reasonable amount of time (of order days). It is thus the method of choice for most hydrogen storage simulations. Important quantities of interest for hydrogen storage are the strength of hydrogen binding (and the associated enthalpy of formation and plateau pressure for hydrogen release) and the kinetics of H₂ release.

The strength of H₂ binding to the material is of importance for hydrogen storage applications. This is defined as

E_Bind(H₂) = E_Tot(X− H₂) + E_Tot(H₂) −E_Tot(X)

where E_Tot(X) is the total energy of the original system, E_Tot(X− H₂) is the total energy of the original system after two hydrogen atoms are removed, and E_Tot(H₂) is the total energy of a H₂ molecule in vacuum, chosen as the reference state. All geometries are relaxed. Often the reference is chosen as the energy of a single atom in vacuum, however as experimentally hydrogen is introduced into these systems in molecular form a more reasonable choice is the molecule reference state. The more positive E_Bind(H₂) the stronger hydrogen is bound. In order to meet the US Department of Energy (DOE) targets for hydrogen storage applications the binding energy (sometimes referred to as the removal energy or equivalently the absorption energy) should be 0.7 eV/H₂≥E_Bind(H₂) ≥ 0.2 eV/H₂.^19,20 This is true regardless of whether the hydrogen is absorbed in a material (enters a material) or adsorbs to a material (binds to the surface). We emphasize that this formulation of the binding energy involves only the electronic contribution to the energy, and thus is explicitly a zero temperature binding energy. However as for large nanostructures it is prohibitively expensive to calculate phonon modes, the electronic binding energy is a suitable measure for hydrogen affinity. Lee et al. have set out a criterion for thermodynamically usable hydrogen capacity f
where ε_l is the zero-point corrected binding energy, l is the number of adsorbed H₂ molecules, k is Boltzmann’s constant, g_l is the multiplicity of a configuration for a given l, and μ is the chemical potential of H₂ in the gas phase at temperature T.²¹ This allows a direct comparison of usable capacity under realistic operating conditions.

For smaller nanostructures, and for periodic solids with primitive unit cells of less than fifty atoms, it is computationally tractable to calculate phonon modes and thus the full enthalpy of formation for successive or complete dehydrogenation steps. The temperature dependent free energy of a solid can be written as²²

G(T) = U(T) −TS(T) + PV≈U₀ + U_vib(T) −TS_vib(T)

where U₀ is the electronic energy, U_vib(T) and S_vib(T) are the vibrational contributions to internal energy and entropy, respectively, while the pressure–volume term is neglected. For the H₂ gas, the translational and rotational contributions to the internal energy (5/2RT) and the PV terms (RT) need to be included. The entropy of the H₂ gas can be taken from standard thermochemical sources, and is usually taken to be 130 J K⁻¹ mol⁻¹ (0.0013 eV K⁻¹) at 300 K and 1 atmosphere.^4,22,23 Experimentally the enthalpy change, ΔH = TΔS, is measured. If the change in entropy is dominated by the H₂ gas, then hydrogen release in the temperature range of 300 to 600 K requires an enthalpy change of between 30 and 78 kJ mol⁻¹ (0.31 to 0.81 eV). If the activity of all the solid phases is taken to be unity, then the equilibrium pressure P of dehydrogenation can be taken as
Plotting these pressures defines the van’t Hoff plot for decomposition.^24–26 Ideally hydrogen release should occur at relatively high pressures, greater than 1 bar.

In addition to favourable thermodynamics it is desirable for hydrogen storage materials to possess quick charge and discharge kinetics, of the order of seconds. In order to determine the kinetics it is necessary to calculate the diffusion rate r through the harmonic transition state expression²⁷

where q_IS and q_TS are the partition functions of the initial (IS) and transition (TS) states, h is Planck’s constant, and ΔE_a is the activation energy (the energy difference between initial and transition states). The most sensitive component of this diffusion rate is the exponential term dependent on barrier height, and catalysts for improving the kinetics of a reaction typically act to reduce the barrier height, e.g.Ti substitutionals in NaAlH₄.²⁸ The harmonic approximation should only work well at low temperatures, but is often applicable at high temperatures as well.^29,30 In principle, as H₂ has a large zero-point energy, quantum effects can still be of importance and might need to be included in the partition function.³¹ In order to calculate the minimum energy pathway, and thus the barrier for activated processes, a commonly used method is the nudged elastic band (NEB) method,^32–34 where multiple intermediate configurations between the initial and final states are connected via fictitious springs and relaxed until all forces perpendicular to the path vanish. Such calculations can be expensive and time-consuming, due to the need for such NEB simulations to be performed on large-scale parallel computers and with convergence tending to be a quite prolonged process. Calculations of the kinetics via the transition state theory are straightforward for systems where the hydrogen is absorbed, but materials where the hydrogen is adsorbed (e.g.metal hydrides) where the stoichiometry changes between initial and final products are too complex for straight forward applications of the theory to apply.

We have thus enumerated four types of quantity of interest for hydrogen storage. These can be divided into:

• Structures and electronic binding energies.

• Phonons, free energies, enthalpies of dehydrogenation and entropies.

• Transition states and activation barriers.

• Diffusion and rate constants.

Typically simulations focus on the first set of items. However, all four sets of information are of potential interest for hydrogen storage.

3. Solid hydride systems

One class of hydrogen storage systems that are of particular interest are those that bind hydrogen in chemical form, e.g.metal hydrides, amides, alanates, and borohydrides. In this section we will report on the current status of theoretical simulations on representative examples of these materials. We will also discuss the concept of thermodynamic destabilisation, an important concept for hydrides research, and computational efforts to explore new reaction paths utilizing this concept.

3.1 Magnesium hydride (MgH₂)

Magnesium hydride (MgH₂) has long been considered as a promising material for hydrogen storage due to substantial thermal stability, low cost, reversibility of the reaction MgH₂↔ Mg + H₂, and significant gravimetric hydrogen storage capacity (7.6 wt%). However, due to the high desorption temperature (between 300 to 400 °C), a low plateau pressure of 10 Pa (or 10⁻⁴ bar) at ambient, and relatively slow absorption/desorption kinetics,³⁵ substantial effort, including simulation work, has been devoted to improving this material for hydrogen storage applications.

At room temperature, rutile MgH₂ has a tetragonal unit cell with lattice constants a = 4.518 and c = 3.021(c) Å, see Fig. 2.³⁶ This results in the (110) surface having the lowest surface energy and the lowest hydrogen binding energy, which DFT-PBE and DFT-PW91 finds to be of magnitude 1.2 to 1.5 eV/H₂.^37–39 Both experimentally⁴⁰ and theoretically⁴¹ the stability of MgH₂ can be affected by particle size, with decreasing particle size decreasing the stability of the hydride. However, in order to get a significant decrease in their thermodynamic stability and hence desorption temperature it is necessary to reduce the grain size to less than 1.3 nm,⁴¹ unlikely to be achievable using mechanical milling. In principle, alternative metastable MgH₂ phases with lower energies of formation, for example an anatase structure with a formation energy of 58.03 kJ mol⁻¹ (0.60 eV) compared with the rutile formation energy of 76.00 kJ mol⁻¹ (0.79 eV) as found in DFT-PBE simulations,⁴² may be utilized, but such metastable phases have not yet been identified in practice. Therefore most approaches to reducing the desorption temperature of MgH₂ have focused on dopants that alter the binding between the Mg and H species. Pozzo et al. have recently published a paper comparing the accuracy of DFT simulations of the formation enthalpy of MgH₂ with quantum Monte Carlo simulations, a very accurate but very computationally expensive technique.⁴³ They find that the GGA functionals give an enthalpy of formation that is in worse agreement with experiment than the LDA. Both types of exchange correlation functional however give enthalpies of formation within ∼0.3 eV f.u.⁻¹ (formula unit) of experiment.

Fig. 2 Crystal structure of MgH₂, where the Mg and H atoms are represented by light blue and dark red, respectively. Figure taken with permission from ref. 43. Copyright (2008) by the American Physical Society (URL: http://link.aps.org/abstract/PRB/v77/e104103).

The differing electronegativities of Mg and H atoms have significant implications for the binding of MgH₂. In principle, Mg should donate charge density to the H atoms forming an ionic bond, see Fig. 3. Simulations indicate that this is partially the case,^44,45 with the Mg atoms having a positive charge of +1.95e. However the H atoms do not accrue a corresponding charge of −1e, instead there is a significant amount of electronic density in the intermediate region,⁴⁴ indicative of a partially covalent bond in agreement with experiment.⁴⁶ Nevertheless this suggests that to weaken the Mg–H bonding one should decrease the ionic interaction. Calculations with the PW91 functional on the effects of vacancies reveals that the hydrogen removal energies are strongly reduced with respect to the bulk (0.73 eV/H₂ compared to the value of 2.90 eV/H₂ for bulk MgH₂). However, with respect to dopant inclusion their formation energies are much larger and thus are less likely to form.⁴⁷ DFT-PBE simulations show that vacancies have small activation energies for self-diffusion from the surface to the bulk (0.45 to 0.70 eV), much smaller than the energies for hydrogen release (1.78 to 2.80 eV/H₂) from the surface, implying that surface desorption is more likely to be the rate limiting step than vacancy diffusion.⁴⁸

Fig. 3 Density difference plot (Δρ = ρ(Mg₄₈H₉₆) −ρ(H) −ρ(Mg₄₈H₉₅)) for the bulk MgH₂, where orange represents charge density accumulation and blue represents charge density depletion. Mg atoms are in green, H atoms in purple. Figure taken from ref. 38.

Metal dopant substitution in order to modify Mg–H bonding has also being investigated by DFT, e.g.Cu,^38,42Nb,⁴⁷Zr,⁴⁹Y,⁴⁹Fe,⁵⁰Mn,⁵⁰Ni,^38,42,51Sc,^49,51 and coupled LiAl codopants.⁵² In all these cases the addition of substitutional dopants reduces the ionicity of the Mg–H binding via the formation of covalent bonds with the dopant. In addition, substitutional dopants can aid hydrogenation of magnesium systems by strongly reducing the activation barrier for H₂ splitting, the rate-limiting step for hydrogenation.³¹ In particular, Ti and V are good catalysts for barrierless H₂ splitting (as found by the NEB method) on Mg(0001).^53–55 Unfortunately they both trap hydrogen near the active site, inhibiting diffusion into the bulk. A DFT-PBE simulation comparing Ni and Ti substitutional dopants in Mg(0001) shows that, while Ni is worse for H₂ splitting due to the presence of a small barrier of height 0.1 eV, it is much better for the subsequent diffusion of hydrogen into the bulk due to the dopant only weakly binding hydrogen compared to the much stronger binding evinced by Ti.⁵⁶ Even with dopant addition MgH₂ does not by itself meet the targets for hydrogen storage applications.⁴⁰ Therefore other solid hydride materials are also being investigated, of which lithium amide (LiNH₂) is of major interest.

3.2 Lithium amide (LiNH₂)

Lithium amide (LiNH₂) has had a long history as a reagent in synthetic organic chemistry,^57,58 however recent research has highlighted its potential for hydrogen storage.⁵⁹ Pure LiNH₂ decomposes to lithium imide (Li₂NH) and NH₃ at temperatures above 573 K, however upon intermixing with lithium hydride (LiH), the mixture produces H₂ at a temperature of 423 K to 723 K.⁵⁹ The reaction is

LiNH₂ + 2LiH ↔ Li₂NH + LiH + H₂↔ Li₃N + 2H₂

where the first reaction liberates 5.5 wt% H₂ and the second 5.2 wt% H₂.^59,60 The proposed mechanism for this lower temperature reaction is a polar interaction involving protic hydrogen (H^δ+) from the amide and hydridic hydrogen (H^δ−) from the hydride.^59,61 However, additional mechanisms involving Li⁺⁶² and NH₃ diffusion^63,64 have also been suggested. The addition of 1 mol% of nanosized Ti or TiCl₃ can reduce the range of hydrogen desorption to 100 K, although without altering the onset temperature.⁶⁵

While LiNH₂–LiH mixtures are of interest due to their high hydrogen content they unfortunately possess several problems that affect their applicability for room temperature applications. These include: (a) high dehydrogenation temperature of 473 K required for significant H₂ plateau pressures, (b) slow release kinetics on the order of several days,⁶² and (c) lack of cyclability.⁶⁶ This final issue is believed to be related to ammonia release from LiNH₂ during reaction with LiH.⁶³

The crystal structure of LiNH₂ as found by powder neutron diffraction is a body centred tetragonal with lattice constants a = 5.034 and c = 10.256 Å, see Fig. 4.⁶⁷ In order to clarify the thermodynamic properties of LiNH₂ and Li₂NH DFT simulations have played a significant role in determining the nature of the bonding in these systems,^68–70 and it was unambiguously found that GGAs are substantially better in matching the experimental enthalpy of formation than the LDA (173.1 kJ mol⁻¹ (experiment), 176 kJ mol⁻¹ (GGA) and 245.8 kJ mol⁻¹ (LDA) for each formula unit).⁶⁸ Electronic structure analysis indicates that the Li atoms do not contribute to the occupied states but have a strong overlap with the unoccupied states, whereas the N and H atoms do have strong contributions to the occupied states. This implies that the Li atoms are ionized to form cations while the NH₂groups form anions, in agreement with experiment.⁶⁷ The bonding between N and H atoms was found to be covalent,⁶⁹ which explains the short experimental N–H distances of 0.94 and 0.98 Å. Confirmation of the nature of the binding between Li and NH₂ ions, and the covalent nature of the N–H bond is provided by a calculation of the electron localisation function (a method of determining electron pair probability in multielectron systems, and thus types of bonds)⁷¹ using the PW91 functional, see Fig. 5.⁶⁸ The unequal bonding between N and H atoms (as found by simulations using the PBE functional) lead to two almost equivalent transient steps upon dehydrogenation: (1) [LiNH]⁻ + H⁺, and (2) Li⁺ + [NH₂]⁻.⁷⁰

Fig. 4 Crystal structure of LiNH₂. Red (large), green (middle), and blue (small) spheres represent Li, N, and H atoms, respectively. Figure taken with permission from ref. 69. Copyright (2005) by the American Physical Society (URL: http://link.aps.org/abstract/PRB/v71/e195109).

Fig. 5 Electron localisation function calculated for LiNH₂ in a plane containing an (NH₂)⁻ complex. Figure taken with permission from ref. 68. Copyright (2005) by the American Physical Society (URL: http://link.aps.org/abstract/PRB/v72/e125120).

Several different crystal structures for Li₂NH have been proposed, see Table 1. Balogh et al.⁷⁷ showed using DFT-PW91 simulations that the disordered models experimentally found by Noritake et al.⁷⁵ and Ohoyama et al.⁷⁶ are equivalent to lower symmetry fully occupied orthorhombic and rhombohedral structures, respectively. Low temperature Pnma⁷⁸ (found with DFT-PW91) and Pbca⁷⁹ (found with DFT-PBE) structures have been predicted by simulations but have not been experimentally confirmed. Recent DFT-PBE simulations find that the lowest energy structure is the cell of symmetry Pbca, however in that work the authors themselves ascribe this to intrinsic errors in DFT and favour the cell of symmetry Ima2.⁷² We suggest that the reason for this conflict is due to the nature of atomic binding in the crystal. As for LiNH₂, the bonding between Li and NH groups is ionic.^73,77 However, there are significant quantum effects in the N–H binding.⁷⁴ The hydrogens are delocalized and thus able to reside over several sites. Thus there is still considerable disagreement about the crystal structure of Li₂NH.

Table 1 Different proposed crystal structures for the Li₂NH system

Symmetry	Lattice parameters/Å	Ref.
Fm3̄m	5.077	75
F4̄3m	5.077	76
Ima2 (T < 400 K)	10.087	77
Fm3̄m (T > 400 K)	5.092	77
Pnma	7.733 (a), 3.600 (b), 4.872 (c)	78
Pbca	5.120 (a), 10.510 (b), 5.270 (c)	79

---

As the binding for both amide and imide systems is ionic an obvious method of decreasing the bond strength is to decrease the charge on the cations, e.g. replacing Li with Mg. The effects of Mg substitutional doping of Li for LiNH₂ supercells were investigated using DFT-PW91 simulations by Jin et al.⁸⁰ They found that the Mg dopant has a significant effect, reducing the Mulliken population overlap in the N–H moiety implying a weaker bond and thus lower temperature H₂ release, in qualitative agreement with experiment.⁸¹ Jin et al. also investigated the electronic structures of several other doped lithium amides, finding that the system with the least population overlap in the N–H moiety is Al-doped LiNH₂. DFT-PBE simulations by Zhang et al. on the effects of Mg and K substitution for Li find that the Mg dopant has a weaker ionic bond with the NH₂group, but also strengthens the bonds of the surrounding Li⁺ ions with their associated [NH₂]⁻groups, suggesting that the stronger Li–NH₂ bonds are the key reason behind the efficacy of Mg doping.⁸² However the energy required to form the Mg dopant is quite substantial (∼1.8 eV) and therefore Zhang et al. suggest that P substitution of N may further reduce hydrogen binding energy. These calculations are of low concentration dopants, as of yet there have been no simulations of the higher concentration (Li,Mg)NH₂ crystal structures, such as the recently determined orthorhombic Li₂Mg(NH₂)₂ crystal.⁸³ Despite controversy behind the reaction mechanism of the combined LiNH₂ + LiH system, no theoretical research has been performed on the mobility of suggested species in the bulk to date, which we regard as an oversight.

3.3 Lithium borohydride (LiBH₄)

Lithium borohydride (LiBH₄) was first synthesized in 1940 and has a history in organic synthesis as a reducing agent. Recently, due to its large theoretical hydrogen storage capacity of ∼18.3 wt% it has become of interest for hydrogen storage applications. Experimentally, lithium borohydride reversibly desorbs hydrogen in two steps, following melting at 560 K

2LiBH₄→ 2LiBH + 3H₂→ 2LiH + 2B + 3H₂

releasing 13.1 wt% hydrogen.^60,84 As this is a multistep process with the associated low kinetics, and as this occurs at heightened temperatures, unmodified LiBH₄ is not ideal for applications.

In comparison with other hydrogen storage materials such as MgH₂ and LiNH₂, the phase diagram of LiBH₄ is more complicated, with a transition from one crystal structure to another occurring at 381 K.⁸⁵ The room temperature crystal structure is uncontroversial, being an orthorhombic unit cell with the space groupPnma and lattice vectors a = 7.179 Å, b = 4.437 Å, and c = 6.803 Å, see Fig. 6.^85,86 Simulations using both LDA and PW91 functionals find that at zero temperature PW91 more accurately reproduces the experimental lattice vectors.⁸⁷

In contrast, there is a significant controversy over the crystal structure of the high temperature phase. Initial experimental results indicated that the high temperature phase is hexagonal with space groupP6₃mc and lattice vectors a = 4.276 Å, and c = 6.948 Å, see Fig. 6.⁸⁵ However, simulation does not find that this is the case. Calculations of the hexagonal phase using the PBE functional find that the lowest energy structure does not match experiment, with the final theoretical crystal structure found to be unstable (e.g. possessing phonon modes with negative frequencies).⁸⁸ Calculations of the Helmholtz free energy using the PW91 functional find that the preferred high temperature crystal structure is a cell analogous to LiAlH₄ of symmetry P2₁/c and lattice vectors a = 7.267 Å, b = 7.174 Å, c = 7.683 Å and β = 148.52°, with the experimental hexagonal structure being vibrationally unstable.⁸⁹ The stability of the P2₁/c phase over the P6₃mc phase using both LDA and PW91 functionals was also verified by Frankcombe et al.⁹⁰ Clearly the description of the binding is significantly different for the PBE and PW91 functionals, with the latter finding that the P6₃mc is a local minimum on the potential energy surface while the former is not. Indeed, recent high quality X-ray diffraction results have reverified the initial hexagonal P6₃mc phase.⁹¹ Two principle factors have been suggested for the failure of DFT to accurately model the high temperature phase. Firstly, the problem of DFT in describing the high temperature structure of LiBH₄ has been related to the anharmonicity of low frequency phonon modes corresponding to the rotation of molecular BH₄ units.^91,92 Secondly, the high temperature phase is found to exhibit significant disorder in the lattice positions of the hydrogens.⁹³ This causes significant issues for simulating these materials in a periodic supercell, as in order to accurately model this orientationally disordered material would require the explicit construction of very large supercells so as to capture the full disorder of the system. This is prohibitively expensive to simulate and thus currently beyond the capability of DFT-based methods. In light of these issues the theoretically proposed Cc phase predicted to form at 520 K should be regarded as highly tentative in nature.⁸⁸

Fig. 6 Crystal structures of LiBH₄ in (a) orthorhombic room temperature phase, (b) high temperature hexagonal phase. Red represents Li, green B, and blue H, respectively. Figure taken with permission from ref. 88. Copyright (2004) by the American Physical Society (URL: http://link.aps.org/abstract/PRB/v69/e245120).

The nature of the binding between Li and BH₄ subunits is similar to that of LiNH₂, with the binding between Li cations and BH₄ anions ionic and the B–H binding covalent in nature.^88,94 Charge density analysis indicates that the Li ion has a charge of +0.68e, while the B ion has a negative charge of −0.51e.⁹⁵ This has significant implications for the stability of the BH₄ subunit. Du et al. performed simulations using the BLYP functional on isolated BH₄ ions, finding that for an ionic charge of −1 the complex is tetrahedral and stable, but with decreasing negative charge the tetrahedron distorts and loses stability, releasing H₂ at a charge of −0.2e.⁹⁶ Although experiment cannot readily identify H sublattice positions, simulations do indicate the presence of distorted BH₄ subunits, with B–H bond variance ranging from 2%⁸⁸ to 23%.⁹⁰

The experimental enthalpy change in LiBH₄ upon dehydrogenation to the hydride is 68.9 kJ mol⁻¹.⁹⁷ The LDA, PW91 and PBE exchange correlation functionals do not accurately reproduce this formation enthalpy, with zero point energy corrected enthalpies of 57 kJ mol⁻¹ (LDA),⁹⁰ 59 kJ mol⁻¹ (PW91),⁸⁷ and 56 kJ mol⁻¹ (PBE).⁸⁸ The addition of finite temperature corrections within the quasiharmonic approximation does improve the calculated heat of formation, but not to chemical accuracy.⁸⁷ A recently developed method whereby the anharmonic corrections to the free energy are estimated from fitting to DFT molecular dynamics calculations at high temperature promises to further improve the determination of thermodynamic properties.⁹²

Simulations of different surfaces of the LiBH₄ crystal indicate that the (010), (100) and (101) faces all have low surface energies of approximately 0.12 J m⁻² and thus should easily be formed from ball milling.⁹⁵ These faces have hydrogen atom vacancy formation energies of 180 to 200 kJ mol⁻¹. We note that as the surface as a whole will be less polarised than the bulk then the BH₄ units possess less negative charge, thus their hydrogen binding energies will also be significantly reduced with respect to the bulk. An alternative way to reduce the charge on the BH₄ unit is to create a neutral Li vacancy. Du et al. performed DFT-PBE simulations on the (010) surface of LiBH₄ to quantify the effects of Li vacancies on hydrogen storage properties.⁹⁸ They found that the hydrogen removal energy from the ideal surface was 1.35 eV/H₂, however in the presence of mono- and di-vacancies this is reduced to 0.08 eV/H₂ and −0.45 eV/H₂, respectively. Vacancies also strongly reduce the activation barriers for hydrogen release (from 3.64 eV for the ideal system to 0.23 eV for the di-vacancy), indicating that the creation of lithium vacancies is essential in order to increase kinetics.

In order to utilise LiBH₄ as a hydrogen storage materials it is necessary to reduce the temperature for hydrogen release. Therefore it is of interest that there is an empirical correlation between increasing Pauling electronegativity of the cation and decreasing temperature of desorption for the borohydride compounds.⁹⁹ Theoretical simulations on the effects of substitutional copper doping into lithium borohydride, where copper has a larger electronegativity than Li, to form (Li_1−xCu_x)BH₄, where x = 0.25, 0.5, 0.75, confirm this. It is found that the heat of formation decreases with increasing copper content meeting application requirements at x = 0.3.¹⁰⁰ Surprisingly, partial substitution of the H atoms by F atoms in the BH₄ anion has also been shown to reduce decomposition enthalpies.¹⁰¹ However, this “functional anion” concept has not been experimentally verified as of yet.

3.4 Thermodynamic destabilisation and design of reaction paths

As seen above, LiNH₂ and LiBH₄ are both materials that by themselves are not suitable hydrogen storage materials but in combination with other materials such as LiH or MgH₂ become more suitable. These are examples of thermodynamic destabilization. The work of Reilly and Wiswall on alloys of magnesium and nickel, and magnesium and copper showed that it is possible to modify the thermodynamics of hydrogenation or dehydrogenation reactions by using additives to form compounds or alloys in the dehydrogenated state that are energetically favourable with respect to the products of the reaction without additives.¹⁰² This reduces the enthalpy and increases the equilibrium pressure. Chen et al. have shown that the addition of LiNH₂ to LiH

LiNH₂ + LiH ↔ Li₂NH + H₂

liberates 6.5 wt% H₂ at lower temperatures.⁵⁹ Vajo et al. have shown that the addition of MgH₂ to LiBH₄

2LiBH₄ + MgH₂↔ MgB₂ + 2LiH + 4H₂

reversibly liberates >9 wt% H₂ at approximately 545 K.¹⁰³ Clearly thermodynamic destabilization is an interesting and fruitful concept, what can simulation do to aid in the discovery of new reactions?

The full calculation of the enthalpy of formation and Gibbs free energy, as discussed in section 2.2, is computationally intensive and laborious. However approximations can be made to speed up calculations, allowing the possibility of rapid pre-screening of selected reactions. Instead of calculating the Gibbs free energy, the enthalpy of the reaction can be calculated instead, thus allowing an estimate of the temperature of the dehydrogenation reaction through ΔH = TΔS. The entropy change can be estimated to be of the order of the entropy of H₂ gas at standard conditions; in practice this is estimated to be 95 ≤ΔS≤ 140 J K⁻¹ mol⁻¹.¹⁰⁴ Therefore to achieve dehydrogenation temperatures between 323 and 423 K requires an enthalpy change of 30 ≤ΔH≤ 60 kJ mol⁻¹. In addition, DFT typically has an error of approximately 10 kJ mol⁻¹, therefore calculations should search for reactions with an enthalpy change of 15 ≤ΔH≤ 75 kJ mol⁻¹.

This approach was initially used by Alapati et al. using DFT-PW91 to calculate the enthalpies of 49 compounds and the enthalpy change involved in 133 reactions based on these compounds.¹⁰⁴ They found five promising reaction schemes. The Gibbs free energy of these desirable reactions was then calculated in order to determine the validity of the pre-screening approach.¹⁰⁵ It was found that the enthalpy change ΔH for these reactions can be estimated using the formula

ΔH≈ΔU(T = 0) + δE

where ΔU(T = 0) is the difference in electronic energies only, and δE is a error factor that is in the range −20 < δE < 0 kJ mol⁻¹. The small empirical variance in δE reflects the fact that the complex metal hydrides share similar bonding motifs.¹⁰⁶ This new expression reduces the time needed to compare enthalpies, with Alapati et al. then using this approach to screen 340 reactions between 163 compounds, finding 12 reactions of interest, although only 6 reactions have significant hydrogen plateau pressures.¹⁰⁷ Recently a new approach has been developed where the grand canonical potential of a mixture of compounds is determined as a function of temperature and pressure (where the only chemical potential that varies as a function of temperature and pressure is that of hydrogen), with the mixture with the lowest potential being favoured. By repeating this calculation at various temperatures, reactions can be detected as changes in the equilibrium state. This approach has been used to evaluate 16 million reactions, from a 14 element set that comprises a set of 212 compounds, finding 43 single step dehydrogenation reactions of interest.¹⁰⁸ We tabulate all reactions of interest from these pre-screening simulations in Table 2. We note that the majority of high hydrogen content paths involve the (BH₄) ligand, either in LiBH₄ or Mg(BH₄)₂ form. Research into the crystal structure, binding, and thermodynamics of these systems is of immense interest and is an ongoing problem.⁹¹

Table 2 List of identified single-step reactions which yield ≥6.5 wt% H₂ and have electronic binding energies in the range of 15 ≤ΔH≤ 75 kJ mol⁻¹. The Table is amalgamated from the calculations of Alapati et al.,^107,108 and Siegel et al.¹⁰⁹ For Li₅N₃Si the unit cell is too large for convenient calculation, and is instead approximated by smaller unit cells with varying stoichiometry. This is indicated by ^a. Calculations of LiBH₄ are taken with respect to the room temperature orthorhombic cell, see Alapati et al. for more information on how this affects enthalpies at heightened temperatures¹⁰⁷

Reaction	Wt% H2	ΔU/kJ mol^−1 H2	Ref.
Mg(BH4)2→ MgB2 + 4H2	14.9	54.0	108
Si + 4Mg(BH4)2→ Mg2Si + 2MgB4 + 16H2	13.2	52.7	108
7LiBH4 + MgH2→ 7LiH + MgB7 + 11.5H2	13.0	71.5–75.5	107
4LiBH4 + MgH2→ 4LiH + MgB4 + 7H2	12.5	69.2	107
LiBH4 + C → LiBC + 2H2	12.0	45.1	107
LiBH4 + 2LiNH2→ Li3BN2 + 4H2	11.9	24.3	107
6LiBH4 + CaH2→ 6LiH + CaB6 + 10H2	11.7	62.1	108
8LiBH4 + MgH2 + BN → 8LiH + MgB9N + 13H2	11.6	66.3	108
MgH2 + 2LiBH4→ 2LiH + MgB2 + 4H2	11.6	66.8	107
3Si + MgSiN2 + 9Mg(BH4)2→ 4Mg2Si + 2MgB9N + 36H2	11.2	48.6	108
BN + 4Mg(BH4)2→ 3MgH2 + MgB9N + 13H2	10.9	51.2	108
NaH + 2Mg(BH4)2→ NaMgH3 + MgB4 + 7H2	10.67	53.2	108
CaH2 + 1.5Si + 3Mg(BH4)2→ CaB6 + 1.5Mg2Si + 13H2	10.6	45.4	108
2C + Mg(BH4)2→ MgB2C2 + 4H2	10.3	43.1	108
CaH2 + 3Mg(BH4)2→ 3MgH2 + CaB6 + 10H2	9.9	47.5	108
8LiBH4 + Mg2Si → 8LiH + Si + 2MgB4 + 12H2	9.6	74.0	108
2LiBH4 + ScH2→ 2LiH + ScB2 + 4H2	8.9	49.7	107
2LiBH4 + TiH2→ 2LiH + TiB2 + 4H2	8.6	22.2	107
2LiBH4 + NaMgH3→ 2LiH + NaH + MgB2 + 4H2	8.6	68.9	108
MgH2 + 2LiBH4 + 2C → MgB2C2 + 2LiH + 4H2	8.6	55.3	107
3NaH + BN + 4Mg(BH4)2→ 3NaMgH3 + MgB9N + 13H2	8.4	48.8	108
2LiBH4 + VH2→ VB2 + 2LiH + 4H2	8.4	24.7	109
LiNH2 + MgH2→ LiMgN + 2H2	8.2	31.9	107
2LiBH4 + Mg(NH2)2→ MgH2 + 2LiH + 2BN + 4H2	8.1	20.6	108
ScH2 + Mg(BH4)2→ MgH2 + ScB2 + 4H2	8.0	37.5	108
ScH2 + 2LiBH4 + C → ScB2C + 2LiH + 4H2	7.9	52.9	107
MgH2→ Mg + H2	7.7	64.7	108
5B + Mg(BH4)2→ MgB7 + 4H2	7.5	41.5	108
2MgH2 + Mg(NH2)2→ Mg3N2 + 4H2	7.4	26.0	108
CaH2 + 2LiBH4 + 2C → CaB2C2 + 2LiH + 4H2	7.4	59.7	107
CaH2 + 3NaH + 3Mg(BH4)2→ 3NaMgH3 + CaB6 + 10H2	7.3	44.3	108
6LiBH4 + 2ScN → 6LiH + 2ScB2 + 2BN + 9H2	7.3	59.5	108
3Si + 8BN + 5Mg(BH4)2→ 2MgB9N + 3MgSiN2 + 20H2	7.3	47.0	108
4LiBH4 + K2MgH4→ 4LiH + MgB4 + 2KH + 7H2	7.3	74.8	108
4LiH + 3Mg(NH2)2 + 2C → 2Li2CN2 + Mg3N2 + 8H2	7.2	47.8	108
3LiNH2 + 2LiH + Si → Li5N3Si + 4H2	7.2	34.2/23.3^a	107
3Ca(AlH4)2 + 4LiH → 4AlLi + Al2Ca + 2CaH2 + 12H2	7.2	31.8	107
6LiBH4 + 2TiN → 6LiH + 2TiB2 + 2BN + 9H2	7.1	35.9	108
2LiNH2 + C → Li2CN2 + 2H2	7.0	31.4	107
ScH2 + 2LiBH4 + 2C → ScB2C2 + 2LiH + 4H2	7.0	52.4	107
Ca(BH4)2 + ScH2→ CaH2 + ScB2 + 4H2	6.9	44.8	109
Ca(AlH4)2 + 2LiH → 2AlLi + CaH2 + 4H2	6.9	33.3	107
Al + MgB9N + 2.5Mg(BH4)2→ AlN + 3.5MgB4 + 10H2	6.8	53.6	108
2LiBH4 + MgB2→ 2LiH + MgB4 + 3H2	6.8	72.5	108
MgB7 + 1.5Mg(BH4)2→ 2.5MgB4 + 6H2	6.7	50.2	108
12LiH + 3Mg(NH2)2 + 4BN → 4Li3BN2 + Mg3N2 + 12H2	6.7	54.2	108
3Ca(AlH4)2 + 2Si → 2Al2Ca + Al2CaSi2 + 12H2	6.7	24.1	107
3Ca(AlH4)2 + 2Si → 4Al + Al2Ca3Si2 + 12H2	6.7	28.4	107
6Ca(AlH4)2 + 17MgH2→ Al12Mg17 + 6CaH2 + 35H2	6.7	35.4	107
Ca(BH4)2 + TiH2→ CaH2 + TiB2 + 4H2	6.7	17.4	109
K2MgH4 + 2Mg(BH4)2→ MgB4 + 2KMgH3 + 7H2	6.6	51.2	108
Ca(BH4)2 + VH2→ CaH2 + VB2 + 4H2	6.6	20.8	109
28LiH + 9Mg(NH2)2 + 4VN → 4Li7N4V + 3Mg3N2 + 32H2	6.5	47.5	108
2ScN + 3Mg(BH4)2→ 3MgH2 + 2ScB2 + 2BN + 9H2	6.5	43.1	108
NaH + ScH2 + Mg(BH4)2→ NaMgH3 + ScB2 + 4H2	6.5	34.8	108

---

General thermodynamic guidelines have been developed from density functional theory simulations to aid in the prediction of destabilisation reactions.¹⁰⁹ In brief:

1. The enthalpy of the proposed destabilization reaction must be less than the decomposition enthalpies of the individual reactant phases.

2. If the proposed reaction involves a reactant that can adsorb hydrogen (such as an elemental metal), the formation enthalpy of the corresponding hydride cannot be greater in magnitude than the enthalpy of the destabilized reaction.

3. In general, it is not possible to tune the thermodynamics of destabilized reactions by adjusting the molar fractions of the reactants. There is only one stoichiometry corresponding to a single step reaction with the lowest possible enthalpy; all other stoichiometries will release H₂ in multistep reactions, where the initial reaction is given by the lowest-enthalpy reaction.

The rapid computational pre-screening of reactions that lead to thermodynamic destabilization is of immense scientific interest. However we emphasize that this is only an initial determination of likely reaction chemistry and suitability for hydrogen storage applications. Firstly, neglecting all entropy contributions apart from that of H₂ gas may strongly affect properties, e.g. vibrational entropy contributions alter the solubility of Al in Si by an order of magnitude.¹¹⁰ Secondly, as mentioned in the general guidelines above, reactions involving almost unstable compounds (e.g. Ca(AlH₄)₂) may lead to direct decomposition into different products.¹¹¹ Thirdly, only thermodynamics are determined not kinetics. For example, Ca(AlH₄)₂ has only a ΔU of 14 kJ mol⁻¹ but is kinetically stable for temperatures up to 473 K.¹⁰⁷ Fourthly, thermodynamics of reactions can only be estimated for known crystal structures, if unknown crystal structures exist for possible intermediate compounds this will affect reaction enthalpies. Fifthly, errors in DFT of the order of 10 kJ mol⁻¹ will lead to errors in the plateau pressure upon dehydrogenation of three orders of magnitude.¹⁰⁷ Regardless, simulation of thermodynamic destabilization is an interesting concept that potentially may prove very fruitful for hydrogen storage research.

4. Carbon systems

Carbon systems (nanotubes, fullerenes, and nanoporous structures) are another group of widely-researched hydrogen storage materials.¹¹² Initial interest was spurred by the early results of Dillon et al. who observed a hydrogen uptake between 5 and 10 wt% for single-walled nanotubes at low temperature (133 K) and very low pressures (0.04 MPa).¹¹³ This was a spectacular result, and along with other spectacular results for other carbon nanostructures¹¹⁴ indicated a promising future for hydrogen storage. Later experiments however revealed much more modest storage capacities of ∼1.55 wt% at 573 K and 0.01 MPa.¹¹⁵ This is due to the weak nature of the Van der Waals interaction between H₂ and carbon. Simulation can play a significant role in understanding and providing insight into two types of carbon systems: graphite sheets which have a sheet structure with a repeating hexagonal motif, and single walled carbon nanotubes (CNTs) which can be envisaged as formed by the rolling of a single graphene sheet.

4.1 Interaction of carbon systems with hydrogen

There have been many simulations to understand the hydrogen adsorption process on graphite systems so as to predict optimal structures, involving both empirical potentials,^116–118 and quantum mechanical calculations.^119–126Hydrogen atoms can bind to carbon nanostructures with a binding energy of magnitude 1 eV within the LDA approximation.¹²⁷ However atomic adsorption is not energetically preferred to molecular adsorption, where DFT finds that molecular binding energy is of the order of 0.1 eV/H₂,^{119,124–126} this is an overestimate with respect to experiment which finds a binding energy of 0.043 eV/H₂.¹²⁸ This is far too weak for applications. A recent calculation that uses highly accurate Møller–Plesset wavefunction-based techniques shows that optimised graphite plates may meet DOE targets if the interlayer distance of the graphite is increased to between 6 and 7 Å.¹²⁹ No mechanism for enabling such a consistent layer separation has been developed however.

Carbon nanotubes have also been investigated for hydrogen storage applications, as due to the curvature of the graphene sheets and varying metallic/semiconducting properties they are expected to present different chemical behaviour to flat graphene sheets. Indeed, the magnitude of the adsorption energy increases with decreasing tube diameter due to the mixing of sp²–sp³ rehybridisation induced by strong curvatures, so a tube of diameter 1.2 nm has a binding energy of 5.9 kJ mol⁻¹/H₂ (0.061 eV/H₂) while a tube of diameter 0.8 nm has a binding energy of 8 to 9 kJ mol⁻¹/H₂ (0.0829 to 0.0933 eV/H₂).^130,131Hydrogen addition to the exterior of a SWNT will be in the form of physisorption with a binding energy of approximately 0.1 eV/H₂, as there is a large barrier of approximately 3 eV in both DFT-LDA and DFT-PW91 simulations for H₂ to dissociate into atomic hydrogen,^127,132 with chirality not affecting binding energy.¹³³ Small diameter nanotubes can sustain 100% coverage but nanotubes over 1 nm in diameter can sustain at best 80% coverage.¹³⁴

DFT-B3LYP simulation shows that excessive atomic hydrogen adsorption can lead to the partial breakdown of tube structure, with a greater than 50% coverage on a (10,0) nanotube creating pores that will allow hydrogen access to the interior of the nanotube.¹³⁵ We emphasize that without a source of readily available hydrogen atoms this particular set of behaviours is very unlikely to happen, as thermodynamically H₂ prefers to adsorb in molecular form. However if the nanotubes are pretreated with an atomic hydrogen beam then a mixture of reversible chemisorption and physisorption interactions has been observed, with a storage capacity of 5.1 ± 1.2 wt%.¹³⁶

Simulation tends to model nanotubes in isolation. However in reality carbon nanotubes tend to exist in bundles with multiple sites for hydrogen binding. These may be mechanically manipulated. DFT has been used to investigate the effects of deformation on reactivity, finding that reactivity increases with increasing bend angle.¹³⁷ Simulations of the solid phase within the PW91 functional indicate that neighbouring carbon nanotubes under a compressive stress of 12.5% dissociates H₂ present in a channel of the interior of the nanotube bundle, see Fig. 7.¹³⁸ This process is reversible upon the release of pressure with an energy barrier of ∼2 eV. Such tuning of adsorption and storage properties are highly promising avenues for future experimental research, with such processes also changing the magnetic properties of the nanotube bundle.¹²⁵

Fig. 7 The structures and transition structures for the dissociative H₂ chemisorption on an array of carbon nanotubes in a solid under high pressure. Three steps are shown: H₂ adsorption and dissociation on a single CNT (A), deformation and displacement of a single H atom through transition state (TS) A so that it chemisorbs on an adjacent CNT (B), rotation of the initial CNT, and rehydrogenation of this CNT through TS B such that final product is a 1,2 addition structure. This causes a noticeable deformation of the CNT. Figure taken with permission from ref. 138. Copyright (2001) by the American Physical Society (URL: http://prola.aps.org/abstract/PRL/v87/i20/e205502).

Although in general GGAs are better functionals than the LDA, an exception occurs for the interaction of H₂ with carbon systems. Comparison of DFT with extremely accurate wavefunction-based methods (Møller–Plesset and diffusion Quantum Monte Carlo) finds that the LDA is more accurate.^120,139 This is largely fortuitous and due to favourable error calculations. Recently a dispersion force corrected density functional theory has been introduced.¹⁴⁰ In this treatment a dispersion energy term to represent the Van der Waals energy is added to the total energy of the system (E_Tot = E_DFT + E_dis) where

where R_ij is distance of the H₂ molecule from the carbon plane, f(R_ij) is a damping function which becomes unity for large R_ij and zero for small R_ij. The C₆ coefficients for C and H are 19.1 and 2.8 Hartrees ×a⁶₀ and are derived from experimental atomic polarisabilities. For diatomic coefficients C^ij₆ can be obtained by the Slater–Kirkwood combination rule.¹⁴¹ This correction gives binding energies and equilibrium distances for H₂ on graphite that are much closer to experiment. Nevertheless most simulations on H₂ adsorption on carbon do not use this methodology.

4.2 Interaction of alkali-metal doped carbon systems with hydrogen

By themselves carbon-based hydrogen storage materials will not meet DOE targets for hydrogen adsorption. This is due to the strength of the Van der Waals interaction between H₂ and the carbon support. To strengthen the binding requires the addition of electrostatic interactions and/or orbital interactions.¹⁴² The positive quadrupole moment of H₂ cannot be modified, therefore in order to increase adsorption energies one needs to charge the carbon system.¹⁴³ In order to cause this charging, the carbon systems could be doped with highly electropositive Group I and II metal nanoparticles, facilitating charge transfer. Simulations show that the insertion of Li and Rb,¹⁴⁴ Cs,¹⁴⁵K,¹⁴⁶ Ca and Sr,¹⁴⁷ and Ni¹⁴⁸ dopants induces a charge on the nanotubes. Alkali metal-doped systems in particular display marked increases in hydrogen storage capacities.¹⁴⁹ Both classical Monte Carlo and DFT-LDA simulations show that Li⁺ ions act as acidic cores for hydrogen binding, see Fig. 8.^150,151 The insertion of K or Cs has been measured to lead to an excess negative charge of −1e per 8 carbon atoms.¹⁴⁵ Recent DFT-PW91 simulations indicate that Ca and Sr dopants can lead to a charge transfer of 2e per dopant atom to a C₆₀fullerene.¹⁴⁷

Fig. 8 Effects of lithium intercalation into graphite, with associated electron-density redistribution for the adsorption of Li on graphene. (a) Top view of surface of constant Δρ = ρ(graphene + Li) −ρ(graphene) −ρ(Li). The red and orange surfaces correspond to +2.5 and −2.5 × 10⁻²e a.u.⁻³, respectively. (b) Side view of constant Δρ. Figure is based on plane wave calculations performed by the authors based on the results and Figure of ref. 133.

Dopants therefore directly boost the H₂ binding energy of the carbon substrate, but in addition they can act as H₂ binding centres. Simulation has aided in the design of these materials,^152–154 however the majority of first-principles simulations consider only the interaction with atomic hydrogen, of limited use for realistic systems where the thermodynamic reference for hydrogen should be the molecule. Froudakis performed a full DFT-B3LYP/MM simulation of molecular adsorption on atomic K dopants adsorbed on a (5,5) SWNT.¹⁵⁴ They found that each K dopant physisorbs 3H₂, with decreasing adsorption energies of 3.4 to 1.8 kcal mol⁻¹/H₂ (0.147 to 0.078 eV/H₂), smaller than the 11.1 kcal mol⁻¹ (0.481 eV/H₂) found by empirical calculations. Positively charged Li ions complexed to a substrate will bind H₂ with binding energies of magnitude 12 kJ mol⁻¹/H₂ (0.124 eV/H₂).¹⁴² Yoon et al. found that a Ca₃₂C₆₀ structure will bind >8.4 wt%, with each dopant binding 5H₂ molecules at a consistent binding energy of 0.4 eV/H₂.¹⁴⁷ In general however direct H₂ binding to an atomic dopant will involve orbital interactions, e.g. H₂σ-bond overlap with a LUMO state or σ*-bond overlap with a HOMO state.¹⁴²

4.3 Interaction of transition-metal doped carbon systems with hydrogen

Recently a new approach has been proposed, transition metal (TM) atom doping.^155,156 This method was initially developed in analogy with organometallic cyclopentadiene (C₅H₅) complexes, as TM atoms interact with a fivefold ring of fullerene C₆₀ in a similar way to cyclopentadiene (Cp) rings, forming a Dewar coordination.¹⁵⁷ This forms a strong bond between the TM atom and the fullerene. Hydrogen molecules then interact with the TM atom via a Kubas interaction, forming σ-bonds with the empty d-states of the TM atom, see Fig. 9.¹⁵⁸ The maximum number of hydrogens that bind to a TM atom is determined by the 18-electron rule; once the s-, p-, and d-orbitals of the TM atom are filled it can no longer bind H₂. The maximum number of hydrogens bound is

N_max = 18 −n_v− 5

where n_v is the number of valence electrons of the TM atom and 5 is the number of electrons contributed by the fivefold ring. Therefore the optimum TM dopant is Sc, as it has the largest number of empty d-states. Significantly, the binding energies for these organometallic systems are in the exact range needed for hydrogenstorage applications.

Fig. 9 Schematic energy diagram and structure model for hybridisation of (a) Cp ring with Sc forming CpSc, (b) CpSc with 2H forming Cp[ScH₂], and (c) Cp[ScH₂] with H₂ forming Cp[Sc(H₂)H₂]. In (a), only the last two π states are shown for Cp. In (c) the two upper states of Cp[ScH₂] are neglected. Arrows indicate level occupation. Figure taken with permission from ref. 155. Copyright (2005) by the American Physical Society (URL: http://link.aps.org/abstract/PRL/v94/e155504).

Using this rational design methodology Zhao et al. found (from DFT-PW91 simulations) that a C₆₀fullerene with each fivefold face decorated by a Sc atom can reversibly adsorb 48 H₂ molecules.¹⁵⁵Boron-doping the fullerene allows an additional withdrawal of charge from the TM atom so that a C₄₈B₁₂ fullerene has optimal hydrogen storage properties, able to form a C₄₈B₁₂[ScH(H₂)₅]₁₂ complex that has a hydrogen storage capacity of 8.77 wt%, is stable at room temperature and binds each H₂ ligand with an adsorption energy of approximately 0.3 eV/H₂. Concurrently another first-principles simulation (DFT-PBE) demonstrated that a single Ti atom dopant adsorbed on an (8,0) SWNT can bind 4H₂ stably at room temperature.¹⁵⁶ If half of the surface hexagons are covered with Ti atoms then the gravimetric density of H₂ is 7.7 wt%, with a binding energy of 0.18 eV/H₂, reduced due to the interaction between Ti atoms. DFT-PBE simulation of a Ni atom dopant on an (8,0) SWNT shows that each atom can bind 5H₂, releasing hydrogen at a temperature of 328 K, allowing in principle a gravimetric storage capacity of ∼10 wt%.¹⁵⁹ Experimentally however nickel nanoparticles on MWNTs only achieve a storage capacity of 2.8 wt%.¹⁵⁹ We note that this general idea is applicable to other hexagonal net systems such as boron nitride (BN), where DFT-PW91 calculations on an idealised system (BN)₃H₆Sc finds a reversible adsorption capacity of 4H₂.¹⁹

Simulations on other transition metal doped carbon-based systems have also been performed. The DFT-PW91 simulation of organic molecules such as C₄H₄, C₅H₅, and C₈H₈ doped with titanium show that these systems can contain up to 9 wt% hydrogen with an adsorption energy of 0.55 eV/H₂.¹⁶⁰ Again maximum hydrogen storage capacity is determined by the 18-electron rule. A DFT-PBE combinatorial search by Lee et al. on metal-decorated polymers found that Ti-decorated cis-polyacetylene is optimal, with gravimetric storage capacities of 7.6 wt% and binding energies from 0.4 to 0.6 eV/H₂ found.²¹ Density functional theory simulations on ethylene complexes reveal that C₂H₄Ti₂ and C₂H₄Zr₂ can both bind 10H₂, with average binding energies of 0.45 eV/H₂ (Ti) and 0.57 eV/H₂ (Zr), respectively and capacities of ∼14 wt% (Ti) and ∼9 wt% (Zr).^161,162 Such complexes have been recently observed experimentally.¹⁶³ Research on transition metal doped polymers has led to additionally research on light-metal doped polymers, for example C₂H₄Li₂ complexes¹⁶² or Li-doped cis-polyacetylene.¹⁶⁴ Due to the lighter mass of the dopant their gravimetric capacities are higher, however the binding energy for H₂ adsorption is less strong than for the TM decorated polymers, and thus these systems would be less suitable for hydrogen storage applications.

A major issue for the polymers is that they are unstable with respect to TM nanoparticle coalescence, consequently storing less hydrogen, see Fig. 10.¹⁶⁵ Indeed, energetically C₂H₄Ti₂ prefers to dimerise to form 2(C₂H₄Ti₂) through the formation of a Ti–Ti bond, reducing gravimetric capacity to ∼7 wt%.^161,164 This process has also been observed for Ti on a C₆₀ surface, hindering hydrogen storage.¹⁶⁶ Therefore in order to maximise hydrogen storage the TM dopants should be pinned strongly to the carbon substrate. One way of doing this is to add dopants to the carbon system such that the binding of the TM atoms is increased, e.g.nitrogen doping of graphenes.¹⁶⁷ Another way would be to firmly emplace the TM atoms in a carbon matrix, either by defecting the support¹⁶⁸ or by forming a metal carbide.^169,170 Simulations on both the former and latter systems find that the maximum hydrogen capacity is less than expected from the 18 electron rule, e.g. for the latter system, the Ti metcar Ti₈C₁₂, Ti atoms on the corners adsorb 3H₂ and Ti atoms in the centre of a face adsorb a single H₂, see Fig. 11. These metcars can still adsorb 6 wt% H₂, be stable at room temperature, and have desirable H₂ binding energies. We note though that these are idealized models, as no one has experimentally realized the formation and subsequent dispersal of these metcar systems. However they (and all of the organometallic systems discussed in this section) are new systems which are of great potential utility for energy applications.

Fig. 10 Two configurations of Ti₁₂C₆₀. The configuration on the left-hand side is 24.8 eV higher in energy than the configuration on the right hand-side. Figure adapted from ref. 166.

Fig. 11 Example of a Ti metcar Ti₈C₁₂ as a potential hydrogen storage media, binding 17H₂ to the molecule. Figure taken with permission from ref. 169. Copyright (2006) by Elsevier Publishing.

5. Conclusions

Hydrogen storage is a scientifically interesting but challenging problem for the worldwide research community. DFT simulations are best used to simulate chemical events that occur over a length scale of nanometres and a timescale of picoseconds. They thus are of great utility for modelling a wide range of behaviours in hydrogen storage materials. However, DFT cannot model truly collective long duration phenomena explicitly, e.g. the transformation in crystal structure upon dehydrogenation of MgH₂. In addition DFT has an inherent uncertainty in the total energy of approximately 0.1 eV, although due to the systematic nature of the errors relative energies are often much more accurate. For most of the ground state properties, errors in the evaluation of energy differences can be factored into results (although this tends not to occur in the literature). However, for calculations of the activation barrier for chemical reactions (and thus reaction rate), the exponential dependence on energy results in a magnification of error, rendering calculated plateau pressures unreliable at best. Finally simulation, by its very nature, models idealized systems. Solid materials will consist of multiple grains and structural defects, dopants will not be ideally dispersed uniformly throughout the material, and nanotubes and nanoparticles will not exist in isolation. Crystals (such as LiBH₄) where there is a substantial amount of orientational disorder need large supercells in order to treat these systems realistically, which is prohibitively expensive to simulate.

Nevertheless DFT simulation is an extremely useful tool in the study of this field, providing insight into the chemistry and basic mechanisms behind hydrogenation and dehydrogenation mechanisms and allowing rapid design and prototyping of novel hydrogen storage materials. They provide insight into the nature of chemical bonding in MgH₂, clarifying that the mechanism by which transition metals reduce the dehydrogenation enthalpy is by lessening the ionicity of the mixed covalent–ionic Mg–H bond. Simulations on LiNH₂ indicate that the strength of the Li–N ionic bond affects the strength of the N–H covalent bond, with partial Mg or Al substitution of the Li ion reducing the temperature at which the N–H bond breaks. Simulations on LiBH₄, another complex metal hydride, find that the binding between Li and BH₄ units is ionic in nature, and that partial Cu substitution of the Li ion or F substitution of a H atom reduces the temperature at which hydrogen desorption occurs. However simulations for this system also reveal the shortfalls inherent in straightforward applications of density functional theory, where without highly sophisticated modifications to the free energy based on semiempirical data the experimentally observed high temperature phase is not observed. DFT simulations can also be used to screen for new reaction paths that may lead to very fruitful avenues for experimental research. We note in passing that the majority of high hydrogen content paths involve the (BH₄) ligand. Thus further fundamental research both experimental and theoretical in nature, should be performed to understand the crystal structure, binding, and thermodynamics of the LiBH₄ and Mg(BH₄)₂ systems. Finally DFT simulations have been used to prototype and rapidly develop new concepts for carbon based hydrogen storage materials. Of particular interest in this field is the recent advent of new transition metal doped carbon systems, which possess great promise in meeting technological targets for hydrogen storage. However experimental verification of these materials is, apart from a few cases, lacking. Simulation, as a single tool in the repertoire of chemical science, is therefore best used in conjunction with, not as a replacement for, experiment.

Acknowledgements

This work was supported by the EPSRC SUPERGEN Initiative under UK-SHEC (GR/S26965/01, EP/E040071/1) and by a Platform Grant (GR/S52636/01, EP/E046193/1). The authors wish to also thank Prof. Dario Alfè, Dr Louis Hector, Prof. Puru Jena, Dr Boutheina Kerkeni, Dr Kazutoshi Miwa, Dr Qiang Sun, Dr Scott Woodley, and Dr Yufeng Zhao for their help and assistance in preparing Figures.

References

1. V. M. Vishnyakov, Vacuum, 2006, 80, 1053.
2. D. K. Ross, Vacuum, 2006, 80, 1084.
3. Y. Song, Z. X. Guo and R. Yang, Phys. Rev. B, 2004, 69, 094205.
4. L. Schlapbach and A. Zuttel, Nature, 2001, 414, 353.
5. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133.
6. P. A. M. Dirac, Proc. Cambridge Philos. Soc., 1930, 26, 376.
7. J. P. Perdew, A. Ruzsinszky, J. Tao, V. N. Staroverov, G. E. Scuseria and G. I. Csonka, J. Chem. Phys., 2005, 123, 066201.
8. J. P. Perdew and Y. Wang, Phys. Rev. B, 1986, 33, 8800.
9. A. D. Becke, Phys. Rev. A, 1988, 38, 3098.
10. C. Lee, W. Yang and G. Parr, Phys. Rev. B, 1988, 37, 785.
11. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B, 1992, 46, 6671.
12. J. P. Perdew, K. Burke and M. Enzerhof, Phys. Rev. Lett., 1966, 77, 3865.
13. F. Hermann, J. P. Van Dyke and I. P. Ortenburger, Phys. Rev. Lett., 1969, 22, 807.
14. R. M. Martin, Electronic Structure: Basic Theory and Practical Methods, Cambridge University Press, 2004.
15. J. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401.
16. A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
17. A. D. Becke, J. Chem. Phys., 1997, 107, 8554.
18. V. N. Staroverov, G. E. Scuseria, J. Tao and J. P. Perdew, J. Chem. Phys., 2003, 119, 12129.
19. S. A. Shevlin and Z. X. Guo, Appl. Phys. Lett., 2006, 89, 153104.
20. S. A. Shevlin and Z. X. Guo, Phys. Rev. B, 2007, 76, 024104.
21. H. Lee, W. I. Choi and J. Ihm, Phys. Rev. Lett., 2006, 97, 056104.
22. S. V. Alapati, J. K. Johnson and D. S. Sholl, Phys. Chem. Chem. Phys., 2007, 9, 1438.
23. NIST Chemistry Webook, NIST Standard Reference Database Number 69, ed. P. J. Lindstrom and W. G. Mallard, National Institute of Standards and Technology, Gaithersburg, MD, 2005, 20899, http://webbook.nist.gov.
24. J. J. Vajo, S. L. Skeith and F. Meters, J. Phys. Chem. B, 2004, 108, 13977.
25. Z. Lodziana and T. Vegge, Phys. Rev. Lett., 2004, 93, 145501.
26. O. M. Løvvik, S. M. Opalka, H. W. Brinks and B. C. Hauback, Phys. Rev. B, 2004, 69, 134117.
27. P. Hänggi, P. Talkner and M. Borkovec, Rev. Mod. Phys., 1990, 62, 251.
28. T. Vegge, Phys. Chem. Chem. Phys., 2006, 8, 4853.
29. T. Itoh, G. R. Bell, A. R. Avery, T. S. Jones, B. A. Joyce and D. D. Vvedensky, Phys. Rev. Lett., 1998, 81, 633.
30. G. Boisvert, N. Mousseau and L. J. Lewis, Phys. Rev. B, 1998, 58, 12667.
31. T. Vegge, Phys. Rev. B, 2004, 70, 035412.
32. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901.
33. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978.
34. D. Sheppard, R. Terrell and G. Henkelman, J. Chem. Phys., 2008, 128, 134106.
35. J. L. Bobet, C. Even, Y. Nakamura, E. Akiba and B. Darriet, J. Alloys Compd., 2000, 298, 279.
36. M. Bortz, B. Bertheville, G. Bottger and K. Yvon, J. Alloys Compd., 1999, 287, L4–L6.
37. A. J. Du, S. C. Smith, X. D. Yao and G. Q. Lu, Surf. Sci., 2006, 600, 1854.
38. W. Yang, S. A. Shevlin, C. X. Shang and Z. X. Guo, in preparation.
39. K. Higuchi, K. Yamamoto, H. Kajioka, K. Toiyama, M. Honda, S. Orimo and H. Fujii, J. Alloys Compd., 2002, 330–332, 526–530.
40. Z. X. Guo, C. Shang and K. F. Aguey-Zinsou, J. Eur. Ceram. Soc., 2008, 28, 1467.
41. R. W. P. Wagemans, J. H. van Lenthe, P. E. de Jongh, A. J. van Dillien and K. P. de Jong, J. Am. Chem. Soc., 2005, 127, 16675.
42. Y. Song and Z. X. Guo, Appl. Phys. Lett., 2006, 89, 111911.
43. M. Pozzo and D. Alfè, Phys. Rev. B, 2008, 77, 104103.
44. R. Yu and P. K. Lam, Phys. Rev. B, 1988, 37, 8730.
45. M. Tsuda, W. A. Diño, H. Nakanishi and H. Kasai, J. Phys. Soc. Jpn., 2004, 73, 2628.
46. T. Noritake, M. Aoki, S. Towata, Y. Seno, Y. Hirose, E. Nishibori, M. Takata and M. Sakata, Appl. Phys. Lett., 2002, 81, 2008.
47. S. Li, P. Jena and R. Ahuja, Phys. Rev. B, 2006, 74, 132106.
48. A. J. Du, S. C. Smith and G. Q. Lu, J. Phys. Chem. C, 2007, 111, 8360.
49. D. Chen, Y. M. Wang, L. Chen, S. Liu, C. X. Ma and L. B. Wang, Acta Mater., 2004, 52, 521.
50. M. Tsuda, W. A. Diño, H. Kasai and H. Nakanishi, Appl. Phys. Lett., 2005, 86, 213109.
51. M. Tsuda, W. A. Diño, H. Kasai, H. Nakanishi and H. Aikawa, Thin Solid Films, 2006, 509, 157.
52. J.-J. Liang, J. Alloys Compd., 2007, 446–447, 72.
53. A. J. Du, S. C. Smith, X. D. Yao and G. Q. Lu, J. Phys. Chem. B, 2005, 109, 18037.
54. X. Yao, C. Wu, A. J. Du, G. Q. Lu, H. Cheng, S. C. Smith, J. Zou and Y. He, J. Phys. Chem. B, 2006, 110, 11697.
55. A. J. Du, S. C. Smith, X. D. Yao, C. H. Sun, L. Li and G. Q. Lu, Appl. Phys. Lett., 2008, 92, 163106.
56. M. Pozzo, D. Alfè, A. Amierio, S. French and A. Pratt, J. Chem. Phys., 2008, 128, 094703.
57. J. H. Titherly, J. Chem. Soc., Trans., 1894, 65, 504.
58. I. A. Kaye, J. Am. Chem. Soc., 1949, 71, 2322.
59. P. Chen, Z. T. Xiong, J. Z. Luo, J. Y. Lin and K. L. Tan, J. Phys. Chem. B, 2003, 107, 10967.
60. S. I. Orimo, Y. Nakamori, J. R. Eliseo, A. Züttel and C. M. Jensen, Chem. Rev., 2007, 107, 4111.
61. J. Lu, Z. Z. Fang and H. Y. Sohn, Inorg. Chem., 2006, 45, 8749.
62. K. F. Aguey-Zinsou, J. Yao and Z. X. Guo, J. Phys. Chem. B, 2007, 111, 12531.
63. T. Ichikawa, N. Hanada, S. Isobe, H. Leng and H. Fujii, J. Phys. Chem. B, 2004, 108, 7887.
64. L. L. Shaw, W. Osborn, T. Markmaitree and X. Wan, J. Power Sources, 2008, 177, 500.
65. T. Ichikawa, N. Hanada, S. Isobe, H. Y. Leng and H. Fujii, J. Alloys Compd., 2005, 404, 435.
66. W. F. Luo and K. Stewart, J. Alloys Compd., 2007, 440, 357.
67. J. B. Yang, X. D. Zhou, Q. Cai, W. J. James and W. B. Yelon, Appl. Phys. Lett., 2006, 88, 41914.
68. J. F. Herbst and L. G. Hector, Jr, Phys. Rev. B, 2005, 72, 125120.
69. K. Miwa, N. Ohba, S.-I. Towata, Y. Nakamori and S.-I. Orimo, Phys. Rev. B, 2005, 71, 195109.
70. Y. Song and Z. X. Guo, Phys. Rev. B, 2006, 74, 195120.
71. A. D. Becke and K. E. Edgecombe, J. Chem. Phys., 1990, 92, 5397.
72. L. G. Hector and J. F. Herbst, J. Phys.: Condens. Matter, 2008, 20, 064229.
73. J. B. Yang, X. D. Zhou, Q. Cai, W. J. James and W. B. Yelon, Appl. Phys. Lett., 2006, 88, 41914.
74. C. Zhang, M. Dyer and A. Alavi, J. Phys. Chem. B, 2005, 109, 22089.
75. T. Noritake, H. Nozaki, M. Aoki, S. Towata, G. Kitahara, Y. Nakamori and S.-I. Orimo, J. Alloys Compd., 2005, 393, 264.
76. K. Ohoyama, Y. Nakamori, S.-I. Orimo and K. Yamada, J. Phys. Soc. Jpn., 2005, 74, 483.
77. M. P. Balogh, C. Y. Jones, J. F. Herbst, L. G. Hector, Jr and M. Kundrat, J. Alloys Compd., 2006, 420, 326.
78. B. Magyari-Köpe, V. Ozoliņš and C. Wolverton, Phys. Rev. B, 2006, 73, 220101.
79. T. Mueller and G. Ceder, Phys. Rev. B, 2006, 74, 134104.
80. H. M. Jin and P. Wu, Appl. Phys. Lett., 2005, 87, 181917.
81. S.-I. Orimo, Y. Nakamori, G. Kitahara, K. Miwa, N. Ohba, T. Noritake and S. Towata, Appl. Phys. A: Mater. Sci. Process., 2004, 79, 1765.
82. C. J. Zhang and A. Alavi, J. Phys. Chem. B, 2006, 110, 7139.
83. J. Rijssenbeek, Y. Gao, J. Hanson, Q. Huang, C. Jones and B. Toby, J. Alloys Compd., 2008, 454, 233.
84. J. K. Kang, S. Y. Kim, Y. S. Han, R. P. Muller and W. A. Goddard III, Appl. Phys. Lett., 2005, 87, 111904.
85. J. P. Soulié, G. Renaudin, R. Černy and K. Yvon, J. Alloys Compd., 2002, 346, 200.
86. A. Züttel, S. Rentsch, P. Fischer, P. Wenger, P. Sudan, P. Mauron and C. Emmenegger, J. Alloys Compd., 2003, 356–357, 515.
87. T. J. Frankcombe and G.-J. Kroes, Phys. Rev. B, 2006, 73, 174302.
88. K. Miwa, N. Ohba, S.-I. Towata, Y. Nakamori and S.-I. Orimo, Phys. Rev. B, 2004, 69, 245120.
89. Z. Łodziana and T. Vegge, Phys. Rev. Lett., 2004, 93, 145501.
90. T. J. Frankcombe, G.-J. Kroes and A. Züttel, Chem. Phys. Lett., 2005, 405, 73.
91. Y. Filinchuk, D. Chernyshov and R. Cerny, J. Phys. Chem. C, 2008, 112, 10579.
92. N. A. Zarkevich and D. D. Johnson, Phys. Rev. Lett., 2008, 100, 040602.
93. Y. Dmitriev, Y. Filinchuk, D. Chernyshov, A. V. Talzyin, A. Dzwilewski, O. Andersson, B. Sundqvist and A. Kurnosov, Phys. Rev. B, 2008, 77, 174112.
94. J. C. Galvez-Ruiz and M. Sanchez, J. Mol. Struct., 2007, 818, 23.
95. Q. Ge, J. Phys. Chem. A, 2004, 108, 8682.
96. A. J. Du, S. C. Smith and G. Q. Lu, Phys. Rev. B, 2006, 74, 193405.
97. M. B. Smith and G. E. Bass, Jr, Phys. Rev. B, 2006, 73, 174302.
98. A. J. Du, S. C. Smith, X. D. Yao and G. Q. Lu, J. Phys. Chem. C, 2007, 111, 12124.
99. S. Orimo, Y. Nakamori and A. Züttel, Mater. Sci. Eng., B, 2004, 108, 51.
100. K. Miwa, N. Ohba, S. Towata, Y. Nakamori and S.-I. Orimo, J. Alloys Compd., 2005, 404–406, 140.
101. L. Yin, P. Wang, Z. Fang and H. Cheng, Chem. Phys. Lett., 2008, 450, 318.
102. J. J. Reilly and R. H. Wiswall, Inorg. Chem., 1967, 6, 2220.
103. J. J. Vajo, S. L. Skeith and F. Meters, J. Phys. Chem. B, 2005, 109, 3719.
104. S. V. Alapati, J. K. Johnson and D. S. Sholl, J. Phys. Chem. B, 2006, 110, 8769.
105. S. V. Alapati, J. K. Johnson and D. S. Sholl, J. Phys. Chem. C, 2007, 111, 1584.
106. T. J. Frankcombe, J. Alloys Compd., 2007, 446–447, 455.
107. S. V. Alapati, J. K. Johnson and D. S. Sholl, Phys. Chem. Chem. Phys., 2007, 9, 1438.
108. S. V. Alapati, J. K. Johnson and D. S. Sholl, J. Phys. Chem. C, 2008, 112, 5258.
109. D. J. Siegel, C. Wolverton and V. Ozoliņš, Phys. Rev. B, 2007, 76, 134102.
110. V. Ozoliņš, B. Sadigh and M. Asta, J. Phys.: Condens. Matter, 2005, 17, 2197.
111. N. Hanada, W. Lohstroh and M. Fichtner, J. Phys. Chem. C, 2008, 112, 131.
112. R. Strobel, J. Garche, P. T. Moseley, L. Jorissen and G. Wolf, J. Power Sources, 2006, 159, 781.
113. A. C. Dillon, K. M. Jones, T. A. Bekkedahl, C. H. Hiang, D. S. Bethune and M. J. Heben, Nature, 1997, 386, 377.
114. A. Chambers, C. Park, R. T. K. Baker and N. M. Rodriguez, J. Phys. Chem. B, 1998, 102, 4253.
115. E. Yoo, T. Habe and J. Nakamura, Sci. Technol. Adv. Mater., 2005, 6, 615.
116. F. Darkrim and D. Levesque, J. Chem. Phys., 1998, 109, 4981.
117. R. F. Cracknell, Mol. Phys., 2002, 100, 2079.
118. P. Piseri, E. Barborini, M. Marino, P. Milani, C. Lenardi, L. Zoppi and L. Colombo, J. Phys. Chem. B, 2004, 108, 5157.
119. J. S. Arellano, L. M. Molina, A. Rubio and J. A. Alonso, J. Chem. Phys., 2000, 112, 8114.
120. Y. Okamoto and Y. Miyamoto, J. Phys. Chem. B, 2001, 105, 3470.
121. F. Tran, J. Weber, T. A. Wesolowski, F. Cheikh, Y. Ellinger and F. Pauzat, J. Phys. Chem. B, 2002, 106, 8689.
122. T. Heine, L. Zhechkov and G. Seifert, Phys. Chem. Chem. Phys., 2004, 6, 980.
123. C. Gu, G. H. Gao and H. Gao, Phys. Chem. Chem. Phys., 2002, 4, 4700.
124. D. Henwood and J. D. Carey, Phys. Rev. B, 2007, 75, 245413.
125. Y. Lei, S. A. Shevlin, W. Zhu and Z. X. Guo, Phys. Rev. B, 2008, 77, 134114.
126. S.-H. Jhi and Y.-K. Kwon, Phys. Rev. B, 2004, 69, 245407.
127. J. S. Arellano, L. M. Molina, A. Rubio, M. J. Lopez and J. A. Alonso, J. Chem. Phys., 2002, 117, 2281.
128. G. Vidali, G. Ihm, H. Y. Kim and M. W. Cole, Surf. Sci. Rep., 1991, 12, 133.
129. S. Patchkovskii, J. S. Tse, S. N. Yurchenko, L. Zhechkov, T. Heine and G. Seifert, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 10439.
130. S. Park, D. Srivastava and K. Cho, Nano Lett., 2003, 3, 1273.
131. J. Li, T. Furuta, H. Goto, T. Ohashi, Y. Fujiwara and S. Yip, J. Chem. Phys., 2003, 119, 2376.
132. E. C. Lee, Y. S. Kim, Y. G. Jin and K. J. Chang, Phys. Rev. B, 2002, 66, 073415.
133. I. Cabria, M. J. Lopez and J. A. Alonso, Comput. Mater. Sci., 2006, 35, 238.
134. X. B. Yang and J. Ni, Phys. Rev. B, 2006, 74, 195437.
135. G. E. Froudakis, Nano Lett., 2001, 1, 179.
136. A. Nikitin, H. Ogasawara, D. Mann, R. Denecke, Z. Zhang, H. Dai, K. Cho and A. Nilsson, Phys. Rev. Lett., 2005, 95, 225507.
137. S. Park, D. Srivastava and K. Cho, Nano Lett., 2003, 3, 1273.
138. S. P. Chan, G. Chen, X. G. Gong and Z. F. Liu, Phys. Rev. Lett., 2001, 87, 205502.
139. Y. H. Kim, Y. F. Zhao, A. Williamson, M. J. Heben and S. B. Zhang, Phys. Rev. Lett., 2006, 96, 016102.
140. A. J. Du and S. C. Smith, Nanotechnology, 2005, 16, 2118.
141. T. A. Halgren, J. Am. Chem. Soc., 1992, 114, 7827.
142. R. C. Lochan and M. Head-Gordon, Phys. Chem. Chem. Phys., 2006, 8, 1357.
143. V. V. Simonyan, P. Diep and J. K. Johnson, J. Chem. Phys., 1999, 111, 9778.
144. G. L. Che, B. B. Lakshmi, C. R. Martin and E. R. Fisher, Langmuir, 1999, 15, 750.
145. L. Grigorian, G. U. Sumanasekera, A. L. Loper, S. Fang, J. L. Allen and P. C. Eklund, Phys. Rev. B, 1998, 58, R4195.
146. G. H. Gao, T. Cagin and W. A. Goddard, Phys. Rev. Lett., 1998, 80, 5556.
147. M. Yoon, S. Yang, C. Hicke, E. Wang, D. Geohegan and Z. Zhang, Phys. Rev. Lett., 2008, 100, 206806.
148. M. M. Shaijumon, N. Bejoy and S. Ramaprabhu, Appl. Surf. Sci., 2005, 242, 192.
149. P. Chen, X. Wu, J. Lin and K. L. Tan, Science, 1999, 285, 91.
150. W. Q. Deng, X. Xu and W. A. Goddard, Phys. Rev. Lett., 2004, 92, 166103.
151. M. Khantha, N. A. Cordero, L. M. Molina, J. A. Alonso and L. A. Girifalco, Phys. Rev. B, 2004, 70, 125422.
152. S. M. Lee and Y. H. Lee, Appl. Phys. Lett., 2000, 76, 2877.
153. P. Dubot and P. Cenedese, Phys. Rev. B, 2001, 63, 241402.
154. G. E. Froudakis, Nano Lett., 2001, 1, 531.
155. Y. F. Zhao, Y. H. Kim, A. C. Dillon, M. J. Heben and S. B. Zhang, Phys. Rev. Lett., 2005, 94, 155504.
156. T. Yildirim and S. Ciraci, Phys. Rev. Lett., 2005, 94, 175501.
157. D. M. P. Mingos, J. Organomet. Chem., 2001, 635, 1.
158. G. J. Kubas, J. Organomet. Chem., 2001, 635, 37.
159. J. W. Lee, H. S. Kim, J. Y. Lee and J. K. Kang, Appl. Phys. Lett., 2006, 88, 143126.
160. B. Kiran, A. K. Kandalam and P. Jena, J. Chem. Phys., 2006, 124, 224703.
161. E. Durgun, S. Ciraci, W. Zhou and T. Yildirim, Phys. Rev. Lett., 2006, 97, 226102.
162. W. Zhou, T. Yildirim, E. Durgun and S. Ciraci, Phys. Rev. B, 2007, 76, 085434.
163. A. B. Phillips and B. S. Shivaram, Phys. Rev. Lett., 2008, 100, 105505.
164. S. Li and P. Jena, Phys. Rev. B, 2008, 77, 193101.
165. S. Li and P. Jena, Phys. Rev. Lett., 2006, 97, 209601.
166. Q. Sun, Q. Wang, P. Jena and Y. Kawazoe, J. Am. Chem. Soc., 2005, 127, 14582.
167. G. Kim, S.-H. Jhi and N. Park, Appl. Phys. Lett., 2008, 92, 13106.
168. S. A. Shevlin and Z. X. Guo, J. Phys. Chem. C, 2008, 112, 17456.
169. Y. F. Zhao, A. C. Dillon, Y.-H. Kim, M. J. Heben and S. B. Zhang, Chem. Phys. Lett., 2006, 425, 273.
170. N. Akman, E. Durgun, T. Yildirim and S. Ciraci, J. Phys.: Condens. Matter, 2006, 18, 9509.

---

Footnote

† Part of the renewable energy theme issue.

---