An anisotropic interlayer force field that describes the interlayer interactions in molybdenum disulfide (MoS₂) is presented. The force field is benchmarked against density functional theory calculations for both bilayer and bulk systems within the Heyd-Scuseria-Ernzerhof hybrid density functional approximation, augmented by a non-local many-body dispersion treatment of long-range correlation. The parametrization yields good agreement with reference calculations of binding energy curves and sliding potential energy surfaces for both bilayer and bulk configurations. Benchmark calculations for the phonon spectra of bulk MoS₂ provide good agreement with experimental data and the calculated bulk modulus falls in the lower part of experimentally measured values. This indicates the accuracy of the interlayer force field near equilibrium. Under external pressures up to 20 GPa, the developed force field provides a good description of compression curves. At higher pressures, deviations from experimental data grow, signifying the validity range of the developed force field.