All DFT calculations were performed with the Vienna ab initio Simulation Package (VASP)[1]}, [2]}, [3]} using projector-augmented wave function (PAW) pseudopotentials[4]}, [5]} and the Perdew-Burke-Ernzerhof exchange correlation functional (PBE)[6]}. We found our results to be converged with a planewave cutoff of 400 eV and an energy convergence criteria of \(10^{-6}\)  eV, both of which were used for the results reported here. Fourth order Methfessel-Paxton smearing[7]} was used with a value of 0.13 eV for all geometry and cell lattice optimizations in order to ensure energy convergence without dependence on the electronic smearing temperature. The Mermin functional[8]} with the same electronic temperature was used for all MD calculations in order to avoid spurious forces due to possible negative occupation numbers from the Methfessel-Paxton approach. Brillouin Zone sampling for all TiH\(_2\)  unit cell calculations was performed with a \(10\times 10\times 10\)  k-point mesh, whereas we used a mesh of \(5\times 5\times 5\)  for 32 formula unit (96 atom) bulk calculations. We used system sizes of 168 atoms/7 layers for the (001) surface, 144 atoms/6 layers for the (011) surface, and 192 atoms/8 layers for the (111) surface, each with a vacuum of 20\(\,\) Å\(\,\)  and a k-point mesh of \(5\times 5\times 1\)  in the direction of the surface.
