An anisotropic interlayer potential (ILP), designed to describe the interlayer interaction in graphene/MX₂ and h-BN/MX₂ (M = Mo, W; X = S, Se) heterostructures, is presented. The ILP is parametrized against density functional theory (DFT) calculations within the Perdew-Burke-Ernzerhof (PBE) generalized-gradient approximation, augmented by nonlocal many-body dispersive (MBD-NL) interactions. The parametrized force field demonstrates excellent agreement with the DFT reference data of binding energy curves and sliding energy surfaces across all heterostructures considered. The transferability of the developed ILP is demonstrated for the phenalenyl (C₁₃H₉)/MoS₂ and B₇N₆H₉/MoS₂ interfaces, which are outside the training set. The force field is then used to study equilibrium interlayer distances, bulk moduli, and phonon spectra by means of molecular dynamics simulations.