# Breaking the H2 chemical bond in a crystalline environment

Miriam Marqués  Centre for Science at Extreme Conditions and School of Physics
and Astronomy, University of Edinburgh, Edinburgh, U.K.  Miriam Peña Alvarez
Centre for Science at Extreme Conditions and School of Physics and Astronomy,
University of Edinburgh, Edinburgh, U.K.  Miguel Martinez-Canales  Centre for
Science at Extreme Conditions and School of Physics and Astronomy, University
of Edinburgh, Edinburgh, U.K.  Graeme J. Ackland  Centre for Science at
Extreme Conditions and School of Physics and Astronomy, University of
Edinburgh, Edinburgh, U.K.

###### Abstract

Through density functional theory and molecular dynamics calculations, we have
analysed various metal polyhydrides to understand whether hydrogen is present
in its molecular or atomic form - tetrahydrides of Ba,Sr,Ra, Cs and La; Ba8H46
and BaH12. We show that, in experimentally reported binary barium hydrides
(BaHx), molecular H2 and atomic H- can coexist with the metallic cations. In
this thorough study of differences between BaH4, higher barium hydrides, and
other binary tetrahydrides we find the number of atomic hydrogens is equal to
the formal charge of the cations. The remaining hydrogen forms molecules in
proportions yielding, e.g. BaH2(H${}_{2})_{x}$, at pressures as high as 200
GPa. At room temperature these are highly dynamic structures with the
hydrogens switching between H- and H2 while retaining the 2:x ratio.

We find some qualitative differences between our static DFT calculations and
previously reported structural and spectroscopic experimental results. Two
factors allow us to resolve such discrepancies: Firstly, in static relaxation
H2 must be regarded as a non-spherical object, which breaks symmetry in a way
invisible to X-rays; Secondly the required number of molecules $x$ may be
incompatible with the experimental space group (e.g. $BaH_{2}(H_{2})_{5}$). In
molecular dynamics, bond-breaking transitions between various structural
symmetry configurations happen on a picosecond timescale via an H${}_{3}^{-}$
intermediate. Rebonding is slow enough to allow a spectroscopic signal but
frequent enough to average out over the lengthscale involved in diffraction.

###### keywords:

High pressure, barium, hydrides, molecular hydrogen

Recent reports of superconductivity approaching room temperature 1, 2, 3, 4, 5
have galvanised attempts to synthesise hydrogen-rich materials. Much effort is
being invested to understand their crystal structures and properties6, 7, 8,
9, 10, 11, 12, 13. However, high-temperature hydride superconductivity is
associated not with a particular crystal structure, but with the stabilisation
of dense atomic hydrogen over molecular dimers 14, 15, 16. In this paper we
address the fundamental issue of whether hydrogen in metal hydrides is atomic
or molecular. 17, 18, 19

Molecular crystals provide a unique challenge for structural studies,
especially under high pressure. From the experimental side, a major challenge
is the very weak X-ray scattering by hydrogen atoms. Spectroscopy can only
give indirect information about the atomic structure, e.g. Raman signal
related to H2 molecules whose shift can be related to their intra-molecular
bond length.

Density Functional Theory (DFT) is the dominant theoretical method applied for
structural determination in these materials. An implicit assumption here is
that the lowest enthalpy static arrangement of the ions will be good
candidates for the the lowest Gibbs free energy state. In static relaxations,
the point symmetry of molecules may constrain the molecular orientation in a
given crystal symmetry. However, in practice, due to quantum or classical
rotation, H2 molecules may behave as spherical rather than linear objects,
which in turn allows high symmetry crystals20.

Molecular hydrogen is understood as two hydrogen atoms joined by a covalent
bond. In compression of pure hydrogen, the H–H bond lengthens from 0.718 at 6
GPa in phase I to 0.84 Å when entering phase IV at 225 GPa21, 22, 23. Although
there is no unique definition of a covalent bond in DFT, many measures can be
indicative. For example, a well-defined interatomic distance which persists
long enough for the bond to oscillate; a high value of electron localization
function (ELF) between atoms; a build up of charge in the region between the
atoms; a narrow Kohn-Sham band of energies which when projected onto localised
basis functions is localised on the bond. We will utilise all these methods to
build a coherent picture of the presence or absence of covalent bond in a
binary hydride solid.

Stable alkaline-earth metals (marked as M hereafter) are known to form binary
hydrides (Fig.1). The lowest hydrogen content is experimentally seen in the
dihydrides, MH2. These satisfy the octet rule, and no molecular hydrogen is
expected24, 25, 26, 27, 28. Tetrahydrides have been experimentally reported
for Ca, Sr and Ba, containing their binary hydride or salt (divalent metal
atom, M2+, anionic hydrogens, H-), and molecular hydrogen which was
interpreted as slightly charged, H${}^{{\delta}-}_{2}$29. Ca30, 31, 18 and
Sr32, 18 tetrahydrides adopt the $I4/mmm$ structure, both with $c/a$ ratios
around 1.6, while BaH4 forms an $I4/mmm$ structure with $c/a$ of 2.2, as well
as an $fcc$ structure 18. Raman spectroscopy supports H–H bonds in these
tetrahydrides30, 29, 18. Moreover, synthesis of $Pm\bar{3}n$ Ba8H46 (BaH5.75)
at 50 GPa33, and $Cmc2_{1}$ BaH12 at 90 GPa19 has also been reported. Neither
compound has an identifiable Raman contribution from the stretching H2 mode.
Therefore, barium hydrides represent an ideal system to study how the change
in stoichiometry alters the H2 molecule within the solid 17, 18, 19. Thus, a
systematic investigation of the mechanisms involved in the hydrogen intra-
molecular bond weakening and the influence of the metallic atomic number is
here undertaken as the route to understand the concepts of bond-weakening and
chemical pre-compression.

In this work we conduct a thorough theoretical analysis of BaH4 and compare
with other binary metal tetrahydrides with $I4/mmm$ symmetry, along with
Ba8H46 and BaH12. We have run an extensive series of DFT calculations.
Electronic structure is explored using the ELF and Bader (Quantum Theory of
Atoms in Molecules, QTAIM34) topologies. Details of the calculations, as well
as further data, may be found in the supplemental material (SM). In this paper
we will start by discussing the DFT approach on molecular hydrogen. We then
discuss quantitative results from static relaxations for all materials,
followed by qualitative results from molecular dynamics.

Breaking the H2 bond may occur due to mechanical pressure or electronic
effects. To illustrate the electronic effects in isolation, we calculated the
hydrogen molecule in a homogeneous electron gas 35, (SM4.3). For low electron
density the band structure comprises a clearly-defined H2 bonding state and a
free electron-like density of states which lies much higher in energy. At
electronic densities above 0.06 $e/\textrm{\AA}^{3}$ ($r_{s}\approx
3.0~{}a_{0}$), the bond breaks spontaneously. This is evidenced by a
discontinuous jump in the H-H separation, the vanishing of the flat Kohn-Sham
band defining the molecule and the appearance of atomic-hydrogen (H-) states
at the bottom of the free-electron band. This suggests an upper limit on the
electron density where molecular hydrogen can be found. Interestinglym the
maximum bond length before the transition is less than 0.9 Å.

The easiest approach to identify a bond is to define some interatomic cutoff
distance $R_{\textrm{bond}}$, and declare that any pair of hydrogen atoms
closer than this are _“bonded”_. Molecular dynamics requires $R_{bond}$ to be
large enough that an oscillating molecule is not marked as “breaking” whenever
the bond-length briefly extends beyond $R_{\textrm{bond}}$ at its largest
extent. Consequently $R_{\textrm{bond}}$ must be set larger than one might
expect a H-H bond to be. Based on the observation of a minimum in radial
distribution functions, $R_{\textrm{bond}}=1~{}$Å is a pragmatic choice. We do
not draw conclusions which are highly sensitive to this choice. We further
refine our definition by counting the number of atoms with precisely one
neighbour within 1 Å , and dividing by 2 to get the number of molecules. In
practice we find that bond-breaking usually occurs via the approach of a third
atom and formation of an intermediate H${}_{3}^{-}$ complex, a process which
our definition counts as one molecule throughout. Nevertheless, we record
every instance where an H-H distance passes through the 1 Å threshold — this
is used for an estimate of bond lifetimes. Such a coarse lifetime estimate is
not quantitatively meaningful - we use it only for comparisons between
structures and not, e.g. for calculating Raman linewidths 36, 37, 38. A
detailed discussion of the technical definition of a bond is given in the SM.

In static relaxation, including from MD snapshots, other measures from
electronic structure can be applied to detect molecules. In $k$-space, a
discrete set of electron bands may represent the bond, and this can be tested
using a projection of the Kohn-Sham bands onto maximally-localised Wannier
functions. On the other hand, in real space, and within Bader topology, H2
molecules are identified by the presence of a covalent bond between the
hydrogen atoms. It is characterized by the high value of the electron density
and the negative sign of the laplacian at the bond critical point (first-order
saddle point of the electron density connecting the two maxima located on the
atoms) (see SM). The charge of the atoms forming a molecule can be
unambiguosly calculated by integration of the electron density within the
topological atoms forming the molecule. They are defined by the union of the
electron density maxima (located on the atoms) with their attraction basins
and delimited by perfectly defined surfaces obeying the zero-flux condition
for the electron density. ELF is a relative measure of the electron
localization with respect to the homogenous electron gas (HEG). In general,
the ELF value approaches 1 in regions of the space where electron pairing
occurs (e.g., atomic shells, bonds and lone pairs). In analogy with QTAIM, a
partition of the space based on the ELF can be performed. It consists of non-
overlapping basins with well-defined chemical interpretation (cores, bonds,
lone pairs). Moreover, the basin populations come from integration of the
electron density within these regions. The signature of a H2 molecule within
the ELF topology is the existence of a high ELF isosurface enclosing both
atoms. Typically, a H2 molecule is considered to exist if the minimum ELF
value between the two hydrogen atoms is above 0.85. We find these methods to
be consistent with one another, and consistent with the much faster bondlength
criterion.

We performed extensive DFT geometry optimisation calculations, fully relaxing
the atomic positions. Important structures are shown in Figure 1, and we will
focus our detailed discussion on 50 GPa. Other structures, as well as details
at other pressures, are described in the SM. We emphasise none of the lowest-
enthalpy structures have the symmetry reported from X-ray experiments.

System | c/a | nearest H–H, Å | volume | o,t ELF basin population | metal
---|---|---|---|---|---
BaH4 $Cmcm$ H2 | 1.716 | 0.806 | 68.00 |  | no
BaH4 $Cmcm$ H3 | 1.469 | 0.944 (H${}_{3}^{-}$) | 68.12 |  | no
BaH4 $C2/c$ | 1.514 | 0.829 | 70.87 |  | no
BaH4 $C2/c$ I | 1.467 | 0.942 (H${}_{3}^{-}$) | 68.17 |  | no
BaH4 $C2/c$ II | 2.000 | 0.822 | 68.756 |  | no
BaH4 $R\bar{3}m$ | - | 0.831 | 72.44 |  | no
BaH4 (low) | 1.763 | 0.803 | 72.28 | 1.047 1.731 | no
BaH4 (high) | 2.323 | 1.495 | 65.96 | 1.367 1.337 | yes
SrH4 (low) | 1.712 | 0.798 | 61.42 | 1.050 1.772 | no
SrH4 (high) | 2.083 | 1.629 | 57.14 | 1.362 1.416 | yes
RaH4 (low) | 1.721 | 0.796 | 78.58 | 1.052 1.773 | no
RaH4 (high) | 2.338 | 1.655 | 72.22 | 1.390 1.362 | yes
CsH4 (low) | 1.642 | 0.768 | 69.38 | 1.029 1.451 | yes
CsH4 (high) | 2.821 | 1.489 | 67.98 | 1.258 1.098 | yes
LaH4 | 1.863 | 1.502 | 64.56 | 1.450 1.604 | yes
Table 1: Calculated PBE properties of $I4/mmm$-XH4 structures at 50GPa, and
five other structures18. c/a and volumes refer to the conventional cells.
“high” and “low” refer to the initial c/a ratio before relaxation. Atomic ELF
population centered on the octahedral site (forming two independent atoms or a
molecule) and located on the tetrahedral site respectively. LaH4 exhibits only
one minimum

