  Here, a magnetic model with an unprecedentedly large number of parameters was
determined from first-principles calculations for transition-metal phosphorus
trisulfides (TMPS$_3$'s), which reproduced the measured magnetic ground states
of bulk TMPS$_3$'s. Our Monte Carlo simulations for the critical temperature,
magnetic susceptibility, and specific heat of bulk and few-layer TMPS$_3$'s
agree well with available experimental data and show that the antiferromagnetic
order of FePS$_3$ and NiPS$_3$ persists down to monolayers. Remarkably, the
orbital polarization, which was neglected in recent first-principles studies,
dramatically enhances the magnetic anisotropy of FePS$_3$ by almost two orders
of magnitude. A recent Raman study [K. Kim et al., Nat. Commun. 10, 345 (2019)]
claimed that magnetic ordering is absent in monolayer NiPS$_3$ but
simultaneously reported a strong two-magnon continuum; we show that the
criterion used to judge magnetic ordering there is invalid in monolayer
NiPS$_3$, thus providing an understanding of the two seemingly contradictory
experimental results. The rich predictions on the magnetic susceptibility and
specific heat of few-layer FePS$_3$ and NiPS$_3$ await immediate experimental
verifications.
