For the CVs defined in each system, the following reference H-X (X = N, C, B, or S) distances (R 0 ) were used: 2.3 bohr for H + @hBN, 2.5 bohr for H@graphite, 2.9 bohr for H@hBN, 3.0 bohr for H@MoS 2 and H + @ MoS 2 . These parameters are slighltly larger (up to 10%) than the average value in the fluctuation of the corresponding H-X distances during the BOMD simulations. Thus, when the H or H + is close to the transition state region, the coordination number of H to X is close to 0. Note that in the present models, WTMetaD setup, and the definition of the CVs, the hopping between the layers can be discussed, however, the direction of the transport is not defined and is governed by Brownian motion. To understand the latter, more CVs would need to be defined, e.g., the so-called path CVs, in several different directions, in order to sample the whole free-energy landscape (this is beyond the scope of this work). We have also investigated the case with three CVs specified, where the third one describes the hopping of H species within the same layer. However, such a process was not observed with the employed number of hills, thus in all cases, the free energy barrier for transport within the same layer was much higher than between the layers. Thus, it is a safe choice to reduce the number of CVs to just the two discussed above. Static calculations: Layered materials exhibit low-energy shearing modes, therefore, we have investigated the potential energy barriers for H + transfer between the layer in two different layer stackings of h-BN, namely AA' (N over B) and A'B stackings (N over N). Such exemplary static calculations reveal strong energy barrier reduction when changing the stacking of the layers. For the AA' stacking, the barrier is as large as 0.65 eV, whereas for the A'B stacking is only 0.17 eV, supporting our conclusion about the shear-assisted transport of H species between layers of the studied systems. For MoS 2 , the shear mode also exists, but the transfer barrier is nevertheless very low and in the static simulations, we have obtained 0.16 eV, see Binding Energies: The binding energies were calculated as follows (explained on the example of h-BN with H + ): E b = E H+@h-BN -E h-BN -E H+ , with E h-BN and E H+ being the energies of the perfect h-BN and proton, respectively, and E H+@h-BN being the total energy of the combined structure of H + @h-BN. Due to the difficulty to calculate total energies of charged solid-state materials, we have used cluster models. For this purpose, we used Gaussian 09 software 10 to calculate the binding energy of proton in h-BN cluster with different flake sizes and number of layers (pyrene and coronene cluster models with 1-4 layers). These calculations were performed using B3LYP functional 11 with 6-311+g (d, p) basis set, both with and without zero point energy correction (see Figure ). We found that the binding energies corrected with zero-point energy converge to -8.55 eV in pyrene and -8.61 eV in coronene cluster models.