$I4/mmm$ is the symmetry reported from X-ray diffraction measurements for the
cation positions in a number of tetrahydrides18, 39. For Sr, Ba, Ra and Cs
$I4/mmm$ tetrahydrides, we find two energy minima in $I4/mmm$ as a function of
c/a ratio (Table 1). Despite having the same symmetry, there is a clear
distinction between the two structures which have very different c/a ratios
and H–H separations. The high c/a version ($c/a\sim$ 2–2.8) is metallic and
has atomic H-. The low c/a version ($c/a\sim$1.7) has molecular H2 and is
semiconducting in alkaline earths, but metallic in Cs and La. Typically, the
atomic version has higher density and becomes stable as higher pressure.

Figure 2 and Tables 1 and S1 show the two $I4/mmm$ structures, ELF populations
and Bader charges. All hydrogens in the high c/a phases have broadly similar
charges, whereas, in the low c/a phases, the molecular hydrogens are
essentially neutral (the integration of the electron density within the H2
units according to both, Bader and ELF topologies gives a value close to 2
electrons) whereas atomic hydrogens have a sizeable negative charge. The
evolution of the ELF charges with pressure is shown in Fig.2b and S4. In low
c/a structures, in all compounds, a narrow band of states around 6 eV below
the Fermi energy is clearly associated with the H2 molecules (Fig.2).
Projecting that band onto maximally localised Wannier functions (Fig 2c) shows
this band to be primarily associated with a H–H $\sigma$ bond. Finally, Figure
2d shows the enthalpy crossover from the low-c/a to the high c/a with
pressure, (zero in LaH4 below 50 GPa where only a single minimum exists).

Figure 1: Crystal structures, with Ba and H atoms represented as blue and pink
spheres, respectively. ELF isosurfaces (ELF=0.85) shown in yellow are
associated with anionic hydrogens and molecules. From top left $P6/mmm$-BaH2,
$I4/mmm$-BaH4 (low c/a), $I4/mmm$-BaH4 (high c/a), $R\bar{3}m$-BaH4,
$Cmcm$-H3-BaH4, $Pc$-BaH5.75, $P2_{1}$-BaH12 at 50 GPa.

Figure 2: Results of DFT calculations for the isostructural forms of
$I4/mmm$-BaH4. Atomic $I4/mmm$-BaH4 has high c/a $\approx$ 2.40 , closest
hydrogen pairs are $\approx$1.5Å apart in the xz plane. Molecular
$I4/mmm$-BaH4 has low c/a $\approx$ 1.8 at 100GPa: bonds are along z-direction
$<0.85$Å (a) Density of electronic states at 10 GPa (b) Charges based on ELF
topology. (c) Maximally Localised Wannier Function whose projection
corresponds to the covalent band (d) Enthalpy differences between high c/a and
low c/a for various $I4/mmm$ compounds, n.b. LaH4 has only one minimum at low
pressure

In BaH4 the static $I4/mmm$ calculation favours the low c/a $\approx 1.76$,
but other BaH4 structures are calculated to be more stable than any $I4/mmm$18
(see SM). Some of these are similar to $I4/mmm$, e.g. $R\bar{3}m$ is simply a
rotation of the hydrogen molecules out of $I4/mmm$ symmetry leaving the Ba in
very similar positions. However, two structures with $Cmcm$ symmetry have
especially low enthalpy. These differ in enthalpy by only a few meV at 50 GPa,
but are chemically very different ($Cmcm$-H2 lattice angle $\gamma=141^{o}$
with H2 molecules and $Cmcm$-H3 $\gamma=114^{o}$ with H${}_{3}^{-}$ units, see
Fig.1). The $Cmcm$-H3 version is more stable with the PBE and BLYP
functionals, whereas LDA favours $Cmcm$-H2. Both $Cmcm$ structures contribute
two e-/f.u. to the “H2” band (Fig.3). In any case, these structures are very
different from either experimentally reported structure for BaH4.18

We now investigate how molecular dynamics (MD) simulations can illustrate the
dynamics of these systems. It would be prohibitively expensive to run accurate
simulations of phase transformations or melting with conventional MD, because
the Ba atoms are 137 times heavier than hydrogen and move too slowly.
Therefore, an efficient way to sample the phase space is to run ab initio
molecular dynamics in CASTEP with fictitious masses, NPT ensemble and ringing-
mode quenching40 (Details SM2.2, Table S3). The masses of H2 and Ba are set to
be equal, so that all atoms are able to move significantly. The MD runs are
analysed using radial distribution functions, lattice parameters, mean squared
displacements, and ab oculo using VESTA and VMD 41, 42. After each MD run,
snapshots from the trajectory were relaxed to 0 K. These provided low-symmetry
structure candidates, and allowed us to validate the bond assignment using ELF
or charge density methods. Full details from the molecular dynamics runs and
their analysis is given in the supplementary materials table S2. Here we
describe features common to all barium hydrides.

The mean squared displacement (MSD) shows that the hydrogens move much further
than Barium, but the MSD quickly flattens off, showing all the materials are
crystalline up to 600 K; not superionic nor liquid. The one exception was
BaH5.75 at 1000 K which went through a melting transformation after 1.2 ps.

Typical radial distribution functions (RDF) are shown in Fig.4. They all have
a well-defined peak between 0.7Å and 1.0Å which defines the molecule.
Integrating the first peak, or counting the number of bonded pairs, gives a
number of bonds in all cases consistent with BaH2(H${}_{2})_{x}$, where $x$ is
the number of bonds required to make up the stoichiometry. This observation
holds for BaH2, all BaH4 structures, BaH5.75 and BaH12. All peaks in the H-H
RDF are in much the same place, regardless of composition. The shortest Ba-H
distance represents atomic size, and the first H-H is as expected for a
molecule, but the second H-H peak being similar is surprising. Under pressure,
the molecular peak moves to larger bondlength and subsequent peaks move to
shorter separation as density increases.

We estimate relative molecular lifetimes between structures by counting how
many times an H-H distance is exactly 1Å(See Table S2). In all MD structures,
molecules have lifetimes of hundreds of femtoseconds, an order of magnitude
larger than their vibrational frequencies. The H$2$ molecular version of
$Cmcm$ has the longest lifetime, while $I4/mmm$ (high) has the shortest. The
bond breaking rate increases with temperature and pressure. Bond breaking is
normally mediated by a third hydrogen atom, which forms an H${}^{-}_{3}$
complex, which then breaks back into a molecule and an isolated atom, which
may originally have been in the molecule. These H${}^{-}_{3}$ units are
usually transient, but are also observed with two equal H-H distances in MD
and static relaxation of e.g. Cmcm-H3, implying that they represent a stable
structural motif. The number of singly-bonded atoms remains remarkably stable
at the $2x$ in all BaH${}_{2}($H${}_{2})_{x}$ compounds. It is robust to the
choice of bond cut-off, whereas the number of unbonded and double-bonded
hydrogens is very system- and time- dependent. Fourfold coordinated hydrogen
is never observed.

Turning to specifics, for $I4/mmm$, the two distinct $c/a$ ratios behave very
differently. In the low $c/a$ the molecules persist at 1 pfu, average
bondlength 0.81Å. They rotate continuously and the $c/a$ drops well below the
static value. To within sampling and finite size error, $c/a=\sqrt{2}$, which
means that in X-ray diffraction from the Ba atoms will be indistinguishable
from the reported fcc structure. At 1000K and 600K the unit cell remains
orthorhombic, but at 300K there is a fluctuation between $\gamma=90^{o}$ and
$100^{o}$. In the high c/a ratio version molecularization occurs between atoms
in the same xy plane, not those sharing the octahedral interstice. Our
algorithm detects 0.84 molecules/fu, but this varies during the run and is
sensitive to cutoff because these molecules are short-lived and long (0.95Å on
average).

The two $Cmcm$ structures are intermediate in density between atomic and
molecular $I4/mmm$. In molecular dynamics, their hydrogen atoms show no
rotation. This offers an immediate explanation for why the $I4/mmm$/$fcc$
structures are more stable at room temperature: there is a significant
rotational contribution to the entropy.

BaH5.75 has been synthesized experimentally and the Ba positions were
identified via X-ray diffraction to correspond to the Weaire-Phelan (A15 or
$\beta$-W) structure. This is a topologically close packed structure
($Pm\bar{3}n$) in which all interstitial sites are tetrahedra. Theory suggests
that the unusual 1:5.75 stoichiometry comes from placing one hydrogen in each
tetrahedron. If the hydrogens are placed in the tetrahedra, and relaxed
preserving $Pm\bar{3}n$ symmetry, the octet rule cannot be satisfied and
precisely one H2 molecule per Ba forms. To satisfy the octet rule (1.875
molecules per Ba) much larger supercells are needed, with at least 8 symmetry-
inequivalent Ba atoms.

The most stable structure we observe has $P_{c}$ symmetry with 30 H2 units
(identified from both RDF and ELF, Fig.1):
$\textrm{BaH}_{2}(\textrm{H}_{2})_{1.875}$. which has 3.75 electron bands per
Ba in ”H2” region of the bandstructure (Fig.3). Throughout the MD, bonds form
between atoms in adjacent tetrahedra: as a consequence the H2 molecules do not
rotate, but the structure is dynamics because bond breaking takes place via
H${}_{3}^{-}$ units.

We performed similar MD calculations on pseudocubic BaH1219, at pressures up
to 150 GPa. We find that any structure consistent with cubic symmetry is, in
DFT, highly unstable compared to H2 bond formation. Metastable structures with
broken symmetry involve at least four formula units (e.g. Ba4H48, P$2_{1}$
symmetry, or several variants with P1). The common theme of these structures
is molecularization, confirmed by ELF analysis (Fig.1): the molecular
distances range from 0.77 to 0.9 Å, with a clear gap to the next H-H distance
(1.02 Å) and the number of molecules follows the octet rule:
BaH2(H${}_{2})_{5}$.

By contrast, BaH2 is an unremarkable ionic insulator, with no close H–H
distances and no electrons in the energy region associated with H2
molecules(Figs.1,3). The ambient cotunnite structure transforms to P62/mmc at
high temperature and above 1.6GPa, and around 50 GPa it transforms to a
metallic P6/mmm 43, 44, 45, 46, 47

Figure 3 shows a remarkable similarity in the electronic density of states
between the BaH4 candidates and other barium hydrides of various compositions.
The extra electrons from hydrogen-rich stoichiometry are all located in the H2
band. Projections of the wavefunctions onto atomic orbitals enables us to
define distinct character to well-defined groups of bands, independent of
crystal structure. For BaH4 distinct bands can be associated with Ba $5p$,
molecular hydrogen and atomic hydrogen. All these solids are small-gap
semiconductors. This suggests that even at 50GPa, all these compounds can be
regarded as BaH2(H${}_{2})_{x}$. This does not hold for $I4/mmm$ with high
$c/a$, which has no molecules and no associated peak in Fig.3, and is highly
unstable in MD.

