Methods. All calculations were performed using density functional theory (DFT) with Perdew-Burke-Ernzerhof (PBE) functional and D3 correction of van der Waals interactions following the Grimme approach as implemented in the CP2K 3.0 package. The Quickstep method was employed, with Goedecker-Teter-Hutter (GTH) pseudopotentials together with DZVP-MOLOPT-GTH-SR basis set for Mo, S, B, C, and N. Hydrogen was treated with the DZVP-all electron basis set and all-electron potential. The plane-wave energy cutoff was set to 360 Ry. All atomic positions and lattice parameters were fully optimized, the optimized parameters for all the systems are given in Table . For more details see Methods in SI.