  Finding the true magnetic structure at given external conditions is crucial
for describing magnetic materials and predicting their properties. This is
especially important for high-throughput screening of potentially good magnets
that without adequate description of the magnetic structure may result in an
enormous waste of computational resources. I introduce a method, which
accurately and efficiently treats magnetic ordering in the course of standard
first-principles calculations. The method is suitable for all temperatures and
is based on Monte Carlo (MC) technique. At high temperatures MC is coupled to
ab initio molecular dynamics, therefore treating atomic vibrations and magnetic
ordering simultaneously. The method takes care of the convergence to the true
ground-state magnetic structure at 0 K, coupling between spins and atomic
vibrations at high temperatures, and spin-spin interactions in non-Heisenberg
systems, i.e. it naturally goes beyond the usual pair-like effective exchange
interactions at all temperatures. I show that the method nicely reproduces the
magnetic structures at low and high temperatures even on relatively small
supercells and provides a ground for truly ab initio calculations of magnetic
systems in high-throughput studies. In its general formulation the MDMC method
deals with non-collinear magnetism, but its collinear version may be more
suitable for fast high-throughput calculations. The examples of the application
of the MDMC method addressed in this paper are $\alpha$-Mn; body-centered cubic
Fe, which is ferromagnetic at low and paramagnetic at high temperatures; and
MnB$_2$W$_2$, a new magnetic compound with non-trivial magnetism.