High ELF values are found at the bond points connecting all pairs of atoms
separated by $<1\textrm{\AA}$ (Fig.1). If the bands localized in energy space
at around $-6$–10 eV below the Fermi energy are projected onto maximally-
localized Wannier functions they are observed to also be localised in real
space.

Figure 3: Density of electron states from statically relaxed Barium hydrides
at 50GPa normalised to BaH2(H${}_{2})_{x}$ calculated using PBE. (a) BaH4
candidates, with four distinct bands associated with their largest projection
onto localised states Ba $5p$, H2, Ba $5d$+H, Ba $6s+5d$. (b) Other
stoichiometries. Cmcm refers to Cmcm-H2, Note that the BaH2 $P6_{3}/mmc$ (red
line) has no H2.

Figure 4:  Analysis of MD at 50 GPa/300 K. I4mmm/fcc indicates simulations in
$I4/mmm$ for which the c/a ratio is close to $\sqrt{2}$ (a) Fraction of bonds
relative to $x$ BaH2(H${}_{2})_{x}$ in as a function of time for various
structures. Each line is displaced by 0.5 for clarity (b) Radial distribution
functions, normalized to number of atoms pfu. The large peak at 2.4 Å
corresponds to BaH pairs.

We also studied the electron density and ELF by calculating “empty” $I4/mmm$
structures comprising just Ba without any hydrogen, and $I4/mmm$ BaH2 without
the hydrogen molecules. We found that for the low $c/a$ ratio where
interstitial hydrogen adopts a molecule form in the octahedral site, there is
a single ELF maximum centered on the octahedral site and elongated along the
direction of the hydrogens forming the molecule. At low pressures, there is no
difference in the ELF value for the associated ELF basin, but at bigger
pressures the ELF value at the ends of the basin is sligtly higher than that
at the octahedral site signalling the stabilisation of the atomic phase at
high pressure. Indeed, for the high $c/a$ phase stable at high pressures and
where hydrogen forms two atoms the original ELF basin is broken in two well
diferentiated ELF maxima with low ELF value in between. Interestingly, these
ELF maxima are not non-nuclear maxima of the electron density. Therefore, the
ability of high-pressure metals to host molecular hydrogen and their
corresponding position should be understood resorting to the ELF topology.

Calculations of substoichiometric $I4/mmm$ showed that the most stable defect
is to remove the hydrogen molecule altogether - replacing H2 with H in the
octahedral site is energetically unfavourable.

Localised H2 molecular units are well defined in five different metrics:
Bondlength, bandstructure, Wannier wavefunction projection, and in both Bader
and ELF topologies. The number of H2 molecules is determined by the octet rule
BaH2(H${}_{2})_{x}$: with x=0,1,1.875,5. This value for $x$ persists in
molecular dynamics of all these materials even after melting. In related
compounds, the octet rule holds and we observed CaH2(H${}_{2})_{x}$,
SrH2(H${}_{2})_{x}$ and RaH2(H${}_{2})_{x}$ but CsH(H${}_{2})_{x}$ and
LaH3(H${}_{2})_{x}$ are formed in molecular dynamics.

Calculated hydrogen positions for stable structures of BaH2(H${}_{2})_{x}$ are
inconsistent with the reported room temperature experimental symmetry. In the
case of BaH${}_{2}($H${}_{2})$ this can be explained by the MD calculations:
molecular rotation in the $I4/mmm$ structure leads to a $c/a$ ratio of
$\sqrt{2}$, for which the Ba positions are at the fcc sites as observed by
diffraction. Rotation means that there is no symmetry-breaking due to molecule
orientation. The rotation also means $I4/mmm$ has high entropy, which
stabilises it against the lower enthalpy $Cmcm$ phases, which are predicted to
be stable at low temperature. For the topologically close packed BaH5.75
compound, molecular rotation cannot occur because each hydrogen atom is
trapped in its own tetrahedra cage. The bonds are continuously being made and
broken - on a picosecond timescale across the face of the tetrahedra while the
Ba atoms remain in positions observed in experiment. Remarkably, at any given
time, the number of such bonds meets the BaH2(H${}_{2})_{1.875}$ octet
criterion. For BaH12 the decomposition into five molecules per unit cell is
inconsistent with the available Wyckoff sites in cubic symmetry. Once again,
molecular dynamics shows that continuous bond breaking and making occurs, such
that a hydrogen is sometimes bonded to each of its neighbours and on average
the high symmetry structure is recovered.

The instability of high symmetry structures to forming electron pairs as bonds
is likely to compete with the superconductivity instability towards Cooper
pair formation. Since our calculations are carried out using Kohn-Sham
wavefunctions, Cooper-pair correlations cannot be observed (see SM7).

Overall, this work reveals the persistence of the H2 bond in high pressure
compounds. Molecular rotation and dynamic bond breaking allow this to occur
without XRD-detectable long-range symmetry breaking. Furthermore, the
requirement from the octet rule to form BaH2(H${}_{2})_{x}$, and the
incompatibility of $x$ with simple Wyckoff sites, leads to frustration in the
bonding network and lowest-energy structures having large unit cells. This
study focused mainly on barium hydrides, but evidence from other hydrides
suggest these principles will be general.

For the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission. MM MPA and GJA would like to acknowledge the support of
the European Research Council (ERC) Grant ”Hecate” reference No. 695527. MPA
acknowledges the support of the UKRI Future Leaders Fellowship Mrc-
Mr/T043733/1.

Supplementary Material

## S1. How to define an H-H covalent bond?

There are many measures which can be used to define a covalent bond in an ab
initio molecular dynamics simulation. A good measure should be a parameter
which can distinguish between bonded and unbonded pairs.

###  0.1 S1.1 Bondlength

The easiest approach is to define some cutoff distance $R_{bond}$, and declare
that any pair of hydrogen atoms closer than this are ”bonded”. This can be
immediately evaluated from integrating the radial distribution function
$g(r)$.

For perfect crystals $g(r)$ is zero between shells of neighbours. In many
cases of md, including solid hydrogen, there of a large range where $g(r)=0$
and the hydrogen bonded pairs can be unambiguously identified, even though the
choice of $R_{bond}$ is ambiguous, because any cutoff with $g(R_{bond})=0$
will give the same results.

In the present work, $g(r)$ does not go to zero. This however there is often a
distinct minimum at around 1 Å. In the static structures we invariably finds
that $g(r=1~{}$Å $)=0$. $R_{bond}=1~{}$Å , therefore a pragmatic choice which
we adopt as part of the definition. However, one must be careful that any
measure of bonding is not highly sensitive to this choice.

The second ingredient of defining bonds from MD comes from investigating
”Bond-breaking” events. In practice we see two types of event.

  * •

A sensible definition requires $R_{bond}$ to be large enough that an
oscillating molecule is not marked as ”breaking” if the bondlength briefly
extends beyond $R_{bond}$ at its largest extent. This is why our $R_{bond}$
must be set larger than one might expect a H-H bond to be.

  * •

Real bond breaking events in the current system come from reactions between
atomic and molecular hydrogen. These involve an H- approaching an H2 molecule,
forming an intermediate state of an H${}_{3}^{+}$ unit, which then breaks so
that the ion is now part of a molecule, and one atom of the molecule is an
ion. We write this process as

H- \+ H${}_{2}\rightarrow$ H${}_{3}^{-}\rightarrow$ H2 \+ H-

where the ordering of the hydrogen atoms in the formula is assumed unchanged.
In isolation, this could be regarded as a double-well potential, but in the
condensed phase each molecule as many option for rebonding.

Our $R_{bond}=1$Å criterion is long enough to flag the H${}_{3}^{-}$ unit as
two bonds, and so we observe that the ”number of bonds” increases during these
reactions in a way that is very sensitive to the choice of $R_{bond}$. Thus we
refine the measure by counting the number of atoms with precisely one
neighbour within 1Å , and dividing by 2 to get the number of bonds. One can
readily verify that during the processes described above this measure give one
covalent bond throughout. It turns out that this measure is insensitive to
choice of $R_{bond}$.

###  0.2 S1.2 Topological analysis of the electron density and ELF

There have been different approaches to analyze the chemical bonding in
crystals. Probably, one of the most straightforward and mathematically
rigorous is the Quantum Theory of Atoms in Molecules developed by Bader.34 It
is based on the topology of the electron density. The topological analysis
first identifies the critical points of the scalar fields, i.e, points where
their gradient vanishes. These are then classified as maxima, first-order
saddle points, second-order saddle points, or minima. The maxima of the
electron density are usually located on the atoms, whereas the first-order
saddle points are associated with chemical bonds (denoted bond critical points
(b.c.p.)). The electron density and the Laplacian at these points are used to
characterize the bond strength and its character. For instance, in a molecule
such as H2 the b.c.p. is characterized by high electron density and negative
laplacian (electrons are locally concentrated). The crystalline space can be
also partitioned in topological atoms defined by the union of the electron
density maxima with their attraction basins and delimited by surfaces obeying
the zero-flux condition for the electron density. Their charge can be
calculated from the integration of the electron density within these regions.
These regions are non-overlapping and additive, recovering the total volume of
the crystal.

The electron localization function was introduced by Becke and Edgecombe in
1990 for the analysis in real space of the electron localization 48, and later
reinterpreted by Savin in terms of the Pauli kinetic energy density ($t_{p}$)
49. By definition, ELF is a relative measure of the electron localization with
respect to the homogeneous electron gas (HEG). In general, the ELF value
approaches 1 in regions of the space where electron pairing occurs (e.g.,
atomic shells, bonds and lone pairs). In analogy with QTAIM, a partition of
the space based on the ELF can be performed. It consists of non-overlapping
basins with well-defined chemical interpretation (cores, bonds, lone pairs).
Moreover, the basin charges come from integration of the electron density
within these regions. For example, a typical molecule will feature an ELF
maxima in the middle of the interatomic distance associated with the bond.
What is more, a region of the space (ELF maxima basin) can be attributed to
it, with an electron population obtained by integration of the electron
density in that region. Hydrogen constitutes an exception. A maximum does not
appear along the H-H bond since the hydrogen molecule only contains two
electrons. In this case, a high ELF isosurface encapsulates both atoms.
Typically, a H2 molecule is predicted to exist if the hydrogen atoms are
connected at a value of ELF higher than 0.85, i.e, if the minimum ELF value in
between the hydrogen atoms is above 0.85.

###  0.3 S1.3 Wavefunctions

A quantum wavefunction can be expanded in any complete set of orthogonal basis
functions. In a DFT calculation we use a plane wave basis set so there are no
molecular orbitals. Plane waves offer the big advantage that one can
systematically approach a complete set by increasing the cutoff, and they are
automatically orthogonal. A basis set of localised orbitals is typically non-
orthogonal and therefore simultaneously overcomplete and undercomplete.

To use the Kohn-Sham wavefunctions to identify molecules one must project them
from plane waves onto localised basis functions. This can be done using
maximally localised Wannier functions. Sup. Figs 5, and 6 show the projection
of the Kohn-Sham bands identified as bonds onto maximally localised Wannier
functions, the localisation in real space and the identification in reciprocal
space are clear.

Figure 5:  (right) Projection of Kohn-Sham band onto Maximally Localised
Wannier function for BaH4 $I4/mmm$ structure showing localisation on a single
H2 unit. (right) Bandstructure of BaH4 with colours showing projection onto H2
covalent bond orbitals.

