Model system. The bulk structure of β-NiOOH used in our simulations is based on our previous study . It is represented as a model system with a 1 × 2 × 1 unit cell, the stability of which compares well with the most stable bulk phases reported in the literature . The material has the ABBCC oxygen stacking sequence , for which hydrogen atoms present in the inter-layer space create a saturated number of hydrogen bonds between the layers, in agreement with infrared spectroscopy measurements . Surface structures were modelled as three-layer slabs of β-NiOOH(0001) with 3 × 4 unit cells. Stoichiometric surface structures include 12 Ni atoms, 24 O atoms, as well as 12 H atoms per layer of the slab for β-NiOOH(0001). These unit cell configurations are commensurate with the addition of explicit water layers. A water coverage of 5/6 ML at NiOOH(0001), initially equilibrated with ab initio molecular dynamics simulation at 140 K for 8 ps (1 fs time-step), was used . The H-down configuration of water was found to be more stable on NiOOH(0001) due to formation of H-bonds between water and unsaturated surface O atoms. The (0001) facet of NiOOH is terminated with 50% of O* and 50% OH*, as shown in Fig. . Surface deprotonation was modelled by successive removal of surface H atoms from the surface, while keeping the unit cell electroneutral at any given degree of deprotonation. Thereafter, adsorption of Li + , Na + , K + and Cs + ions on the fully deprotonated surface was considered, as shown in Fig. (b). To find the most stable adsorption site for cations, all possible surface O adsorption sites on the deprotonated surface -including all atop, bridge and the three-fold sites -were examined, from which the cations were found to be stabilized only on the four possible three-fold sites. Among these sites, site α, as shown in Fig. , was found to be the most stable site for all cations, which is on top of the surface Ni atom and involves two of the deprotonated O atoms, as shown in Fig. . The relative stability of these four sites is provided in Tab. 1. DFT+U calculations. Spin-polarized DFT-based calculations were performed with the Vienna Ab initio Simulation Package (VASP) . The ionic cores were represented by projector augmented waves (PAW) . Kohn-Sham one-electron wave functions were expanded in a plane wave basis set with an energy cutoff of 500 eV. Electronic exchange and correlation effects were incorporated within the generalized gradient approximation (GGA), using the PBE exchange-correlation functional . The DFT+U method has been employed to account for the strong correlations of d-electrons of Ni atoms. We used the approach by Dudarev et al. with an effective Hubbard U eff parameter (U-J) 105 of 5.5 eV, which was calculated by Li and Selloni using linear response theory . A geometry relaxation for bulk NiOOH was performed with a 12×12×8