  In this study, we use Hubbard-Corrected density functional theory (DFT+$U$)
to derive spin model Hamiltonians consisting of Heisenberg exchange
interactions up to the fourth nearest neighbors and bi-quadratic interactions.
We map the DFT+$U$ results of several magnetic configurations to the Heisenberg
spin model Hamiltonian to estimate Heisenberg exchanges. We demonstrate that
the number of magnetic configurations should be at least twice the number of
exchange parameters to estimate exchange parameters correctly. To calculate
biquadratic interaction, we propose specific non-collinear magnetic
configurations that do not change the energy of the Heisenberg spin model. We
use classical Monte Carlo (MC) simulations to evaluate DFT+$U$ results. We
obtain the temperature dependence of magnetic susceptibility and specific heat
to determine the Curie-Weiss and N\'eel temperatures. The MC simulations reveal
that although the biquadratic interaction can not change the N\'eel
temperature, it modifies the order parameter. We indicate that for a fair
comparison between classical MC simulations and experiments, we need to
consider the quantum effect by applying $(S+1)/S$ correction in classical MC
simulations.