Figure 6: (right) Projection of Kohn-Sham band onto Maximally Localised
Wannier function for the $Cmcm$-H3 structure. (a) Band structure: bands are
enumerated from lowest energy upward, including the Ba 6s (b) Projection of
MLWF onto Bloch states. (c) Projection onto Ba 5p MLWF states. (d) Projection
onto atomic H 1s MLWF states (e) Projection onto MLWF states for H${}_{3}^{-}$
(f) Isosurface of an H3 MLWF

Some 20% spillage of the low-lying Ba 6s states projected onto the Ba-p
orbitals give some indication of the ambiguity in using the projection: it is
good for qualitative measure. This identified localised states on the Ba and
the H+ as one expects for ionic materials (Sup. Fig.6c-d).

The H2 unit can be considered in terms of its two molecular orbitals, the
bonding $1\sigma_{g}$ corresponds to the band some 5eV below the Fermi energy,
while the antibonding $1\sigma_{u}$ state is above the Fermi energy. The H3
unit can also be considered in terms of molecular orbitals. Since it has three
atoms, H3 has precisely three molecular $s$-orbitals, most simply
written111Sometimes the three names ”bonding, nonbonding and antibonding” are
applied to these states. as 1$\sigma_{g}$, 1$\sigma_{u}$ and 2$\sigma_{g}$,
whose wavefunctions have zero, one and three nodes. H${}_{3}^{-}$ implies that
1$\sigma_{g}$ and 1$\sigma_{u}$ are occupied.

The Wannier approach was also applied to the Cmcm-H3 structure, which has
H${}_{3}^{-}$ units (Fig.6). The H3 1$\sigma_{g}$, 1$\sigma_{u}$ electronic
states are located in the same region of space. This enables the MLWF to
generate combination of bonding and nonbonding states localised on one end of
the H3 unit. Such a state is preferred by the MLWF procedure because it is
more localised than either 1$\sigma_{g}$ or 1$\sigma_{u}$ states which are
spread across three atoms. Projection of this MLWF state onto bandstructure
shows it occupies both the 5eV region where the H2 bands lie in other
compounds, and the DoS peak just below the Fermi Energy degenerate with he H-
states (Sup. Fig.6e-f). We propose that the in the molecule orbital picture
1$\sigma_{g}$ and 1$\sigma_{u}$ can be assigned to these energy regions, while
2$\sigma_{u}$ lies above the Fermi energy.

##  1 S2 Density functional and ab initio molecular dynamics calculations

###  1.1 S2.1 Codes and settings

For static relaxations we carried out density functional theory calculations
using the CASTEP code, with ultrasoft pseudopotentials and the PBE
functional50, 51. We used a plane-wave basis with a 700 eV cutoff and a
k-point spacing better than 0.04 Å-1. Key results were repeated using LDA and
BLYP functionals52, 53, 54, 55, which give similar qualitative results but
have some minor effect on transition pressures. Molecular dynamics were run
with a a cutoff energy of 463 eV and k-point grids up to 263 in CASTEP, with a
timestep of 0.5fs in the NPT ensemble, and were followed by geometry
optimisations of snapshots to zero temperature.

ELF and Bader analyses were carried out using the CRITIC2 program56, and
applied to all single unit cell static structures

Wannier analysis was carried out using the wannier90 code57

###  1.2 S2.2 Fictitious Mass Molecular Dynamics

The partition function in the NPT ensemble:

$Z=\sum_{i}\exp(-H_{i}/kT)=\sum_{i}\exp(-U_{i}/kT)\exp(-PV_{i}/kT)\exp(-KE/kT)$

Where $U_{i}$ is the potential energy of the microstate, $V_{i}$ is the volume
and KE is the Kinetic energy. In the classical limit, the kinetic energy is
proportional to temperature, so the final term factors out. We observe that
the first two terms are independent of the atomic masses, so we are free to
choose them as we wish.

In the case of barium hydride, the two orders of magnitude difference in mass
means that in conventional MD the bariums would scarcely move. However, the
fictitious mass approach enables us to correctly sample the partition
function, so phase transformations involving the Ba atoms can be observed.
Likewise a simulation run above the melting temperature will move into the
phase-space region of the melt. The partition function is correctly sampled,
so equilibrium thermodynamic properties are correct. However the dynamics of
the melting process, or the vibrational spectra, will not be correct.

The most challenging thing to determine in this partition function is the
definition of $i$ \- this indicates ”Microstates corresponding to the I4mmm
structure”. Of course, finite temperature microstates have no symmetry, and
NPT simulation in designed to allow for a phase change. In practice, we
monitor the box dimensions and the fraction of molecules to determine distinct
phases.

We derived this method independently, but it seems such an obvious idea that
we assume it is not new.

###  1.3 S2.3 Hydrogen in jellium

Breaking the H2 bond can be due to physical confinement in an interstial site:
”chemical pressure”, or by electronic effects ”antibonding”.

The embedded atom method is one of the most enduring and successful techniques
in materials modelling. It was introduced with the hypothesis that the binding
energy of hydrogen in a metal would depend on the electron density of the
other atoms, and there was an optimal value for this. Hence the electronic
contribution to bond breaking can be calculated in isolation by embedding the
H2 molecule in a homogeneous electron gas. This system has been studied in
detail35, and here we repeat that work with the methods employed in this
paper. For low electron density the band structure comprises a free electron-
like density of states which lies much higher in energy than a clearly-defined
H2 bond. Above 0.06 s $e$ Å-3: the bond breaks spontaneously, evidenced by a
discontinuous jump in the H-H separation, the vanishing of the flat Kohn-Sham
band defining the molecule and the appearance of an atomic-hydrogen state at
the bottom of the free-electron band (Sup. Sup. Fig. 7). This sets an upper
limit on where molecular hydrogen can be found.

Figure 7:  The maximum electron density into which an H2 molecule can be
embedded is 0.07 Å -3. a) Plot of bondlength vs DFT Calculation using PBE
hydrogen embedded in a jellium background at the shown electron density.
Calculations were done with one hydrogen molecule in a 343 Å 3 box with a 163
k–point grid. Similar calculations using BLYP give the bond breaking at the
same density b) Bandstructure calculations (density of states) for the system
described above with a total of 23 or 24 electrons in the system. The covalent
bond in the 23 electron system lies some 10 eV below the Fermi energy, or 4.5
eV below the bottom of the free electron band. After bond-breaking, the atomic
states lie at the bottom of the free electron band. Other oscillations in the
free-electron band structure come from the cubic symmetry of the supercell and
are not relevant.

##  2 S3 $I4/mmm$ structures

Figure 8:  Evolution of $c/a$ with pressure for the $I4/mmm$ structures. Solid
and open symbols correspond to the low and high $c/a$ structures. The $I4/mmm$
structures with the lowest enthalpies at each pressure are indicated by solid
lines.  Figure 9:  Evolution of $a$ lattice parameter with pressure for the
$I4/mmm$ structures. Solid and open symbols correspond to the low and high
$c/a$ structures. The $I4/mmm$ structures with the lowest enthalpies at each
pressure are indicated by solid lines.

Figure 10: a) Charges as a function of pressure in $I4/mmm$ compounds based on
the ELF topology for the metal and hydrogen atoms located on 4$e$ and 4$d$
sites

###  2.1 S3.1 Comparing I4/mmm materials

We considered $I4/mmm$ structures for SrH4, LaH4, CsH4 and RaH4 to determine
trends across the periodic table making no definitive claim about the
stability of the calculated compounds. C/a ratio, lattice parameters and ELF-
basin populations are shown in Sup. Figs. 8, 9 and 10.

BaH4 is discussed in detail in the main text. It is unusual in that the low
c/a structure calculated to be stable in DFT is different from the high c/a
ratio structure seen in the experiment. Nevertheless, both structures are
local energy minima. Many other materials exhibit the same theoretical
coexistence of two energy minima in the same I4/mmm structure.

SrH4: Strontium is above Barium on Group 2A, and calculations show the atomic
and molecular forms do adopt separate enthalpy minima, with the atomic form
only becoming stable above 185GPa. Experiments on SrH4 have not reached these
pressures. The molecular form is non-metallic, e.g. at 50GPa there is a
bandgap of about 2eV (Sup. Fig.11)

CsH4: Caesium lies beside Barium, and calculations show the atomic and
molecular forms do adopt separate enthalpy minima, but the atomic form is
unstable up to 200GPa. This is despite the fact that, being monovalent, the
material is calculated to be metallic with the Fermi level in the atomic
hydrogen band. Being so chemically unfavoured suggests that this structure
will not be stable, probably decomposing into CsH 58 and hydrogen or a
polyhydride.

RaH4: Radium lies below Barium in the periodic table, and is highly
radioactive so the nature of its high-pressure hydrides is unlikely ever to be
tested. Calculations suggest that at all pressures up to at least 200GPa the
low c/a molecular structure is stable. The band structure is similar to BaH4
with somewhat broader bands and a DFT/PBE band gap which closes at 50GPa,

LaH4: Lanthanum is next to Ba in the periodic table, and typically forms
trivalent ions. Previous work 39 studied the high-pressure form of LaH4. The
pressure dependence of the bond in the LaH4 shows a trend towards steady
lengthening with pressure, with no distinct high and low variants. At 50GPa
two short H-H separations are already 1.5Åapart in the density of states, and
the characteristic “molecular” peak has vanished (Sup. Fig.11) The LaH4 is
always metallic, even in the molecular case, with the Fermi level lying in the
La $sd$ conduction band. Given this, it is likely that LaH4 is unstable with
respect to higher polyhydrides.

The band structures reveal significant similarities between the different
materials , and differences between atomic and molecular forms (see Density of
states at 50GPa). For BaH4, we calculate five distinct groups of state: low-
lying Ba 5s (-22eV) and 5p orbitals (-10eV). Then at -4eV come the H2 (0.79Å )
molecular bonding orbitals. Closest to the Fermi energy are the atomic
hydrogen H- hybridized with Ba 5d: population analysis in atomic basis sets
suggests that the Ba 6s are unoccupied.

In the Sr compound, this allocation is even clearer. In Ra, the electronic
bands are broader but energy gaps in the density of state mean there is no
hybridization. At 50GPa the H2 bond in both Sr and Ra is 0.79(1)Å

Overall, the comparison of the adjacent materials indicates that Ba is unusual
in showing the atomic-molecular transformation at the lowest pressure. The
band structures are remarkably similar across all materials, with the position
of the Fermi energy determined by the cation valence. We have not done
extensive structure searching: our molecular dynamics simulations indicate
that the molecular hydrogens are orientationally disordered at room
temperature, while quantum nuclear effects are likely to be important at zero
temperature. Given the emerging picture that these materials can be thought of
as dihydrides with interstitial H2 it is unlikely that the metallic compounds
with monovalent and trivalent cations are stable.

