An anisotropic interlayer force field that describes the interlayer interactions in homogeneous and heterogeneous interfaces of group-VI transition metal dichalcogenides (MX₂, where M = Mo, W, and X = S, Se) is presented. The force field is benchmarked against density functional theory calculations for bilayer systems within the Heyd-Scuseria-Ernzerhof hybrid density functional approximation, augmented by a nonlocal many-body dispersion treatment of long-range correlation. The parametrization yields good agreement with the reference calculations of binding energy curves and sliding potential energy surfaces. It is found to be transferable to transition metal dichalcogenide (TMD) junctions outside of the training set that contain the same atom types. Calculated bulk moduli agree with most previous dispersion-corrected density functional theory predictions, which underestimate the available experimental values. Calculated phonon spectra of the various junctions under consideration demonstrate the importance of appropriately treating the anisotropic nature of the layered interfaces. Considering our previous parametrization for MoS₂, the anisotropic interlayer potential enables accurate and efficient large-scale simulations of the dynamical, tribological, and thermal transport properties of a large set of homogeneous and heterogeneous TMD interfaces.