  Thermal properties of molybdenum disulfide (MoS$_2$) have recently attracted
attention related to fundamentals of heat propagation in strongly anisotropic
materials, and in the context of potential applications to optoelectronics and
thermoelectrics. Multiple empirical potentials have been developed for
classical molecular dynamics (MD) simulations of this material, but it has been
unclear which provides the most realistic results. Here, we calculate lattice
thermal conductivity of single- and multi-layer pristine MoS$_2$ by employing
three different thermal transport MD methods: equilibrium, nonequilibrium, and
homogeneous nonequilibrium ones. These methods allow us to verify the
consistency of our results and also facilitate comparisons with previous works,
where different schemes have been adopted. Our results using variants of the
Stillinger-Weber potential are at odds with some previous ones and we analyze
the possible origins of the discrepancies in detail. We show that, among the
potentials considered here, the reactive empirical bond order (REBO) potential
gives the most reasonable predictions of thermal transport properties as
compared to experimental data. With the REBO potential, we further find that
isotope scattering has only a small effect on thermal conduction in MoS$_2$ and
the in-plane thermal conductivity decreases with increasing layer number and
saturates beyond about three layers. We identify the REBO potential as a
transferable empirical potential for MD simulations of MoS$_2$ which can be
used to study thermal transport properties in more complicated situations such
as in systems containing defects or engineered nanoscale features. This work
establishes a firm foundation for understanding heat transport properties of
MoS$_2$ using MD simulations.