System | Bader charge | ELF basin population
---|---|---
| on M (e) | on Ho (e) | on Ht (e) | on M (e) | on Ho (e) | on Ht (e)
$I4/mmm$-BaH4 (low) | +1.174 | -0.041 | -0.546 | 8.446 | 1.047 | 1.731
$I4/mmm$-BaH4 (high) | +1.036 | -0.253 | -0.265 | 8.591 | 1.367 | 1.337
$I4/mmm$-SrH4 (low) | +1.247 | -0.049 | -0.574 | 8.355 | 1.050 | 1.772
$I4/mmm$-SrH4 (high) | +1.143 | -0.250 | -0.321 | 8.444 | 1.362 | 1.416
$I4/mmm$-RaH4 (low) | +1.278 | -0.047 | -0.591 | 8.350 | 1.052 | 1.773
$I4/mmm$-RaH4 (high) | +1.136 | -0.277 | -0.291 | 8.497 | 1.390 | 1.362
$I4/mmm$-CsH4 (low) | +0.691 | -0.031 | -0.315 | 8.040 | 1.029 | 1.451
$I4/mmm$-CsH4 (high) | +0.538 | -0.194 | -0.075 | 8.290 | 1.258 | 1.098
$I4/mmm$-LaH4 | +1.536 | -0.322 | -0.446 | 8.893 | 1.604 | 1.450
Table 2: Bader charges and ELF basin populations on the metal atoms (M),
hydrogens on 4$e$ sites (forming hydrogen pairs) and hydrogens on 4$d$ sites
(tetrahedral sites) for the $I4/mmm$ structures at 50 GPa  Figure 11: Band
structures at 50GPa in SrH4, BaH4, LaH4, CsH4 and RaH4. The zero of energy has
been offset to place the peak of the H2 ”band” together, as indicated by the
arrow. The calculated peaks contain 2 (not shown), 6 (double peak, widely
separated e.g. in Cs), 2 and 4 (double peak) electrons per formula unit
respectively. In BaH4 these can be interpreted chemically as Ba 5s, Ba 5p,
molecular H${}_{2}\sigma$, and atomic H 1s. The Fermi energy lies in the gap
at +6 eV and the conduction band has 5d and 6s character.

###  2.2 S3.2 Can the atomic and molecular forms coexist in BaH4?

When two different forms of hydrogen can be found in a single structure, it
opens the question of whether the real material should contain a mixture of
each. This could naively be treated as a mean field two-state system, with the
probability of finding a given state being $e^{-\Delta G/kT}/[1+e^{-\Delta
G/kT}]$ with $\Delta G$ being the free energy difference. At room temperature
this would be close to 50%, and no discontinuous phase transformation could be
observed.

However, the strong coupling between the atomic/molecular transition and the
strain suggests that this view is overly naive. Stabilization of the atomic
form occurs only in concert with a large enlongation of the octahedral
interstice, which is incompatible with a molecular form in an adjacent site. A
more appropriate model is that of Bragg and Williams (BW) which explicitly
includes an interaction enthalpy between adjacent sites. BW considered the
interaction to be mainly a short-ranged bonding effect, however here it seems
that strain provides the major contribution.

To test this we ran a series of molecular dynamics type calculations on a
number of cells. The purpose of these was to sample the available phase space,
rather than study dynamical trajectories. The partition function is
independent of the atomic masses, we set the H2 and Ba masses to be equal. The
results were monitored by tracking all dihydrogen pairs, defined geometrically
as atoms less than 1Å apart. We record each occasion when these connections
change from one timestep to the next as an indication of how often the bonds
break. We also calculated the radial distribution functions, the mean squared
displacements, by visualisation of the trajectories using vmd, and by
relaxation of selected snapshots.

Long MD runs of single unit cells, just 10 atoms, followed by static
relaxation was used to track dynamic instability against $\Gamma$-point
phonons, it was found that starting from the low-c/a:

  * •

Molecular I4/mmm 50GPa with 300K rotated molecules and rebonded, relaxation
gives structures with H${}_{3}^{-}$ units.

  * •

Molecular I4/mmm 50GPa with 600K transformed to a Cmcm-H3-like structure by
breaking and reforming bonds rather than rotating molecules.

  * •

Molecular I4/mmm 200GPa 300K transforms to the atomic phase, with molecules
diffusing in the plane by a bond-making and breaking mechanism.

These small simulations with unphysical masses allow us to span the phase
space. We also ran larger cells with 80 atoms and correct masses (137:1) to
investigate the dynamics. We ran NPT simulations with a Parrinello-Rahman
barostat which allows the cell to fluctuate but only permits small relative
movements of the Ba atoms. At 50GPa MD shows that the hydrogen molecules
retain their identity and rotate rather than remain orientated along the $c$
axis. A consequence of that is further reduction in c/a so that the mean
lattice parameters are a=3.6 (3.38), b=3.8 (3.38) c=5.2 (6.07), with 0K values
in brackets. At 600K a similar scenario plays out, with molecules breaking and
reforming on a sub-picosecond timescale.

Constant features across all MD on BaH4 are the creation of close to one H2
pair per formula unit, and the non-metallic bandstructure of the relaxed
structures.

Figure 12: Structures from $I4/mmm$ MD simulations. (left: MD snapshot
50GPa/300K; centre: snapshot after relaxation with H${}_{3}^{-}$ units; right:
200GPa snapshot after relaxation with H2 units)

Figure 13:  At room temperature both high and low c/a structures of I4/mmm
BaH4 contain molecules Snapshots from PBE MD simulation of (a) high
$<a>=5.38,<b>=5.40,<c>=12.29$ and (b) low MD simulations at 200GPa and 300K
$<a>=5.807,<b>=5.50,<c>=11.2$. Analysis of the number of dihydrogen pairs at
less than 1Å shows a distinct difference. The low c/a version averages 0.99
molecules per formula unit across the simulation. The high c/a has ”only” 0.85
pfu, however this is very high considering that the relaxed structure which
has zero. Much of the difference arises from the more rapid breaking and
remaking of bonds in the high c/a ratio case. These bonds typically last for,
of order, 0.1ps, consistent with 300cm-1 lifetime broadening in Raman peaks.
The different nature of the two structures is evident from the pictures, the
high c/a looks more like alternating Ba and H layers, which the low c/a has
hydrogen primarily in interstitial locations.

Figure 14:  Crystal Structures for BaH4, I4/mmm low and high c/a (top);
Cmcm-H2 and Cmcm-H3 (bottom) Bonds are shown for Ba-H less than 2.4Å, H-H less
than 1Å.

###  2.3 S3.3 Stability of $I4/mmm$ BaH4

Molecular dynamics runs (Sup. Figs. 12 and 13) show that the hydrogens in
I4/mmm are dynamically disordered and only obey I4/mmm symmetry on average.
Moreover, geometry optimisation from molecular dynamics runs finds different
symmetries, including the Cmcm-H3 which appears to be to low temperature
ground state (Sup. Fig.14).

###  2.4 S3.4 Substoichiometric $I4/mmm$ materials

Figure 15:  Density of states of $I4/mmm$ BaH4 at 50GPa with 50% (red) and
100% (black) of H2 units converted to H. The covalent states around -4 eV are
entirely absent when all molecules are removed. The added H atoms contribute 2
states and 1 electron to the H- bands, leading to metallic structures.

In experiments, BaH4 is formed from BaH2 in a hydrogen-rich environment. We
therefore investigated substoichiometric versions of the $I4/mmm$ structure,
replacing the H2 unit with a single H molecule.

With all H2 units changed to H, the stoichiometry is BaH3 and the c/a ratio
relaxed to $\sqrt{2}$. the hydrogens remain in the high symmetry tetrahedral
and octahedral sites, such that the overall symmetry of this structure is
Fm-3m and the Ba are located on the sites of an fcc lattice. With 50% of H2
changed to H, we find the c/a at 1.6.

The band structure of these off-stoichiometric structures have a clear
signature of the disappearance of the H2 bands and a broadening of the H- band
leading to closure of the bandgap (Sup. Fig.15).

Nevertheless the fcc all-atomic BaH3 structure in DFT is highly unstable to
molecule formation - by about 0.7eV per formula unit, whether one compares
decomposition into cotunnite BaH2 and BaH4, or a BaH3 structure with one
molecule and one empty octahedral site.

In molecular dynamics simulation with a single ”missing” H we see that the
hydrogen can diffuse on a sub-picosecond timescale. The mechanism for this
diffusion is the transfer of a hydrogen atom to molecule forming an
H${}_{3}^{-}$ complex, with one of the original pair moving away to become a
free atom, leaving a new pairing. We observed this migration mechanism by
inspection in vmd and by monitoring the graph of connections $<$1Å, between
molecules in the MD.

##  3 S4 Empty lattices

We considered the ELF in the BaH4 structures with the H2 removed at a range of
pressures (Sup. Figs.16 and 17 ).

Figure 16: ELF value at the 2$b$ site (oct), 4$d$ site (tet) and the ELF
maxima on 4$e$ sites (maxoct) with the 2$b$ site as middle point for the empty
Ba lattices of low and high $I4/mmm$-BaH4.

Figure 17: a) ELF isosurface (ELF=0.32) for the pure $Pm\bar{3}n$ Ba
sublattice of BaH5.75 at 50 GPa (in yellow). Ba atoms represented as blue
spheres and ELF maxima as green, pink and grey spheres. b) ELF isosurfaces,
ELF=0.95 (in orange) and ELF=0.65 (in yellow) for $Pm\bar{3}n$ BaH5.75 at 50
GPa. Ba atoms represented as blue spheres and the 3 non-equivalent H atoms as
green, pink and grey spheres, respectively.

##  4 S5 Other BaHx compounds

$Cmcm$-H3 is an interesting structure where we identify H3 units (H (8$f$)-H
(4$c$)1-H (8$f$)). In this structure, Ba donates electrons to both, H3 (4$c$)2
and the H3 units, that hold similarly high negative charges. This structure
should be understood as BaH(H$3$). We also calculated the ELF for the high-
symmetry structures BaH5.75 and its empty-lattice equivalent, and for the
high-pressure superconductor BaH12 (Sup. Fig. 17, 18).

Atom | x | y | z | Bader charge | ELF charge
---|---|---|---|---|---
Ba (4$c$) | 0.000 | 0.117 | 0.250 | +1.164 | +1.485
H (4$c$)1 | 0.000 | 0.402 | 0.250 | +0.085 | +0.320
H (8$f$) | 0.000 | 0.610 | 0.565 | -0.354 | -0.540
H (4$c$)2 | 0.000 | 0.810 | 0.250 | -0.541 | -0.732
H3 |  |  |  | -0.623 | -0.754
Table 3: Bader charges and ELF basin populations for the $Cmcm$-H3 structure
at 50 GPa (a=3.462 Å. b=7.740 Å, c= 5.080 Å).

Figure 18: ELF value at the middle H-H point for H-H distances up to 1.5 Åat
50 GPa. a) BaH5.75, b) BaH12

##  5 S6 Details of Molecular Dynamics of specific materials

###  5.1 S6.1 MD of $I4/mmm$ BaH4

We ran a series of MD simulations on BaH4 with high and low c/a ratios, at
50GPa and 200GPa, and at 300K, 600K and 1000K.

Figure 19: Radial distribution functions from molecular dynamics simulations
in $I4/mmm$ BaH4 at various conditions. Each graph shows total g(r) normalised
to total atomic density.

At 50GPa/300K, (Sup. Fig.12) the I4/mmm simulation box transformed to a
distorted structure, which after relaxation has a cell ( a = 7.078 , $\alpha$
= 86.9, b = 7.593 , $\beta$ = 89.6, c = 10.624 , $\gamma$ = 101.7) this
contrasts with molecular dynamics with ”correct” Ba mass where the Ba motion
was too slow to allow for the transformation and rotation of the hydrogen
molecules was observed. A distinctive feature of this structure is the
$H_{3}^{-}$ units, as shown in figure 12. This type of “polyhydrogen” chain
has been associated with insufficient k-point sampling59, but here it is
stable in grids up to 12x12x6 in an 80-atom supercell. In the RDF (Sup.
Fig.19) these units manifest as a shoulder on the first H2 peak. In figure 12
we show that these H${}_{3}^{-}$ units in the relaxed structure are largely
parallel, but including two which have different orientations. The
H${}_{3}^{-}$ units were previously noted in a metastable $C2c$ structure, so
the relaxed MD structure appears to be this structure with misoriented
defects. The bond-breaking analysis detects many bond breaking events, but on
closer inspection these are almost all within the H3 triplets. This means of
bond-count is remarkably stable since both H${}_{3}^{-}$ units and H2
molecules contain 2 atoms with one neighbour, so if a third atom approaches a
molecule, or an H${}_{3}^{-}$ unit breaks, the measure is unchanged.

