## 2 Methodology  
Born-Oppenheimer (BO)-MD simulations were performed under periodic boundary conditions (PBC) using a vdW-DF dispersion functional based on generalized gradient approximation (GGA) employing Dion-Rydberg-Schroder-Langreth-Lundqvist (DRSLL) non-local correlation correction [56; 57]. This functional was recently shown to perform exceptionally well for liquid water [58; 59; 60; 61; 62; 63] and has also yielded realistic results for the H\({}_{2}\)-hopping free-energy barrier in H\({}_{2}\) hydrates [3] and Raman vibrational properties vis-a-vis experimental results [45]. The simulations were performed with a time step of 0.2 fs in a canonical (NVT) ensemble. The system temperature was imposed at 258 K using a Nose-Hoover thermostat [64; 65] with fictitious mass 120 Ry.fs2. The charge density, described in the DZP basis set, converged to an accuracy of 10\({}^{-4}\) a.u at each MD step. The calculations were performed using the electronic structure software package SIESTA [66]. The ab-initio (AI)-MD trajectories were run for 5 ps at 258 K on a single unit cell in the NVT ensemble for large cage occupation numbers of \(n\) = 1,..., 4 and single occupation of small cages for each and every cage. In addition, the case of an almost-empty lattice with a solitary H\({}_{2}\) molecule in just one large cage was also simulated at 258 K, with the lattice both entirely free and held rigid, so as to determine if H-H peak-splitting occurs for just a single H\({}_{2}\) molecule (eliminating H\({}_{2}\)-H\({}_{2}\) interactions in multiply occupied cages, and even those weaker ones from neighboring cages) and the role, if any, that a mobile, vibrating lattice may have on this. Each structure was first run in the NPT ensemble at 258 K and 2 kbar beforehand to equilibrate to experimentally relevant conditions [35], producing equilibrium cell lengths of 16.5-16.6 A. In the present study, all calculations, including those with empirical models (vide infra), were performed with a 16.6 A cell length, consistent with DFT results. Footnote 2: The \(\hat{H}\)-H\({}_{2}\) molecule is a \(\hat{H}\) molecule, with  
where \(\chi_{i}\) terms have analytical solutions. Ideally, the reference Hamiltonian is chosen to be not too dissimilar from the original Hamiltonian \(\hat{H}\). We can use the states of the reference Hamiltonian as an orthonormal basis in which to expand  
\[\sum_{j}\langle\chi_{j}|\hat{H}|\chi_{j}\rangle\ \langle\chi_{j}|\psi_{\rangle} \rangle=\epsilon\langle\chi_{j}|\psi_{\rangle}\rangle \tag{3}\]  
Indeed, defining  
\[H_{i,j}=\left\langle\chi_{j}|\hat{H}|\chi_{j}\right\rangle\text{ and }a_{n,i}=\left\langle\chi_{i}|\psi_{ \rangle}\right\rangle\]  
we have  
\[\sum_{j}H_{i,j}a_{n,j}=\epsilon\langle a_{n,j} \tag{4}\]  
Here, the eigenvalues are the energy levels; the eigenvectors can be used to reconstruct the eigenfunctions associated with each energy level. In practice, the \(H_{i,j}\) matrix elements are calculated up to a reasonable maximum value of \(i,j\), which can then be diagonalized to obtain the eigenvectors and eigenvalues. For H-H stretch modes, a natural basis choice is provided by the eigenfunctions of the quantum harmonic oscillator, which, in the \(\chi\) representation, are given by  
\[\chi_{n}(\Delta x)=\frac{1}{2^{n}n!}\Big{(}\frac{\alpha}{n}\Big{)}^{\frac{1}{ 2}}\exp\!\left(-\frac{\alpha\Delta x^{2}}{2}\right)H_{n}\Big{(}\sqrt{\alpha} \Delta x\Big{)}. \tag{5}\]  
For \(n=0\), \(1\), \(2\),... with associated eigenvalues \(\chi_{n}=(n+1/2)h\omega_{0}\), where \(\Delta x\) is the displacement from the minimum \((\Delta x=x-x_{0})\) and \(\alpha=h\omega_{0}/m\star\), where m* is the reduced mass and \(\omega_{0}\) is the harmonic oscillator frequency, chosen at the outset to be close to the curvature of the potential minimum.