For same 50GPa/300K conditions starting with high c/a=2.323, the c/a ratio
remained roughly constant. Despite this, the radial distribution function is
remarkably similar to the low c/a version, dominated by H2 molecules, with an
H${}_{3}^{-}$ shoulder below 1 Åand the Ba-H peak at around 2.3 Å.

At 200GPa/300K the molecules, initiated along the z-direction, rotated into
the x-direction and rebonded, the simulation box remained close to
orthorhombic. The H2 molecules retained their integrity, although in the
relaxed final configuration each H2 bond (0.87 Å , Mulliken bond population
$\equiv 1$) is associated with a third hydrogen (1.02 Å , Mulliken bond
population $\equiv 0.5$). These objects seem to be intermediate between the
symmetric H3 units and independent H2 and H-. At these pressures the molecular
rotation is strongly suppressed, giving rise to a much more structured radial
distribution function. The high c/a version RDF is completely different the H2
peak is poorly defined and now at 0.95 Å. Bond lifetimes are half that of the
low c/a version

At 50GPa/600K the molecules started to rotate and c/a ratio dropped abruptly
from the static value of 1.76 to 1.32, with $b<a$ by about 10%. Taking the
mean of 2c/(a+b) across the run gives 1.41, the ratio at which the Ba atoms
lie on an fcc lattice.

Figure 20:  Equal-mass MD can be used to demonstrate melting : molten Ba
Details from the ”equal-mass” MD simulation in BaH5.75 which melts at 1000K
100GPa a) Snapshot of the molten configuration (doubled cell, periodic
boundary conditions suppressed). b) Mean squared displacements of the ”light”
Ba (red) and Hydrogen (black, also x,y,z components showing homogeneous
motion)

Figure 21: BaH5.75 as BaH${}_{2}($H${}_{2})_{x}$ (left) Radial distribution
function of the 100GPa/300K simulation, compared with those at 1000K in liquid
and solid form. (right) Pre- and post-melting radial distribution functions,
and the breakdown by species. Analysis of the bonding shows that the MD in
typically 94% fully bonded, compared to the low symmetry, relaxed structure
BaH2(H${}_{2})_{x}$ (227 compared to the ideal 240/368 hydrogens in bonds).
The lifetime of these bonds, as measured by counting changes in the bond list,
is given in Table 3.

Figure 22: The most stable found structure (Pc - Ba16H92) of BaH5.75 relaxed at 50GPa (left) Showing 30 H2 bonds pfu in range 0.813-0.936Å. (i.e. BaH${}_{2}($H${}_{2})_{x}$ \- next shortest is 1.049Å. (right) Showing Ba-H bonds less than 3Å, illustrating that each hydrogen has precisely 4 Ba neighbours (i.e. is confined to a tetrahedron).  System | Time (ps) | T (K) | P (GPa) | Mass | bonds | Events N,$ps^{-1},$ | Bond
---|---|---|---|---|---|---|---
$Cmcm$-H2 (Ba16H64) | 2.5 | 600 | 50 | 4 | 32,0 | 715 $\implies$ 4.5 | 0.80
$Cmcm$-H2 (Ba16H64) | 10 | 300 | 50 | 4 | 32,0 | 12 $\implies$ 0.01 | 0.79
$Cmcm$-H3 (Ba16H64) | 5 | 300 | 50 | 4 | 32,7 | 16317 $\implies$ 25.5 | 0.91
$I4/mmm$ (high, Ba16H64) | 5.6 | 300 | 200 | 2 | 27,2 | 21058 $\implies$ 69.2 | 0.95
$I4/mmm$ (high, Ba16H64) | 5.0 | 300 | 50 | 2 | 32,2 | 6080$\implies$19.0 | 0.81
$I4/mmm$ (low, Ba16H64) | 5.0 | 300 | 200 | 2 | 32,6 | 12154 $\implies$ 38.0 | 0.84
$I4/mmm$ (low, Ba16H64) | 4.4 | 1000 | 50 | 2 | 32,1 | 8901 $\implies$15.8 | 0.80
$I4/mmm$ (low, Ba16H64) | 3.9 | 600 | 50 | 2 | 32,6 | 4590$\implies$ 18.4 | 0.82
$I4/mmm$ (low, Ba16H64) | 7.0 | 300 | 50 | 2 | 32,0 | 8019$\implies$ 15.9 | 0.81
$I4/mmm$ (low, Ba16H63) | 4.8 | 300 | 50 | 2 | 31,1 | 5805 $\implies$19.5 | 0.81
WP(Ba64H368) | 6.5 | 300 | 100 | 137 | 227,10 | 75761$\implies$ 25.7 | 0.85
WP(Ba64H368) | 0.5 | 700 | 100 | 2 | 227,10 | 4462$\implies$ 19.7 | 0.80
WP(Ba64H368) | 0.3 | 700 | 0 | 2 | 238, 0 | 1538 $\implies$ 10.8 | 0.76
WP(Ba64H368) | 1.2 | 1000 | 100 | 2 | 175,10 | 36685$\implies$ 67.4 | 0.89
liquid (Ba64H368) | 4.5 | 1000 | 100 | 2 | 202,17 | 80194 $\implies$ 44.1 | 0.83
WP(Ba16H92) | 5.0 | 300 | 50 | 2 | 59,1 | 7514 $\implies$ 12.8 | 0.82
Ba8H96 | 2.1 | 300 | 50 | 2 | 80,1 | 1855 $\implies$ 5.5 | 0.77
Ba8H96 | 2.1 | 300 | 100 | 2 | 80,4 | 3526 $\implies$ 10.6 | 0.78
Ba8H96 | 5.0 | 300 | 150 | 2 | 80,8 | 10387 $\implies$ 13.0 | 0.78
Ba8H96 | 7.0 | 300 | 250 | 2 | 77,8 | 26254 $\implies$ 24.4 | 0.80
Table 4: Details of molecular dynamics simulations. WP is the Weaire-Phelan
structure, Mass refers to the cation, bonds is the average number of hydrogen
atoms within 1Å of precisely (one,two) other atoms. Broke is the number of
times and atoms move into or out of the 1Å range of another, and this is
normalise with time and number of bonds to give a number of such ”events” per
bond per ps. Relating this to the lifetime of the bond is complicated. e.g.
harmonic bond oscillating beyond 1Å will record four ”events” (two atoms,
making and breaking) without being a genuine bond breaking. The common proton
exchange process also records four ”events” (one atom makes a bond, one breaks
a bond, the central atom both makes and breaks). However in most of the
recorded H3 creation events, the ”reaction” is ”unsuccessful” and the
initially-atomic hydrogen leaves as an atomic hydrogen. Taken together, a
rough estimate is that bond breaking occurs with a frequency one order of
magnitude less than the ”Event” rate. Bondlength is the position of the first
peak in the radial distribution function.

###  5.2 S6.1 MD of other BaH4 candidates

MD simulation in Cmcm-H2 at 50GPa shows this structure remaining stable for at
least 5ps, with well-defined H2 and H- (Sup. Fig. 23) . Nevertheless, there
are hints of the asymmetric H${}_{3}^{-}$, with about 10% of H2 dimers being
associated with a third H within the first minimum of the radial distribution
function. Likewise, the Cmcm-H3, started in the stable state with
H${}_{3}^{-}$ units remains very close to the original structure, with 32
molecule throughout. Our bond-breaking analysis marks this structure as having
frequent bond breaking, but this is largely an artefact of fructuations in the
H${}_{3}^{-}$ units, which have mean bondlengths above 0.9Å and are therefore
susceptible to H-H distances greater than 1Å appearing in thermal
fluctuations.

Figure 23: Radial distribution functions from molecular dynamics simulations
starting in the $Cmcm$H2 structure. Each graph shows total g(r) normalised to
atomic density, and contributions from each bond type.

###  5.3 S6.2 MD of BaH5.75

The experimentally-observed BaH5.75 $Pm\bar{3}n$ structure is stabilized by
temperature, static relaxation calculations suggest a symmetry breaking
distortion is which the ”atoms” pair up to become molecular. At 50GPa, the
minimum energy structure has only Pc symmetry, with a clear separation between
the 15 H-H separations below 0.94Å and with bond population above 0.73e-, and
the 16th closest dihydrogen pair (1.05Å /0.5e-). These numbers demonstrate
that in conventional chemistry terms the Ba8H46 structure could be regarded as
a solid solution 8BaH2+15H2, rather than a Barium hydride. Further evidence
for this comes from the electronic density of states which strongly resembles
that in molecular I4mmm BaH4, with distinctive bands associated with low-lying
Ba 5s (-22eV) and 5p orbitals (-10eV), the H2 molecular bonding orbitals.
Closest to the Fermi energy are the hybridised bands of the Ba 5d and atomic
hydrogen: population analysis in atomic basis set suggests that the Ba 6s are
unoccupied.

The Weaire-Phelan topologically close-packed BaH5.75 phase was run in an MD
simulation with 64 formula units (Ba64H368 at 100GPa at 300K with full-mass Ba
and 100GPa at 700K and 1000K with equal-mass MD. The 700K run was continued at
0 GPa. Molecular dynamics of the Ba64H368 supercell (100GPa/300K) reveals a
different picture from the static case (Sup. Figs.20 and 21). The
molecularization to approximately Ba64H128(H${}_{2})_{120}$ occurs, but the
covalent bonds are continuously breaking and reforming. The topological close
packing means that all interstitial regions are tetrahedral, and there is
precisely one hydrogen atom in each tetrahedron; However, the hydrogens are
able to pair up with covalent bonds forming through the triangular faces of
the tetrahedra, and analysis shows that on average 114 bonds exist at any
time. This number follows a skewed distribution with a maximum at 120: as
expected for $x$ from the formula (BaH2)(H${}_{2})_{x}$. The radial
distribution function does not drop to zero, and so this is somewhat sensitive
to Rcut. With 368 hydrogen atoms in the cell, on average there are 133.5
unbonded H units compared to 128 expected in the relaxed ground state
structure. The simulation ran for 6.5ps during which time all hydrogen atoms
in the cell spent some time as molecules, most forming bonds at some point
with more than one of the four atoms in neighbouring tetrahedra. No hydrogen
atom ”escaped” from its original tetrahedron.

Heating to 700K with the equal mass sampling gives a similar result, while
continuing the simulation at 0GPa/700K reduces the rebonding and increases the
molecule count to an average of 119.

A longer run with 100GPa/1000K melted after a lengthy equilibration period.
The melt structure has a rather similar first peak to the solid H-H radial
distribution, showing that the majority of hydrogen is still in molecular
form: the exact number depends of the radius chosen for a molecules. A
distinct Ba-H rdf peak is evident around 2.2Å in both liquid and solid, but it
is less distinct in the liquid. The Ba-Ba liquid rdf has an unusual double
peak at 2.8Å and 3.9Å This is quite unlike the solid. We can associate the
first peak with adjacent Ba2+ ions (textbook ionic radius 1.35Å ), and the
second peak with two Ba’s separated by hydrogen. It is notable that, despite
the low Ba mass, the mean-squared displacement of the Ba is still much lower
than that of hydrogen, indicating the role of size in diffusion.

We attempted to determine stable ground-state structure by relaxing MD
snapshots. Relaxed structures were always significantly more stable than the
$Pm\bar{3}n$ structure reported from X-ray diffraction, although the Ba atoms
are always close to that structure (Fig.22).

Relaxing a snapshot from the large MD at 100GPa to 0GPa gave precisely 120
bonds by the sub 1Åcriterion, with a gap in observed bondlengths between 0.92
and 1.26Å, coinciding with a sharp drop in the Mulliken bond charge from 0.8
to 0.2.

Relaxing to 50GPa had a similar effect, with slightly longer bonds and a less
sharp division between bonded and unbonded. Again, the Mulliken charge drops
from 0.8 to 0.2 between 0.9Å and 1.2Å, but now there are several bonds with
distances (and charges) between.

Other snapshots showed similar results, but no unique structure emerged. It is
evident from the relaxed structures that the H2 molecules are oriented
perpendicular to the vector the short BaH vector. This follows from the fact
that the bond links hydrogens in different tetrahedra, and therefore goes
through a triangle of three bariums.

These unit cells are too large for lattice dynamics calculation of Raman
spectra, but we can get some indication of the vibrations in the system from
the velocity autocorrelation function (Sup. Fig.24 This shows a large bump
between 300 and 1800 cm-1 corresponding to lattice and librational motion, and
a broad shoulder extending up to 3300 cm-1. This breadth is due to both the
range of H2 environments, giving rise to multiple vibrational frequencies, and
their short lifetime which broadens each peak.

Figure 24:  Fourier transform of the velocity autocorrelation function from MD
of the BaH5.75 structure (100GPa/300K). The H2 vibrations show a broad but
unresolvable region between 2000 and 3000 wavenumbers.

###  5.4 S6.3 MD of BaH12

BaH12 was reported by Chen et al. to be a molecular metallic barium hydride,
identified and pseudocubic, and that ab initio calculations show that newly
discovered semimetallic BaH12 contains H2 and H3 molecular units and detached
H chains which are formed as a result of a Peierls-type distortion of the
cubic cage structure.

Figure 25: Radial distribution functions from molecular dynamics on BaH12.
Solid lines are total RDF, dashed lines count only the H-H distances. The
sharp first peak at 0.8Å indicates the preferred molecular bondlength seen in
these simulations and is almost independent of pressure. The compression with
pressure is evident only in the second peak - the intermolecular distance. The
second largest peak is the BaH nearest-neighbour separation, which can be seen
to reduce with pressure. The overall crystal structure can be described as
pseudocubic.  Figure 26: Radial distribution functions at all pressures
considered (300K) from molecular dynamics on BaH12.

Figure 27: BaH12 as BaH${}_{2}($H${}_{2})_{5}$ Relaxed P1 BaH12 structures at
50, 100 and 150GPa (left to right) Note the 100GPa case two Ba atoms have
moved across the periodic boundary conditions. Analysis of molecular dynamics
simulations from this structure reveals that the number of ”H2 bonds”
(criterion: atoms with one neighbour within 1Å ) is 5 per formula unit. This
conclusion is remarkably robust to reasonable choice of the bondlength. By
contrast, the more-easily calculated bond criterion ”number of H-H distances
less than 1Å ” is very sensitive to the choice of cutoff. This is because
there is a significant proton exchange via the
H2+H${}^{-}\rightarrow$H${}_{3}^{-}\rightarrow$H${}^{-}+$H2 mechanism.

We ran molecular dynamics at 300K/250GPa 300K/150GPa and 300K/50GPa and found
it to be stable (Sup. Figs. 25, 26 and 27 We used a supercell containing 8
formula units which in principle allows up to 48 molecules to form. For this
high hydrogen concentration bond making and breaking occurred frequently,
invariably via an H3 transient. Once again, the material is revealed as
BaH2(H${}_{2})_{x}$. Typical calculations with Ba8H96 generate precisely 40
molecules. The sharp H-H molecular peak is present at all pressures, becoming
notably less sharp, and developing a pronounced tail at higher pressures. The
second peak, corresponding to non-bonded H-H, is distinctive and moves to
shorter distances with pressure, as does the BaH separation.

There is some slow but significant diffusion of hydrogen (MSD at $\equiv 2\AA$
in 5ps). Interestingly, the diffusion is slowest at 50GPa, and fastest at
150GPa, with 250GPa being intermediate suggesting that both the rate of bond
breaking and the space available for rotation are important.

The mean squared displacement plateaus, as expected for a solid. The Ba atoms
move very little from their lattice sites (around 0.1Å ), despite the reduced
mass. However theBarium hydrogens are fairly mobile, with a MSD around 1.5Å .
This is consistent with making and breaking of bonds, plus molecular
rotations. Detailed analysis of the trajectories using vmd and inspection of
lists of neighbours shows a slow diffusion of hydrogens at a timescale of tens
of picoseconds associated with
H+H${}_{2}\rightarrow$H${}_{3}\rightarrow$H${}_{2}+$H processes followed by
molecular rotation.

We used the ’Atoms with precisely one neighbour within Rcut=1Å ’ criterion to
define molecules. If molecules were forming randomly, one would expect to see
a Gaussian distribution about some average value, with strong dependence on
Rcut.

In 10000 timesteps of MD with Ba8H96 at 150GPa/300K, 91% of snapshots showed
40 molecules, 8% showed 39 and 1% showed 38. Only at one timestep (0.01%) was
a 41 molecule configuration reported. The sub-40 molecule configurations
typically lasted for a few fs, indicating that they are a fast physical
process rather than noise. The precise number of non-40 molecule snapshots is
sensitive to the choice of Rcut, but the massive preponderance is of
BaH2(H${}_{2})_{5}$ configurations.

At 250GPa, the bond-breaking is more frequent, with 40 molecules appearing
only 16% of the time, but still 99% of snapshots show between between 36 and
40 neighbours.

In Chen’s paper they present results for fcc, I4/mmm, and Cmc21 structures,
although the most stable of their presented structure is pseudocubic P21. In
molecular dynamics calculation we find many snapshots with lower enthalpy, all
with P1 symmetry and BaH2(H${}_{2})_{5}$ character. H${}_{3}^{+}$ units are
observed transiently at high T, and can be captured in relaxed structures. We
find no evidence for longer chains or hydrogen molecular units with more than
three atoms.

##  6 S7 Relation to superconductivity

There has been a tremendous interest in high temperature superconductivity in
polyhydrides. This follows from Ashcroft’s observation that metallic hydrogen
could become a high-Tc superconductor14. The pressures required to make
metallic hydrogen mean this has yet to be demonstrated. However, many
hydrogen-rich compounds have shown high Tc superconductivity at achievable
pressures. When hydrogen is mixed with a cation which donates electrons either
to an atomic hydrogen conduction band (favourable) or a molecular antibonding
orbital (unfavourable). In this way the molecular-atomic transition moves to
lower pressures. However, as in BaH4, the atomic form is much denser than the
molecular form, so pressure is an essential stabilising factor via the PV term
in the Gibbs free energy.

The rationale for high Tc polyhydride superconductors is that the temperature
required to break the electron-phonon coupling is related to the phonon
energy, and light hydrogen atoms leads to high phonon frequency and energy.

This leads to a dichotomy at low pressures. The hydrogens need to be close
enough to form an extended phonon, but not so close that they form H2
molecules. Another is that the atomic hydrogens in the conduction band are
spaced by the hydrogens in molecules. These molecules nevertheless play a
significant role in the Cooper-pair forming lattice phonons. BaH4 illustrates
this mechanism: here the molecules are H2 but other spacer molecules are
possible provided their vibration couples to the lattice phonon.

A final consideration is that the molecule formation is an alternative
electron-pairing mechanism to Cooper-pair formation. In both cases the
associated band gap opening stabilises the structure. In our Born-Oppenheimer
MD there is no quantum coupling between electronic and nuclear degrees of
freedom, so if the superconducting state was lower in energy than the
molecularised one, we would not know.

## References

  * Drozdov et al. 2015 Drozdov, A. P.; Eremets, M. I.; Troyan, I. A.; Ksenofontov, V.; Shylin, S. I. Conventional superconductivity at 203 K at high pressures. _Nature_ 2015, _525_ , 1–19
  * Drozdov et al. 2018 Drozdov, A.; Minkov, V.; Besedin, S.; Kong, P.; Kuzovnikov, M.; Knyazev, D.; Eremets, M. Superconductivity at 215 K in lanthanum hydride at high pressures. _arXiv preprint arXiv:1808.07039_ 2018,
  * Drozdov et al. 2019 Drozdov, A.; Kong, P.; Minkov, V.; Besedin, S.; Kuzovnikov, M.; Mozaffari, S.; Balicas, L.; Balakirev, F.; Graf, D.; Prakapenka, V., et al. Superconductivity at 250 K in lanthanum hydride under high pressures. _Nature_ 2019, _569_ , 528
  * Nakao et al. 2019 Nakao, H.; Einaga, M.; Sakata, M.; Kitagaki, M.; Shimizu, K.; Kawaguchi, S.; Hirao, N.; Ohishi, Y. Superconductivity of pure H3S synthesized from elemental sulfur and hydrogen. _Journal of the Physical Society of Japan_ 2019, _88_ , 123701
  * Osmond et al. 2022 Osmond, I.; Moulding, O.; Cross, S.; Muramatsu, T.; Brooks, A.; Lord, O.; Fedotenko, T.; Buhot, J.; Friedemann, S. Clean-limit superconductivity in I m 3¯ m H 3 S synthesized from sulfur and hydrogen donor ammonia borane. _Physical Review B_ 2022, _105_ , L220502
  * Roman et al. 2022 Roman, L.; Lilia, B.; Christoph, H., et al. In-silico synthesis of lowest-pressure high-T c ternary superhydrides. _NPJ Computational Materials_ 2022, _8_
  * Errea 2022 Errea, I. Superconducting hydrides on a quantum landscape. _Journal of Physics: Condensed Matter_ 2022, _34_ , 231501
  * Errea et al. 2020 Errea, I.; Belli, F.; Monacelli, L.; Sanna, A.; Koretsune, T.; Tadano, T.; Bianco, R.; Calandra, M.; Arita, R.; Mauri, F., et al. Quantum crystal structure in the 250-kelvin superconducting lanthanum hydride. _Nature_ 2020, _578_ , 66–69
  * Pickard et al. 2020 Pickard, C. J.; Errea, I.; Eremets, M. I. Superconducting hydrides under pressure. _Annual Review of Condensed Matter Physics_ 2020, _11_ , 57–76
  * Zhang et al. 2022 Zhang, X.; Zhao, Y.; Yang, G. Superconducting ternary hydrides under high pressure. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2022, _12_ , e1582
  * Duan et al. 2017 Duan, D.; Liu, Y.; Ma, Y.; Shao, Z.; Liu, B.; Cui, T. Structure and superconductivity of hydrides at high pressures. _National Science Review_ 2017, _4_ , 121–135
  * Flores-Livas et al. 2020 Flores-Livas, J. A.; Boeri, L.; Sanna, A.; Profeta, G.; Arita, R.; Eremets, M. A perspective on conventional high-temperature superconductors at high pressure: Methods and materials. _Physics Reports_ 2020, _856_ , 1–78
  * Sun et al. 2020 Sun, W.; Kuang, X.; Keen, H. D.; Lu, C.; Hermann, A. Second group of high-pressure high-temperature lanthanide polyhydride superconductors. _Physical Review B_ 2020, _102_ , 144524
  * Ashcroft 1968 Ashcroft, N. W. Metallic hydrogen: A high-temperature superconductor? _Phys. Rev. Lett._ 1968, _21_ , 1748
  * Ashcroft 2004 Ashcroft, N. W. Hydrogen Dominant Metallic Alloys : High Temperature Superconductors ? 2004, 1–4
  * Zurek et al. 2009 Zurek, E.; Hoffmann, R.; Ashcroft, N.; Oganov, A. R.; Lyakhov, A. O. A little bit of lithium does a lot for hydrogen. _Proceedings of the National Academy of Sciences_ 2009, _106_ , 17640–17643
  * Pena-Alvarez et al. 2020 Pena-Alvarez, M.; Li, B.; Kelsall, L. C.; Binns, J.; Dalladay-Simpson, P.; Hermann, A.; Howie, R. T.; Gregoryanz, E. Synthesis of Superconducting Cobalt Trihydride. _The journal of physical chemistry letters_ 2020, _11_ , 6420–6425
  * Peña-Alvarez et al. 2022 Peña-Alvarez, M.; Binns, J.; Marqués, M.; Kuzovnikov, M. A.; Dalladay-Simpson, P.; Pickard, C. J.; Ackland, G. J.; Gregoryanz, E.; Howie, R. T. Chemically Assisted Precompression of Hydrogen Molecules in Alkaline-Earth Tetrahydrides. _The Journal of Physical Chemistry Letters_ 2022, _13_ , 8447–8454
  * Chen et al. 2021 Chen, W.; Semenok, D. V.; Kvashnin, A. G.; Huang, X.; Kruglov, I. A.; Galasso, M.; Song, H.; Duan, D.; Goncharov, A. F.; Prakapenka, V. B., et al. Synthesis of molecular metallic barium superhydride: pseudocubic BaH12. _Nature Commun._ 2021, _12_ , 1–9
  * Van Kranendonk 1983 Van Kranendonk, J. _Solid Hydrogen: Theory of the Properties of solid H2, HD, and D2_ ; Springer, 1983
  * Magdău and Ackland 2013 Magdău, I. B.; Ackland, G. J. Identification of high-pressure phases III and IV in hydrogen: Simulating Raman spectra using molecular dynamics. _Phys. Rev. B_ 2013, _87_ , 174110
  * Ackland and Loveday 2020 Ackland, G. J.; Loveday, J. S. Structures of solid hydrogen at 300 K. _Phys. Rev. B_ 2020, _101_ , 094104
  * Howie et al. 2012 Howie, R. T.; Guillaume, C. L.; Scheler, T.; Goncharov, A. F.; Gregoryanz, E. Mixed molecular and atomic phase of dense hydrogen. _Physical Review Letters_ 2012, _108_ , 125501
  * Tse et al. 2007 Tse, J.; Klug, D.; Desgreniers, S.; Smith, J.; Flacau, R.; Liu, Z.; Hu, J.; Chen, N.; Jiang, D. Structural phase transition in CaH2 at high pressures. _Phys. Rev. B_ 2007, _75_ , 134108
  * Smith et al. 2009 Smith, J. S.; Desgreniers, S.; Klug, D. D.; John, S. T. High-density strontium hydride: An experimental and theoretical study. _Solid state communications_ 2009, _149_ , 830–834
  * Kinoshita et al. 2007 Kinoshita, K.; Nishimura, M.; Akahama, Y.; Kawamura, H. Pressure-induced phase transition of BaH2: Post Ni2In phase. _Solid state communications_ 2007, _141_ , 69–72
  * Smith et al. 2007 Smith, J. S.; Desgreniers, S.; Tse, J. S.; Klug, D. D.; Smith, J. S.; Desgreniers, S.; Tse, J. S.; Klug, D. D. High-pressure phase transition observed in barium hydride. _Journal of Applied Physics_ 2007, _043520_
  * John et al. 2009 John, S. T.; Song, Z.; Yao, Y.; Smith, J. S.; Desgreniers, S.; Klug, D. D. Structure and electronic properties of BaH2 at high pressure. _Solid State Communications_ 2009, _149_ , 1944–1946
  * Bi and Zurek 2021 Bi, T.; Zurek, E. Electronic structure and superconductivity of compressed metal tetrahydrides. _Chem. Eur. J._ 2021, _27_ , 14858–14870
  * Mishra et al. 2018 Mishra, A. K.; Muramatsu, T.; Liu, H.; Geballe, Z. M.; Somayazulu, M.; Ahart, M.; Baldini, M.; Meng, Y.; Zurek, E.; Hemley, R. J. New Calcium Hydrides with Mixed Atomic and Molecular Hydrogen. _The Journal of Physical Chemistry C_ 2018, _122_ , 19370–19378
  * Wu et al. 2019 Wu, G.; Huang, X.; Xie, H.; Li, X.; Liu, M.; Liang, Y.; Duan, D.; Li, F.; Liu, B.; Cui, T. Unexpected calcium polyhydride CaH4 : A possible route to dissociation of hydrogen molecules Unexpected calcium polyhydride CaH 4 : A possible route to dissociation of hydrogen molecules. _J. Chem. Phys._ 2019, _150_ , 044507
  * Hooper et al. 2014 Hooper, J.; Terpstra, T.; Shamp, A.; Zurek, E. Composition and constitution of compressed strontium polyhydrides. _The Journal of Physical Chemistry C_ 2014, _118_ , 6433–6447
  * Peña-Alvarez et al. 2021 Peña-Alvarez, M.; Binns, J.; Martinez-Canales, M.; Monserrat, B.; Ackland, G. J.; Dalladay-Simpson, P.; Howie, R. T.; Pickard, C. J.; Gregoryanz, E. Synthesis of Weaire–Phelan Barium Polyhydride. _The journal of physical chemistry letters_ 2021, _12_ , 4910–4916
  * Bader 1994 Bader, R. F. W. Atoms in Molecules: A Quantum Thory. 1994,
  * Bonev and Ashcroft 2001 Bonev, S.; Ashcroft, N. Hydrogen in jellium: First-principles pair interactions. _Physical Review B_ 2001, _64_ , 224112
  * Ackland and Magdau 2014 Ackland, G. J.; Magdau, I. B. Efficacious calculation of raman spectra in high pressure hydrogen. _High Pressure Research_ 2014, _34_ , 198–204
  * Cooke et al. 2020 Cooke, P. I.; Magdău, I. B.; Peña-Alvarez, M.; Afonina, V.; Dalladay-Simpson, P.; Liu, X.-D.; Howie, R. T.; Gregoryanz, E.; Ackland, G. J. Raman signal from a hindered hydrogen rotor. _Phys. Rev. B_ 2020, _102_ , 064102
  * Cooke et al. 2022 Cooke, P. I.; Magdău, I. B.; Ackland, G. J. Calculating the Raman Signal Beyond Perturbation Theory for a Diatomic Molecular Crystal. _arXiv preprint arXiv:2202.09604_ 2022,
  * Bi and Zurek 2021 Bi, T.; Zurek, E. Electronic structure and superconductivity of compressed metal tetrahydrides. _Chemistry–A European Journal_ 2021, _27_ , 14858–14870
  * Ackland 2016 Ackland, G. J. Rapid Equilibration by algorithmic quenching the ringing mode in molecular dynamics. _MRS Advances_ 2016, _1_ , 2857–2865
  * Momma and Izumi 2011 Momma, K.; Izumi, F. VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data. _Journal of applied crystallography_ 2011, _44_ , 1272–1276
  * Humphrey et al. 1996 Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. _Journal of molecular graphics_ 1996, _14_ , 33–38
  * Kinoshita et al. 2007 Kinoshita, K.; Nishimura, M.; Akahama, Y.; Kawamura, H. Pressure-induced phase transition of BaH2: Post Ni2In phase. _Solid state communications_ 2007, _141_ , 69–72
  * Smith et al. 2007 Smith, J. S.; Desgreniers, S.; Tse, J. S.; Klug, D. D. High-pressure phase transition observed in barium hydride. _Journal of Applied Physics_ 2007, _102_ , 043520
  * Zhang et al. 2018 Zhang, X.; Wang, X.; Wang, Q.; Ma, X.; Liu, C.; Li, P.; Liu, C.; Han, Y.; Ma, Y.; Gao, C. Hydride ion (H-) transport behavior in barium hydride under high pressure. _Physical Chemistry Chemical Physics_ 2018, _20_ , 8917–8923
  * Novak et al. 2020 Novak, E.; Haberl, B.; Daemen, L.; Molaison, J.; Egami, T.; Jalarvo, N. Pressure-induced phase transition in barium hydride studied with neutron scattering. _Applied Physics Letters_ 2020, _117_ , 051902
  * Novak et al. 2021 Novak, E.; Daemen, L.; Page, K.; Neuefeind, J.; Everett, M.; Egami, T.; Jalarvo, N. Temperature dependent local atomic structure and vibrational dynamics of barium hydride and calcium hydride. _The Journal of Physical Chemistry C_ 2021, _125_ , 24328–24339
  * Becke and Edgecombe 1990 Becke, A. D.; Edgecombe, K. E. A simple measure of electron localization in atomic and molecular systems. _J. Chem. Phys._ 1990, _92_ , 5397–5403
  * Savin et al. 1997 Savin, A.; Nesper, R.; Wengert, S.; Fässler, T. F. ELF: The electron localization function. _Angewandte Chemie International Edition in English_ 1997, _36_ , 1808–1832
  * Clark et al. 2005 Clark, S. J.; Segall, M. D.; Pickard, C. J.; Hasnip, P. J.; Probert, M. I.; Refson, K.; Payne, M. C. First principles methods using CASTEP. _Zeitschrift für kristallographie-crystalline materials_ 2005, _220_ , 567–570
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. . _Phys. Rev. Lett._ 1996, _77_ , 3865
  * Perdew and Zunger 1981 Perdew, J. P.; Zunger, A. Self-interaction correction to density-functional approximations for many-electron systems. _Physical Review B_ 1981, _23_ , 5048
  * Perdew and Zunger 1981 Perdew, J. P.; Zunger, A. Self-interaction correction to density-functional approximations for many-electron systems. _Physical Review B_ 1981, _23_ , 5048
  * Becke 1988 Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. _Physical review A_ 1988, _38_ , 3098
  * Lee et al. 1988 Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. _Physical review B_ 1988, _37_ , 785
  * Otero-de-la Roza et al. 2014 Otero-de-la Roza, A.; Johnson, E.; Luana, V. Critic2: A program for real-space analysis of quantum chemical interactions in solids. _Comput. Phys. Commun._ 2014, _185_ , 1007
  * Mostofi et al. 2014 Mostofi, A. A.; Yates, J. R.; Pizzi, G.; Lee, Y.-S.; Souza, I.; Vanderbilt, D.; Marzari, N. An updated version of wannier90: A tool for obtaining maximally-localised Wannier functions. _Computer Physics Communications_ 2014, _185_ , 2309–2310
  * Ahuja et al. 1998 Ahuja, R.; Eriksson, O.; Wills, J. M.; Johansson, B. Theoretical high-pressure studies of caesium hydride. _Journal of Physics: Condensed Matter_ 1998, _10_ , L153
  * Magdău and Ackland 2017 Magdău, I. B.; Ackland, G. J. Charge density wave in hydrogen at high pressure. Journal of Physics: Conference Series. 2017; p 042058