Jaime Garrido Aldea
ab,
Dorye L. Esteras
*a,
Stephan Roche
ac and
Jose H. Garcia
a
aCatalan Institute of Nanoscience and Nanotechnology (ICN2)CSIC and BIST Campus UAB, Bellaterra, Barcelona 08193, Spain. E-mail: jaime.garrido@icn2.cat
bUniversitat Autànoma de Barcelona (UAB), 08193 Barcelona, Spain
cInstitució Catalana de Recerca i Estudis Avançats (ICREA), Spain
First published on 7th July 2025
The discovery of intrinsic magnetism in two-dimensional (2D) materials has opened new frontiers in material science and technology. This review offers a detailed guide to modeling 2D magnetic materials using Density Functional Theory (DFT), focusing on both fundamental concepts and practical methodologies. Starting with the principles of magnetism, it examines the unique challenges of 2D systems, including the effects of anisotropy in stabilizing magnetic order, the limitations imposed by the Mermin–Wagner theorem, and the critical role of exchange interactions. The review introduces DFT basics, highlighting approaches to address electron delocalization through methods like DFT+U and hybrid functionals, and emphasizes the importance of incorporating van der Waals corrections for layered systems. Strategies for determining ground-state spin configurations for both collinear and non-collinear arrangements, are discussed, alongside advanced techniques like spin-constrained DFT and the Generalized Bloch Theorem for spin-spiral states. Methods for extracting magnetic exchange parameters and estimating critical temperatures from first-principles calculations are comprehensively covered. Practical insights are provided for applying these techniques to explore material databases and identify 2D magnets with promising properties for room-temperature applications. This review serves as a resource for theoretical and computational studies of 2D magnetic materials.
2D magnetic materials have remained so elusive all this time as a consequence of the Mermin–Wagner theorem,4 known since 1966. This theorem states that any 2D material with infinite size and with a continuous symmetry cannot present long range magnetic order at nonzero temperature when considering an isotropic Heisenberg model with finite range exchange interactions. Fortunately, this theorem does not contemplate anisotropy, and we know as of today that introducing single ion anisotropy or exchange anisotropy can stabilize magnetic order in 2D at nonzero temperature.5,6
These groundbreaking discoveries have positioned 2D magnetic materials as pivotal components for emerging technologies, including spintronics,7,8 magnetic sensors,9 energy harvesting systems and green energy applications10,11 and nonvolatile magnetic memories.12 Moreover, their 2D nature introduces unique opportunities, such as stacking and twisting, to tailor properties for specific applications.13
For practical device applications, achieving magnetic stability above room temperature remains a crucial challenge. Efforts are underway to identify 2D magnets with high intrinsic Tc.5,14–19 While various techniques exist to enhance magnetic stability,20 finding materials with intrinsic stability is paramount, as many enhancement methods can complicate fabrication or alter other key properties.
The rapid proliferation of novel 2D materials21 and the development of extensive databases22,23 have made individual analysis impractical. Consequently, research has shifted toward designing high-throughput workflows to efficiently model 2D magnetic materials and identify candidates with promising Tc values.14–19,24,25
Density Functional Theory (DFT)26,27 remains a cornerstone for modeling 2D magnetic materials, allowing the determination of magnetic ground states.28,29 However, the reduced dimensionality of 2D systems poses unique challenges, as critical magnitudes are often weak, ranging from meV to μeV. This energy scale necessitates careful parameter tuning and introduces difficulties in capturing the coupled nature of electronic structure and spin order.30
Beyond determining magnetic ground states, DFT enables the extraction of exchange parameters through total energy calculations30 or related methods.31 These parameters can be integrated into complementary approaches like spin-wave theory,32 Metropolis Monte Carlo simulations,33 or Green's function methods34 to estimate critical temperatures.
In this review, we aim to equip readers with the theoretical and practical tools necessary for modeling 2D magnets using DFT. We begin by introducing basic principles of magnetism, followed by an overview of DFT and its associated challenges, such as electron delocalization and spin-order determination. Subsequently, we describe techniques for deriving parameters for magnetic Hamiltonians and estimating critical temperatures, providing a comprehensive comparison of methods tailored to different material anisotropies. This review serves as a foundational guide for researchers navigating the complex landscape of 2D magnetism.
As an initial step, we begin by recalling that a generic Hamiltonian for electrons in a crystalline array of atoms can be written, within the Born approximation, considering three main contributions:
Ĥ = Hkin + Hlat + He–e | (1) |
![]() | (2) |
The Coulomb interaction can be further split into three contributions:
• The direct or Hartree term VHij = Vi,j,i,j, which, after some manipulation, can be expressed as the classical repulsive field between electron densities located at sites i and j of the lattice. Due to its classical origin, this term modulates the kinetic energy of the electron and alters the single-particle band structure of the system.
• The exchange term VExij = Vi,j,j,i, which arises due to the Pauli exclusion principle and favors lower energy for electrons with the same spin, owing to the antisymmetry of their orbital wavefunctions. This term is the essential source of ferromagnetism in materials.
• Correlation terms , which represent the dynamical response of an electron's transition from site i to site j to all the remaining electrons in the system. Correlations can substantially modify the ground state and induce phase transitions, such as superconductivity.
Dealing with the full many-body Hamiltonian is a formidable task. As a result, many alternative approaches and approximations have been developed to handle different aspects of the problem. The exchange interaction can be treated exactly within Hartree–Fock mean field theory, while the former plus correlation effects are approximated in different ways within density functional theory.26,40,41
In the subsequent sections we will briefly discuss some of the effective Hamiltonians which are used to model magnetism and its relation to the all electron Hamiltonian, since this will help to understand the proper way of modelling it.
![]() | (3) |
Using the number operator i,σ = ĉ†i,σĉi,σ, we can define the local spin density operators at site i:
![]() | (4) |
These operators represent the z-component of the spin and the spin–flip processes, respectively. Defining 2Ŝxi = Ŝi+ + Ŝi− and 2iŜyi = Ŝi+ − Ŝi− we can rewrite the exchange contribution to the Hamiltonian as:
![]() | (5) |
For systems with tightly bound d- or f-electrons, the kinetic energy is significantly reduced due to the strong localization of the wavefunctions. Under these circumstances, the magnetic interaction is dominated by the exchange interaction between nearest neighbors. Moreover, the quantum spins in the Heisenberg model can be approximated as classical spins for large spin quantum numbers (S ≫ 1) given that quantum fluctuations become negligible. Similarly, at finite temperatures where thermal fluctuations dominate, or in systems with long-range magnetic order, the spin deviations are small enough to justify a classical description. Additionally, in Density Functional Theory codes, the magnetic moment of the atom is usually calculated by integrating the continuous spin density in a region centered around the atom, motivating even more the adoption of a classical approximation in which the spin operator is approximated by a continuous magnetic moment vector. This assignment of a continuous magnetic moment localized around an atom is what makes this kind of model and approximation receive the name of atomistic magnetic models. These approximations simplify the magnetic models and make it a versatile tool for studying magnetic properties in many systems.
Then, in practice, the magnetism in the material is studied with a parametric model, the most simple being the Heisenberg model:36
![]() | (6) |
Depending on the symmetry of the system and the anisotropy of the exchange interaction, the Heisenberg model can be reduced to simpler forms, the most prototypical being:13,37
• Ising model: the Heisenberg model transitions to the Ising model42 when there is strong anisotropy along one spin direction (e.g., z-axis), making the spin components perpendicular to this axis negligible. The resulting Hamiltonian becomes:
![]() | (7) |
• XY model: the Heisenberg model reduces to the XY model when the spin interactions are restricted to the x–y plane, with negligible contributions from the z-component. The Hamiltonian in this case is given by:
![]() | (8) |
With an exchange constant Jij that is independent of the direction x or y. In this case, the magnetism is said to be of the easy-plane type.
This rationale can be extended to more complex interactions, giving birth to a variety of atomistic magnetic models in which the interactions are modeled by means of parameters that can be calculated from first principles. Again, the adjective “atomistic” is used in this context to highlight the fact that this construction assumes the localization of the magnetic moments around the atoms. In addition to the isotropic exchange interaction discussed above, some important effects to keep in mind are those arising from spin–orbit coupling interactions, which introduce anisotropy in the system. Some of these are:13,37
• Single ion anisotropy (SIA):
![]() | (9) |
• Anisotropic exchange:
![]() | (10) |
![]() | (11) |
These are not the only interactions frequently used in atomistic magnetic Hamiltonians. Some other less frequent but still worth mentioning are:46
• Biquadratic exchange:
![]() | (12) |
Which differs from the isotropic exchange in the square dependence on the dot product of the spins. Therefore, when Bij > 0, the biquadratic exchange favours a collinear alignment independently on whether it is ferromagnetic or antiferromagnetic. However, its dependence on ij can help to lift the degeneracy of states that would be degenerate with a Heisenberg Hamiltonian. Biquadratic exchange is usually expected to be less important than its bilinear counterpart (isotropic and anisotropic exchange). Nevertheless the work of Kartsev et al.47–49 showed in detail how it can have a value of an important fraction of the bilinear exchange. Additionally, magnetic properties such as the Curie temperature were calculated for different 2D magnets with and without biquadratic exchange; and the inclusion of this interaction was found to lead to the best agreement with experimental measurements. The inclusion of the biquadratic exchange also showed a big impact on the shape of the spin-wave spectra.
• Four-spin interaction:
![]() | (13) |
Similarly to the Biquadratic exchange, the four-spin interaction serves to break the degeneracy on the Heisenberg Hamiltonian.
• Dipole–dipole interactions:
![]() | (14) |
As previously mentioned the Kinetic and lattice contribution can be combined into a single particle non-interactive Hamiltonian
![]() | (15) |
![]() | (16) |
![]() | (17) |
Here, 〈i,σ〉 represents the average occupation for spin σ. The difference in spin populations defines the magnetization:
M = 〈![]() ![]() | (18) |
The competition between the kinetic energy, which favors equal spin populations, and the exchange interaction, which lowers the energy for unequal populations, leads to the Stoner criterion for ferromagnetic instability:
UD(EF) > 1, | (19) |
Stoner magnetism arises because the exchange interaction reduces the energy for electrons with parallel spins, lowering the total energy of the system when the spin populations are unequal. This mechanism is distinct from the localized spins in the Heisenberg model, as it depends on the delocalized nature of the electronic wavefunctions.
The Stoner model provides a simple and intuitive picture of itinerant magnetism, particularly in metallic systems such as ferromagnetic transition metals (e.g., Fe, Co, Ni). However, it neglects electron correlation effects beyond the mean-field approximation, which can significantly influence magnetic properties, particularly in strongly correlated systems. Extensions to the Stoner model, such as dynamical mean-field theory (DMFT), address these limitations and provide a more complete description of itinerant magnetism.
• Ferromagnetic order. Is characterized by the parallel alignment of magnetic moments, resulting in a net macroscopic magnetization. The key condition for ferromagnetic order is that the exchange interaction favors parallel spin alignment. In localized electron systems, ferromagnetism arises due to direct exchange or superexchange mechanisms while for itinerant systems, ferromagnetism is driven by the Stoner criterion.
• Antiferromagnetism order. Is defined by an alternating spin alignment that results in no net magnetization. This order arises when the exchange interaction favors antiparallel spin alignment (Jij < 0). In localized systems, antiferromagnetism is often stabilized by superexchange interactions. For example, in Mott insulators, virtual hopping processes between neighboring sites lower the energy when spins are antiparallel. In itinerant systems, antiferromagnetism can emerge due to Fermi surface nesting, where certain wavevectors q connect regions of the Fermi surface. This enhances the susceptibility at q, leading to spin-density waves (SDWs) with periodic modulation of spin density. In some cases, Jij < 0 produces an antiparallel alignment of the spins but with a nonzero net magnetization because one of the atoms of the two sublaticces has a greater magnetic dipole moment. These cases lead to the so-called ferrimagnetic order.53
• Non-collinear magnetic orders. Non-collinear systems are arrangements of magnetic moments that are not oriented in the same direction. The most fundamental case of non-collinear magnets are the in-plane systems, where the spins are contained in an easy plane, with a hard perpendicular axis. Non-collinear magnetism can also present exotic configurations giving rise to phenomena such as skyrmions or spin spirals that arise when the spins form angles with respect to each other rather than aligning parallel or antiparallel. These exotic orders tend to be stabilized by other kind of effects, that compete with the isotropic exchange interaction such as the DMI induced by spin–orbit coupling, magnetic frustration due to competing exchange interactions or the shape anisotropy originated by the magnetic dipoles.
A special case among the possible magnetic configurations is the case of spin-spirals. These states have spins that rotate with respect to the initial alignment a certain angle over a specific direction given by the spin-spiral wavevector q. Spin spirals have a wavelength separating two sites with the same spin direction (phase). A spin spiral showing these features is sketched in Fig. 1. Spin spirals are special because incommensurate spin-spirals such as the ones in 3D γ-iron54–56 or long wavelength spin-spirals cannot be handled in any exploratory supercell approach. Additionally, collinear magnetic order can be regarded as a particular case of spin spiral order when the wavevector lies on the Γ point of the first Brillouin zone or at some points of its high symmetry path.29
![]() | ||
Fig. 1 Sketch of a spin-spiral. The precession axis is taken to be the z-axis. The cone angle θ is specified as the angle between the spin direction and the precession axis. The spin-spiral moves along the direction of the q vector and the difference in phase between two consecutive spins is Δφ = q·R where R is the vector connecting the two sites. Figure extracted from ref. 57. Reprinted figure with permission from S. Mankovsky, G. H. Fecher and H. Ebert, Phys. Rev. B, 2011, 83, 144401. Copyright 2025 by the American Physical Society. |
What we know as of today is that anisotropy can stabilize magnetic order at nonvanishing temperature as shown for CrI3 experimentally1 in 2017. This material presents ferromagnetic behaviour up to 45 K that is stabilized owing to the presence of exchange anisotropy in the z-direction, causing an out of plane spin orientation that produces ferromagnetism.6 Monolayer Fe3GeTe2 presents a similar feature with a ferromagnetic phase driven by strong out of plane anisotropy up to 130 K.3 Hence, in 2D magnets, the isotropic exchange usually dominates the parallel/antiparallel alignment of the spins but the long range magnetic order is stabilized by a source of anisotropy.
Another special feature of 2D magnets is that the diverse composition of 2D van der Waals materials tends to present different species that do not directly participate in the exchange. The presence of these ligands foments indirect mechanisms that allow interactions between neighbouring metallic centers that are too far to interact directly. The overlapping between metal–ligand–metal orbital connections creates new channels of interactions that are called indirect or super-exchange interactions.37 Moreover, more complex overlapping of orbitals can be present in these systems originating super–super or even super–super–super exchange.
In the end, the special keys about 2D magnetism are the important role anisotropies play and the variety of possible exchange paths. Consequently, the core of its ab initio studies focuses on the identification of the many possible mechanisms that produce that anisotropy/exchange and the calculation of its strength. When doing so in the framework of Density Functional Theory (DFT), some difficulties appear due to two main factors:
• The intrinsic problems within DFT† when it comes to the accurate modelling of certain interactions such as exchange, correlations or long range van der Waals interactions.
• The reduced dimensionality of 2D materials makes the energy scales of the anisotropies very low: from meV to μeV. For example, the isotropic exchange constant, which is usually the greatest in magnitude, is usually on the order of the meV where as in 3D materials can be of the order of 100 meV. The calculation of these parameters is strongly influenced by the structure and Hubbard parameters and can also be importantly influenced by more fundamental computational details, such as the pseudopotentials or the approximation used to describe the exchange correlation functional, which is especially important in the case of itinerant magnets.
In the following section, we give a brief introduction to Density Functional Theory and provide a compilation of the necessary tools to simulate 2D magnetic materials within DFT.
The first attempts to establish an alternative formulation to the Schrödinger equation based on the use of total electronic density n(r) and the density functional were introduced in 1927–1930 by Thomas, Fermi and Dirac.61–63 This formulation laid the origins of DFT, with the so-called Thomas–Fermi (or Thomas–Fermi–Dirac) model that introduced the LDA approximation to describe the kinetic energy of electrons. The TFD model paved the way for significant advances in the history of electronic structure, although this nascent model was not able to provide the required quantitative accuracy (Fig. 2) and did not provide a formal and complete theory of electronic structure.
![]() | ||
Fig. 2 Comparison of the electron density of Argon using the LDA-Thomas–Fermi model and Hartree–Fock. The TF approach results in an overall good description of the electron distribution, but importantly fails in the description of the peak structure. Figure extracted from ref. 64. Reprinted figure with permission from W. Yang, Physical Review A, 1986, 34, 4575. Copyright 2025 by the American Physical Society. |
It was in the 1960s, when the ideas behind an alternative formulation of the electronic problems, based on the total charge density of the system, were formally written, giving birth to DFT in 1964–1965 in the hands of Hohenberg, Kohn and Sham.26,59,65 The solid principles of DFT are in the present expressed in terms of the Hohenberg–Kohn theorems (HK). The first HK theorem proves that for any system of electrons interacting via Coulomb interactions in an external potential, the external potential is fully determined (up to a constant shift) by the ground state electron density n0(r). The second theorem states that a functional of the energy E[n] exists for any external potential, that the ground state energy is its global minimum and the electron density that minimizes it is n0(r). Unfortunately, the theorem does not give any hint about how to obtain such functional E[n], but we can express the total energy without loss of generality as:27
![]() | (20) |
The important consequence of these two theorems comes after the realization that the wavefunction of the system is determined once the external potential is known (since the form of the Coulomb repulsion is known already). But the first theorem states that the external potential is determined by n0(r). As a result, knowing the ground-state electronic density n0(r) is equivalent to knowing the wavefunction of the system. This establishes the electronic density as the fundamental variable of the system, determining all its properties. And this result is extremely convenient since the electronic density is a function of three variables where as the wavefunction depends on 3N variables.
The next key result came in 1965 by Kohn and Sham.26 They proposed to assume that there is a system of non-interacting particles, called auxiliary system, with the same electronic ground-state density as the real system. Then, to obtain n0(r), we can work with the auxiliary system and the many-body effects can be incorporated via an external exchange–correlation potential expressed as a functional of the density. This way, the total energy of the auxiliary system can be expressed by rewritting (20) as:
![]() | (21) |
![]() | (22) |
Is the classical interaction energy of the charge density with itself. Since (21) is just (20) rewritten, by making them equal and solving for Exc[n]:
Exc[n] = T[n] − Ts[n] + Eint[n] − EHartree[n] | (23) |
Expression (23) shows explicitly how the exchange–correlation functional aims to capture all the many-body effects of the real system in an external potential so that the auxiliary system can be of non-interacting particles. Unfortunately, there is no way to know the exact form of the exchange–correlation functional and DFT ends up being exact in theory but approximate in practice. Most approximations are based on the exchange–correlation energy density of the uniform electron gas26,27,66,67 and then, the success of DFT comes from how extremely simple approximations to the exchange–correlation functional give very good results for many systems making DFT extremely popular.68
The practical development of the Kohn–Sham approach leads to a set of equations:
![]() | (24) |
That must be solved self-consistently in the given order since the Hartree term and the exchange–correlation term depend on the electronic density. The eigenstates ψσi are the so called Kohn–Sham states. These states are in principle nothing else but the eigenstates of the auxiliary system but the success of DFT has made standard the description of the systems in terms of these single particle states. Nevertheless, we remark that these states do not necessarily have any connection with reality, they simply provide a useful and simple language via single particle states that serve to describe a complex system. In practice, the Kohn–Sham states are expanded in a basis set either made by plane waves, localized atomic orbitals or a smart combination of both.
To achieve an accurate electronic structure using practical implementations of DFT, several considerations must be taken into account. These include the choice of an appropriate exchange–correlation functional, the selection of pseudopotentials, adequate convergence parameters, the inclusion of noncollinear spin configurations69 and spin–orbit coupling,70,71 as well as smearing techniques.72–74 Incorporating van der Waals (vdW) interactions can be essential for modeling 2D materials, van der Waals heterostructures, and layered systems, as these weak forces strongly influence structural stability and electronic properties.75–77 They are crucial for capturing interlayer coupling and stacking-dependent behaviors, making them indispensable for such simulations.
While most of these settings are standard features of any DFT tutorial78 and will not be elaborated here, we consider it critically important to address the delocalization problem inherent to all approximations of the exchange–correlation functional79–82 and discuss potential improvements due to their inherent relation to magnetic materials as we will discuss below.
The delocalization error is a well-documented limitation of standard functionals and becomes particularly problematic in systems with strongly localized d and f electrons.83–85 These electrons play a key role in the properties of magnetic materials, as they directly influence exchange interaction, the fundamental mechanism driving magnetic ordering. Furthermore, the Coulomb interaction, which underpins these exchange effects, also determines the degree of electron localization, closely linking it to magnetic phenomena.86,87
![]() | (25) |
This aspect results crucially important to consider in the modelling of materials that present strongly localized electrons, such as the ones present in d or f orbitals. Given these kind of orbitals tend to present natural open shells with unpaired electrons, they are directly behind the origin of magnetism and thus, a correct description of the electron localization is fundamental for the computation of magnetic properties of 2D materials.
From the different methodologies available, Hubbard corrections are among the most significant methodologies to improve the approximation of the electron behaviour in the density functional and have demonstrated to provide an excellent compromise between improvement and computational efficiency,86 being implemented in most of the DFT packages (see Table 1).
DFT code | License | Basis set | sc-DFT | GBT | vdW corrections | DFT+U | U obtention | DFT+U+V | V obtention |
---|---|---|---|---|---|---|---|---|---|
VASP150,151 | Commercial | Plane-waves (pseudopotentials) | Yes166 | Yes | DFT-D2167 | Yes | Yes100 | No | No |
DFT-D3168,169 | |||||||||
DFT-D4170 | |||||||||
Method of libMBD171 | |||||||||
Tkatchenko–Scheffler172–174 | |||||||||
MBD@rSC175,176 | |||||||||
MBD@rSC/FI177,178 | |||||||||
SIESTA179 | GPL (academic) | Localized atomic orbitals (pseudopotentials) | Yes134 (but not public) | Yes (but not tested) | DRSLL180 | Yes71 | No | No | No |
LMKLL181 | |||||||||
KBM182 | |||||||||
C09183 | |||||||||
BH184 | |||||||||
VV185 | |||||||||
FLEUR152,153 | GPL (academic) | FP-LAPW (all electron) | Yes | Yes186 | DFT-D3168,169 | Yes187 | No | No | No |
DRSLL180 | |||||||||
Quantum ESPRESSO114,118 | GPL (academic) | Plane-waves (pseudopotentials) | Yes | No | DFT-D2167,188 | Yes | Yes100,105 | Yes87 | Yes100,105 |
DFT-D3168 | |||||||||
Tkatchenko–Scheffler172 | |||||||||
MBD176 | |||||||||
XDM189,190 | |||||||||
OPENMX154,155 | GPL (academic) | Localized atomic orbitals (pseudopotentials) | Yes | Yes157,158 | DFT-D2167 | Yes191,192 | No | No | No |
DFT-D3168,169 | |||||||||
GPAW146,147 | GPL (academic) | Plane-waves (pseudopotentials) | Yes | Yes | DRSLL180,193 | Yes | No | No | No |
LMKLL181 | |||||||||
BH184 | |||||||||
KBM182 | |||||||||
C09183 | |||||||||
ABINIT194 | GPL (academic) | Plane-waves (pseudopotentials) | Yes195 | No | Yes196 | Yes197 | Yes100,108 | No | No |
DRSLL180 | |||||||||
DFT-D2167 | |||||||||
DFT-D3168 | |||||||||
CASTEP198 | Commercial | Plane-waves (pseudopotentials) | Yes | Yes | Tkatchenko–Scheffler172 | Yes | Yes | No | No |
MBD@rSC175,176 | |||||||||
DFT-D2167 | |||||||||
DFT-D3168 | |||||||||
OBS199 | |||||||||
XDM200 | |||||||||
ELK201 | GPL (academic) | FP-LAPW (all electron) | Yes | Yes | No | Yes | No | No | No |
WIEN2K202 | Academic (paid) | FP-LAPW (all electron) | No | No | DRSLL180 | Yes | No | No | No |
LMKLL181 | |||||||||
SCAN+rVV10203 |
Inspired by the role of the Hubbard model in the description of strongly correlated electrons, DFT+U initially aspired to improve the description of the electrons present in strongly localized orbitals. However, the active implementation of the Hubbard U parameter in the exploration of new materials, has revealed to the community that DFT+U can have a broader impact, by restoring the piecewise-linear relationship between the total energy and the electronic occupations, thereby helping to mitigate the errors arising from the violation of this constraint.
In this sense, DFT+U does not directly address localization and strong correlations. Instead, it helps to solve a more fundamental issue: the piecewise-linear constraint of the exact exchange–correlation functional. The violation of this constrain heavily affects electron localization which is especially relevant in these cases and more importantly, in 2D magnetism in which electron localization and magnetism are extremely correlated.
The Hubbard corrective term and the double counting term can take different expressions depending on the formulation,86,96,97 but we want to highlight how one of the most simple formulations resembles the Hubbard model and the energy penalty for partial occupation of the orbitals within the same l-shell:
![]() | (26) |
We can see a similarity between the product and the Hubbard term
of the Hubbard model. Here
is the sum of all the occupation numbers of the localized states of atom I. The second and third terms of the right hand side of the equality correspond to the energy correction and the double-counting correction respectively. This equation shows that the new exchange–correlation functional has a set of parameters, the set of {UI}, one for each atom with localized orbital states. These parameters are called the Hubbard parameters or simply the U parameters. They are nonnegative and they represent the mean coulomb repulsion energy between two electrons in the same atom and belonging to the same l-shell. As one can see from the second term in the right hand side of (26), the correction term introduces an energy penalty if the occupation numbers are not 0 or 1, favouring localization. The energy penalty has a value of UI independently of the value of the quantum number m. This assumption is justified when the localized orbitals retain atomic-like symmetry i.e. spherical symmetry.86 Therefore, the DFT+U approach has a worse performance when there are physical features removing the approximate equivalency of the orbitals with the same l such as strong SOC or crystal field.86
DFT+U schemes are parametrized by a U input parameter for every single l-shell to which the correction is applied, removing the ab initio character of the simulation. This introduces a problematic situation from a first-principles perspective. Over several years of intense research, different alternatives to deal with this parameter have been explored by the community. Its consistent obtention is indeed of extreme relevance since the value of U can have great influence on the observed results.25,98,99
One of the most extended approaches involves empirically extracting the Hubbard U through the fit to specific properties computable within the DFT framework. Usually one of the most representative target parameters that can be found in literature, is the gap of an electronic band structure. Despite the correction introduced by DFT+U has proved to play an important role in the improvement of the electronic structure description,86,100 DFT was never designed to predict spectroscopic properties. The solely use of a band gap to benchmark the correct DFT+U modelling of a material is a disregarded approach that in addition, strongly limits the predictive capabilities of the first principles approach. As it will be discussed below, there are many other properties that result more reliable and that can considered in addition to the band gap to ensure a correct modelling of the materials under the framework of DFT+U.
A more complete approach consists of performing an extended exploration of the Hubbard U, broadening the window to the computation of different important properties in addition to the electronic structure, such as magnetic moments, exchange parameters, energies and shape of magnon dispersions to name a few. A very robust addition to this approach, is to scan also other external conditions such as, the geometrical distortions or the charge modifications induced by strain and electrostatic doping simulations, respectively. The computation of different sophisticated properties and their evolution under these external stimuli, easily accessible from the DFT calculations, provides a very valuable dataset of information. The critical analysis of this information is an effective technique for assessing the DFT+U model and determining the Hubbard parameter from a broad perspective that accounts for both the model's performance and limitations.
Fortunately, the U parameter represents the average Coulomb repulsion energy between electrons within the same l-shell and its physical meaning can be exploited in order to design ways to estimate it. Hence, a third approach to determine it is simply calculating it with a systematic method. The most famous methods are the Hartree–Fock (HF) approach to calculate it,101,102 the linear response approach,86,100,103–105 the constrained random phase approximation106–108 (cRPA) and some machine learning approaches109,110 that include Bayesian optimization.111,112 All these methods propose a fully systematic manner to calculate U, making the DFT+U flavour a zero free parameters approach and going back to the first principles philosophy.
The HF approach basically calculates the average Coulomb repulsion energy from unrestricted HF calculations. The linear response approach aims to tune the U parameter so that the total energy as a function of the number of electrons gets to be linear-piecewise, following the constraint of the exact functional.88 This linear response approach has the advantage of offering a fully self-consistent way of calculating the parameters within DFT. Moreover, it has been reformulated recently in Density Functional Perturbation Theory (DFPT),103,104,113 offering better performance and accuracy and even a fully automated package to compute it called HP105 available as part of the DFT code QUANTUMESPRESSO.114 This process is shown in Fig. 3. The structure is initially relaxed for an initial value of U = Uin, which can be zero. Then, the ground state is obtained and a value of U = Uout is calculated for that ground state. Afterwards, a new ground state is calculated with the obtained Uin = Uout. For this new ground state, a new value of Uout is calculated. The process is then repeated until the difference between the input U and the output U is less than a threshold.
The DFT+U approach serves to describe more localized regimes. Following a similar idea, DFT+U+V was born:87 a different approach capable of describing more general localization regimes by introducing an energy term that favours partial occupation of orbitals of different atoms that are hybridized. This additional parameter adds more flexibility to force the linear-piecewise constraint by allowing a more general case of localization to be corrected and improved. The DFT+U+V correction without the double counting correction is:86,87
![]() | (27) |
The DFT+U+V approach was used for NiO, Si and GaAs in its original publication.87 For NiO, it slightly improved the resulting bandgap and it greatly improved the description of the DOS. In the case of Si and GaAs, it improved the resulting values of different magnitudes, including the bandgap, with respect to DFT+U. The authors conclude that this approach performs better in Si and GaAs because they are more isotropic. Recall that the DFT+U+V assumes orbital independent electron–electron interactions and by construction it will perform better in highly isotropic systems.
The DFT+U+V approach is widely used110,115–117 and the obtention of V is fully automated thanks to the code HP105 which is part of QUANTUM ESPRESSO114,118,119 and its integration into the automation infrastructure AiiDA.120–122
Despite having ways to computes the Hubbard parameters, the approach has intrinsic limitations that could make it fail when describing some systems. The introduction of V parameter clearly illustrates how the approach can be generalized to consider more general regimes of localization, leaving then room for further improvements. Another possibility involves the implementation of the orbital-resolved Hubbard123 U and V to more correctly describe the Hubbard manifold, thereby circumventing the limitations of the current shell-averaged approximation.
As a general conclusion, the self-consistent computation of the Hubbard parameters is a promising way to go, which recovers the predictive capabilities of DFT and results in an overall good modelling of the materials. However, a self-consistent Hubbard U does not always guarantee a good description of the properties of the system124 and thus, performing a complete analysis of various properties and external conditions provides a complementary methodology to determine the computational details of the DFT+U simulations.
Magnetic properties are heavily influenced by correlation effects47 and the overlap between electronic orbitals. The adequate description of these correlation effects can only be done with DFT+U methodologies and a good estimation of the Hubbard parameter. Additionally the U value affects electron localization. For all these reasons, there is a strong dependence of magnetic properties such as magnetic anisotropy energy (MAE) or exchange constants on the value of U.25,30,47,124–126 These magnetic properties determine the critical temperature of the material (as explained in section 4) and therefore, following this logic, the U parameter affects the final calculated critical temperature as well.
Ehybrxc = a0Eexactx + (1 − a0)EDFTx + EDFTc, 0 ≤ a0 ≤ 1. | (28) |
![]() | (29) |
Hybrid functionals can alleviate some of the deficiencies of standard DFT with the caveat of a higher computational cost due to the computation of (29) and the need of tuning the a0 parameter. Hence, this approach removes the ab initio character of the simulation. This parameter is usually fitted semiempirically and its value is usually set around 0.2–0.3.86
It is important to remark that a higher amount of exact exchange (increasing the a0 parameter) does not necessarily imply a better result. This is because the exchange contribution from DFT, EDFTx has also a contribution for the correlation, the so called nondynamic correlation (more details in ref. 41). Therefore, eliminating part of the correlation by setting a0 = 1 can be critical.
The mixing scheme of hybrid functionals can be extended even further by letting the mixing coefficient be spatially dependent:
![]() | (30) |
![]() | (31) |
Hybrid functionals can help describing magnetic materials in a similar manner as DFT+U: improving the description of the electronic density by improving the description of the localization of the electrons. Hybrid functionals have shown to improve the results obtained with standard DFT functionals.41,127–129 Specially for systems with strongly correlated electrons in which the usual functionals fail. Unfortunately, they have the caveat that its usage introduces an additional input parameter in charge of balancing how much exchange is obtained from HF or from the usual functional. Additionally, the calculation of the exact exchange makes the usage way more expensive than with a standard functional.
In order to compensate for this inaccuracy, several dispersion corrections for DFT exchange–correlation functionals have been designed in the last 20 years. Some of them are included in Table 1. These corrections calculate the dispersion correction in a semiclassical manner. In practice, the total energy computed within the self-consistent process is’:131
EDFT+vdW = EDFT + Edisp | (32) |
![]() | (33) |
The details of the calculations of the terms Cn,AB are beyond the scope of this review, but we invite the interested reader to start with published reviews in the topic.130,131
In the particular case of 2D magnetic materials can be stacked to form bilayers or multilayers that are bonded by van der Waals forces. Additionally, the reduced dimensionality in monolayers makes subtle effects more important since they are not as easily quenched by other interactions as it happens in 3D materials. Therefore, when working with 2D magnetic materials, specially when working with more than one layer, it is a good practice to include any of these corrections to properly account for the long range correlations.
In simple terms, finding the magnetic ground state is, in all cases, finding the spin arrangement on the lattice that minimizes the total energy of the system. The most widely used approach to find the magnetic ground state is the exploration of different magnetic configurations of the material, being the most frequently simulated magnetic configurations the FM and Néel AFM. However, a good practice is to extend this analysis to more complex arrangements of spins, such as the zigzag and stripy AFM (Fig. 4), that usually require the simulation of supercells. On the other hand, bulk materials demand the exploration of the FM/AFM coupling between adjacent layers. In addition to this conventional approach, the community has also explored alternative methods like crystal orbital Hamilton population (COHP) analysis to investigate the role of magnetic ordering in the structural stabilization of quasi-two-dimensional transition metal compounds.132
![]() | ||
Fig. 4 Sketch of the principal magnetic configurations considered in the exploration of the magnetic ground state in honeycomb materials. Black (white) spheres represent spin up (down). The illustrated configurations are: (a) ferromagnetic, (b) Néel antiferromagnetic, (c) zigzag antiferromagnetic, and (d) stripy antiferromagnetic. Figure extracted from ref. 133. Reprinted figure with permission from B. L. Chittari, Y. Park, D. Lee, M. Han, A. H. MacDonald, E. Hwang and J. Jung, Phys. Rev. B, 2016, 94, 184428. Copyright 2025 by the American Physical Society. |
Some spin configurations such as spin-spirals or spin-canted might be difficult to simulate with a supercell approach in standard DFT due to the potential possible change in both the magnitude and the direction of the spins during the self-consistency algorithm. In those cases in which we need an special effort to maintain the desired spin arrangement, we can rely on spin-constrained DFT134 (sc-DFT) which allows to constrain the spins of the system, both in direction and magnitude. We highlight that one rule of thumb when initializing a spin arrangement is to always introduce as an input a value for the spin of the atom slightly higher to the one we expect. This is because DFT implementations tend to decrease this value during the process while the cases in which it increases are less common. So in summary, finding the ground state is about exploring all possible configurations to determine the one minimizing the energy.
The problem of this approach of exploring the landscape of possible spin arrangements is that the number of possible configurations is computationally unreachable,7,13,30,135 even when considering collinear configurations since the magnetic ground state can show a magnetic order with a periodicity that goes beyond the primitive unit cell, requiring the usage of supercells that can dramatically raise the cost of the calculation. This difficult downside can be partially solved by establishing a smart criteria to decide in advance which magnetic configurations are more likely to be the ground state while skipping the calculations of those configurations that are very unlikely to be the ground state. This idea resembles the idea of Bayesian optimization111,112 of balancing exploration and exploitation, being in this case exploitation the idea of skipping some configurations while progressing with those that seem to be close to the actual magnetic ground state.
This idea is put into practice in published workflows such as ref. 136 and 137. In ref. 136, a genetic evolution algorithm is designed in such a way only those configurations with low energy survive. The next generation of configurations is obtained from their ancestors, in such a way the new magnetic configurations to try will inherit partially the order of the parents. In principle, this would maintain the “likelihood” of a configuration to be the actual ground state. This workflow is named Magnene and it is designed to work for both collinear and noncollinear spin configurations.
In ref. 137, the workflow is designed only for collinear spin configurations. They set a ranking of most common experimentally found magnetic ground states and they set the likelihood of a new magnetic configuration based on this ranking. This way, the workflow starts calculating the most probable configurations leading to less time consumption most of the times.
With the guidelines given above, the only way to simulate spin-spirals would rely on a supercell approach and potentially using sc-DFT. This would make the calculations extremely expensive and potentially unreachable for long wavelength spin-spirals. Fortunately, there is an approach capable of considering these cases within a simple unit cell: the Generalized Bloch Theorem (GBT).28,29,54,138–141 This approach has one huge advantage and is the fact that it can model a spin spiral within a primitive unit cell. The GBT takes into account not only the translational symmetry of the lattice but also collective spin rotations over the same axis and along a given direction. The GBT extends the usual Bloch theorem stating that, considering the spin symmetry, the one electron wavefunctions can be written as:140
![]() | (34) |
![]() | (35) |
Using the GBT produces a spin-spiral spectra along the high symmetry path of the BZ. The inspection of this energy spectra gives the wavevector q that minimizes the energy and consequently, with this vector we obtain the corresponding magnetic order that minimizes the energy for the given initial spin alignment.
It is important to remark that the GBT requires an initial spin alignment in the computational cell. Therefore, the GBT should be applied for several different initial spin alignments. Additionally, some collinear configurations can be regarded as limiting cases of spin spirals28 and therefore, the GBT can be applied to the study of collinear magnetic ground states as well.
The GBT+DFT has been already used in previous works such as ref. 28, 29, 144 and 145 where it is applied for several 2D materials and for 3D γ-iron in ref. 56 and 140. The GBT within DFT is already available in codes such as GPAW,146–149 VASP150,151 and FLEUR.152,153 The localized atomic orbitals DFT code OPENMX154–156 has also a GBT implementation.157,158
Unfortunately, exploring that many spin configurations is an unreachable task. That is why it is common to see works that explore just a small number of collinear spin configurations and the ground state is then chosen among them. This approach can easily be implemented in high-throughput workflows such as19,159 that look for 2D ferromagnets from well established databases of materials.
When experimental information about the magnetic order is known beforehand, the reach of the magnetic order with DFT usually turns more into a verification rather than an exploration of all the possibilities. Instead of exploring multiple spin configurations until the energy minimum is found, one usually simply verifies that the experimentally obtained magnetic order is more stable than other similar magnetic configurations.
The realm of 2D materials provides a plethora of different magnetic scenarios, making this characteristic one of its most attractive features. There are many examples of 2D materials with the magnetic orders previously introduced, from the most basic to the most exotic. During a large part of the history of 2D magnetic materials, the scientific community has focused in the exploration of ferromagnetic materials such as the CrX3 (X = I, Cl, Br) or CrSBr monolayers. Other materials are well known by their intralayer-AFM such as the MPS3 (M: Mn, Fe, Co, Ni). Fe, Co and Ni compounds are examples of the zigzag AFM introduced in Fig. 4, by the other side MnPS3 presents an example of a Néel AFM. Also these materials provide a good example of different AFM (FePS3) and FM (CoPS3, NiPS3 and MnPS3) bulk phases. The exploration or verification of the magnetic orders from a simple collinear point of view, is typically applied in 2D materials, and tends to easily find agreement with the experimental findings.133,160,161
Noncollinear magnetic orders add more complexity to the determination of the ground state. The addition of SOC is usually enough to determine the noncollinear behaviour of a material. However, there are important exceptions such as the CrCl3 or CrSBr, where SOC is very weak and the in-plane nature of the spins is a consequence of the exchange and shape anisotropy and thus, the magnetic dipoles.
For example, the monolayer of CrI3 is an out of plane ferromagnet up to 45 K.1 The addition of SOC is able to verify the out of plane easy axis of this material and correctly estimate the difference in energy between the in-plane and out of plane spin alignment.6 In contrast, the ground state magnetic order of monolayer CrSBr is described as an in-plane ferromagnet in the presence of SOC.51,162 The stability hierarchy between the different possible in-plane spin directions agrees with the experimental results163 only when shape anisotropy is considered.160
More complex ground states are represented by the spin-spirals, where the spins are not aligned uniformly, but instead, they rotate continuously in space, originating helical patterns in the spin orientation across the material. Spin spirals are a manifestation of complex magnetic interactions and are often found in 2D materials with competing magnetic exchange interactions or broken inversion symmetry. Particularly well studied cases of spirals in 2D materials are the metal dihalides (MX2) such as the Ni or Co compounds.
In Fig. 6 we show one of the results of ref. 29 for CoI2. The figure shows the energy spectra obtained after doing DFT calculations with the generalized Bloch theorem. This material presents its energy minima between the M and Γ points. Therefore, under the approximations of the GBT, its magnetic ground state is a non-collinear spin-spiral. Ideally, Fig. 6 should be repeated for different initial spin alignments within the computation unit cell. However, that would require even more exploration of configurations and therefore, more time and resources. It is always necessary then to stop the exploration of more configurations based on the quality of the results and commit to what you already have.
![]() | ||
Fig. 6 Spin-spiral spectra of CoI2. Figure extracted from ref. 29 without changes under Creative Commons Attribution 4.0 International License https://creativecommons.org/licenses/by/4.0/. |
Other 2D materials present an strong itinerant magnetism, which complicates their modelling and simulation, such as Fe3GeTe2 and Fe3GaTe2. These materials have attracted the attention of the community, given their very high critical temperature.3,164,165
In conclusion, finding the magnetic order of 2D magnets with DFT involves exploring several different configurations i.e. a lot of time and computational effort. However, in practice, experimental knowledge or simply common sense limit that huge exploration to simply the calculation of a small set of spin arrangements from which the magnetic order will be chosen.
![]() | (36) |
Now, consider a different spin configuration: χj. Then: and since the spin configuration is known beforehand, the last equality can be expressed in terms of the Hamiltonian parameters. With this approach, for a Hamiltonian with n parameters, one needs n energy differences i.e. n + 1 DFT results. This way, we create a system of n equations with n variables that we can solve. Solving the system will give the expression of the magnetic parameters. In practice, the final expressions depend on the type of spin lattice under consideration and therefore they are geometry-dependent. Most of the times, it will be necessary to use a bigger computational unit cell in order to capture more exchange interactions such as nearest or second nearest neighbours depending on the specific material. This is one big caveat since it can increase a lot the computational cost of the method. This method can be applied for both collinear30 and noncollinear configurations.47
The important limitation of this method is that it frequently requires to construct supercells to consider the different magnetic configurations, a task that strongly affects the computational efficiency of this method. In different materials, the computation of supercells imposes limitations in the convergence threshold to converge the charge density, compromising both the convergence effort and the quality of the calculations. Moreover, this method importantly relies on the extraction of n + 1 different energies, and often configurations can reach apparently correct local minima with importantly wrong energies, seriously affecting the extraction of the sensitive meV energies of the exchange interactions. Last but not least, in order to achieve meV–μeV accuracy, the convergence parameters required are usually very large. In conjunction with SOC to capture anisotropy effects, the calculation becomes very expensive. In the end, obtaining exchange parameters by total energies analysis requires both a vast amount of computational resources and a good control over the inputs and parameters so that the DFT calculations can achieve the necessary accuracy.
It is important to remark that this method has the implicit assumption that all spin configurations are eigenstates of the magnetic Hamiltonian . This way, the energy from DFT for a certain spin configuration is mapped to an eigenstate of the same spin symmetry. This puts a significant constraint on the magnetic configurations that can be run in DFT. Hence, most of the times the spins in
are treated as classical vectors so that all spin configurations are regarded as eigenstates of
. This approximation works best when quantum effects are not important i.e. for high values of S and/or high temperatures. This issue has been discussed already in ref. 208. In this work, they highlight that the AFM configuration is not an eigenstate of the Heisenberg Hamiltonian and therefore, it should not be considered for rigorous energy mapping purposes. Instead, one should use the non-interacting magnon state (NIM).32 Considering this state and the FM state for the energy mapping, the resulting expression of the exchange parameter differs from the one considering FM and AFM states:
![]() | (37) |
Another relevant remark in ref. 208 is that the value of β is bigger in 2D than in 3D materials. This hints that the accuracy of the energy mapping would likely be worse in 2D materials.
One relevant remark is that since the GBT does not consider SOC explicitly but perturbatively, the fitting of the spin-spiral spectra must be done only with those interactions that are present without SOC. For those that arise with SOC such as SIA or DMI,37 a similar approach can be done but with the spectra that contains only the perturbative contribution to the total energy.
One advantage of this method is that it allows to verify if all the interactions included in the magnetic Hamiltonian are enough to describe the system. If the resulting fit is not capable of describing some features of the spectra, it is likely that one necessary interaction has not been considered. The main disadvantages are those of using the GBT.
In 1987, LKAG obtained the Heisenberg exchange constants by considering small changes in the total energy of the system due to small perturbations of the spins of the ground state.31 When the perturbation is small enough, the energy variation can be calculated using the so called magnetic force theorem:31
![]() | (38) |
Jx is another code based on the same LKAG and magnetic force theorem approach.221 However, the code calculates only the symmetric exchange. Jx is compatible with any output of Wannier90 and directly compatible with OpenMX.
Last but not least, we have Nojij,222 capable of calculating the symmetric exchange from Siesta calculations. Nojij is being extended to a package named Grogu223 which is capable of calculating symmetric exchange and SIA. Unfortunately, to our best knowledge, Grogu is not yet publicly available.
One of the most important steps when modeling the magnetism of the material is the selection of the atomistic Hamiltonian and the included interactions. All the examples in the paragraph above were capable of giving such a prediction because those interactions were included in their model magnetic Hamiltonian. It is common to see in literature how certain interactions are assumed to vanish. For example, there are many studies in which the two-ion anisotropy is assumed to be zero along the x and y directions. For obvious reasons, this assumption is incompatible with the potential prediction of an in plane ferromagnetism driven by anisotropic exchange. Analogously, the SIA is usually taken to be along the z axis i.e. out of plane direction. However, this assumption can be wrong, and the clear example is monolayer CrSBr which has an in-plane easy axis and in-plane ferromagnetism driven by SIA.
In conclusion, choosing an adequate amount of parameters to include in the model is equally important as their calculation. As the widespread quote says: “you can only be as good as your model”.
In a lattice of identical magnetic ions, the critical temperature with MFT is as follows:226
![]() | (39) |
In two-dimensional materials, the presence of magnetic order necessitates anisotropy. However, eqn (39) does not account for either anisotropy or dimensional effects, rendering it inapplicable for precise critical temperature calculations in low-dimensional systems. Nevertheless, this equation can still serve as a rough estimate or an upper bound for the magnetization. Studies indicate that eqn (39) tends to yield a critical temperature roughly twice that obtained from Metropolis Monte Carlo simulations.5,14 Such simulations generally require an as input parameter the maximum temperature to be considered. For this use case, eqn (39) provides a convenient and systematic method to establish a maximum temperature to be simulated in Monte Carlo simulation or equivalently, an upper bound for the critical temperature.
The Mermin–Wagner theorem4 implies that magnetic anisotropy is required for long-range magnetic order in two-dimensional materials with infinite system size. The Ising model, however, implicitly assumes an infinitely strong single-ion anisotropy, which can lead to an overestimation of the critical temperature in two-dimensional systems, as noted in previous studies.15,205,229 For instance, in ref. 205, the calculated Ising model critical temperature for both bulk and monolayer CrI3 was found to exceed experimental values significantly – more than twice the experimental value for monolayer CrI3 and nearly four times that of the bulk material. Similarly, ref. 230 reported a critical temperature of 161 K for CrI3 using the Ising model, far above the experimentally observed 45 K.1
Despite these discrepancies, the Ising model's critical temperature remains a practical upper bound when the lattice structure is one of the studied cases, and the exchange parameters are known. This upper bound is useful for setting a maximum temperature in simulations of magnetization versus temperature, such as Metropolis Monte Carlo methods, providing a systematic approach for simulation limits.
In spin-wave theory, the traditional spin operators in the atomistic magnetic Hamiltonian are substituted by the so called Holstein–Primakoff (HP) spin operators. These are:231,232
![]() | (40) |
The physical intuition of these operators is that the spin state |S,mS = S 〉 is mapped to the vacuum bosonic state. Then, the states with lower mS are obtained as excitations of the vacuum bosonic state i.e. acting with b†i creates a boson and decreases mS by one. As a consequence of the (2S + 1) possible values of mS, there can be at most 2S bosons when considering a single magnetic ion. The appearance of the square root term has two functions. The first one is to maintain the commutation relations of the spin operators i.e. [S+S−] = 2Sz and [SzS±] = ±S±. The second one, to keep the number of bosons to be 2S at most. Consider the state |S,mS = −S〉. In the bosonic picture, this state corresponds to |2S〉B i.e. 2S bosons. If we want to add one more boson, equivalently, applying the S− operator, we get S−|2S〉B = 0 since the term in the square root vanishes. This way, the number of possible bosonic states is kept to (2S + 1) and the map between the spin states and the bosonic states is biyective. With this approach, the Holstein–Primakoff operators are equal to the traditional spin operators in the sense that they have the same matrix elements.
In the bosonic picture we have described before, the bosons are interpreted as magnons which break the magnetic order.
In practice, these operators are substituted in the desired magnetic Hamiltonian and then, expanded with their first order Taylor term:
![]() | (41) |
Operating this expansion leads to a sum of products of bosonic terms times powers of . Products of four or more bosonic operators which correspond to magnon interactions are usually neglected. This is called linear spin-wave theory because only the linear terms are included. This approximation is reasonable at low temperatures, when these interactions are not dominant6,15,52 or at high values of S, where the multiplying factor
makes these contributions negligible. When interactions are included, the process is referred as nonlinear spin-wave theory. The main problem of nonlinear spin-wave theory is that despite including interactions, usually they have to be treated in a Hartree–Fock-like or mean field manner6,52,232,239 and the final results end up being equally inaccurate at high temperatures as shown in ref. 15 for both the linear and nonlinear case.
Either approach and the usage of Bloch's theorem240 leads to the energy dispersion relation E(k) that can be used to count the total number of magnons. Using linear spin-wave theory on ferromagnets leads to a dispersion relation for low |k| as:
E(k) ≈ Δ0 + ρk2 | (42) |
In this context, Δ0 denotes the spin-wave gap, while ρ represents the spin-wave stiffness. The spin-wave gap Δ0 depends on both single-ion anisotropy and exchange anisotropy.6,15,241 However, the precise functional form of this dependence is determined by the specific magnetic Hamiltonian employed. This dependency has a physically intuitive basis: in the absence of anisotropy, the system exhibits isotropy in spin orientation, implying that all spin directions are energetically equivalent. Consequently, it is possible to rotate all spins in the system without altering the total energy. Under these conditions, a spin may be excited with minimal energy to deviate from its initial state, thereby enabling the generation of magnons at very low energies. In two dimensions, this facilitates the destabilization of magnetic order due to low-energy magnon excitations. This interpretation is consistent with the fact that the critical temperature increases with increasing spin-wave bandgap.6
Nonlinear spin-wave approaches lead to an energy dispersion that depends on the total number of magnons:15,32,225,232
E(k) = E(k;Nmagnons) | (43) |
But since magnons are spin 1 bosons, the obey Bose–Einstein statistics and the total number follows the Bose–Einstein distribution:15,232
![]() | (44) |
![]() | (45) |
Since magnons are spin 1 bosons, we can calculate the total magnetization (assuming a single magnetic ion per unit cell) as:
![]() | (46) |
![]() | (47) |
And in the end, obtaining the magnetization with temperature is about solving (46) or (47) numerically as a function of temperature. When nonlinear spin-wave theory is applied, the solution of (43) and (46) or (47) must be done self-consistently. Afterwards, the critical temperature can be detected by inspecting where the magnetization or the spin-wave gap vanish.
Even though spin-wave theory fails to describe properties at high temperatures, in particular the critical temperature, we include it here for the sake of completeness and because its accuracy at low temperatures can be very useful when aiming to calculate magnitudes that depend on the magnetization and its derivatives.
However, as noted in ref. 237, this truncation also breaks the commutation relations of the spin operators and eliminates the restriction S−|2S〉B = 0 since now S− is just a truncation of what it originally was. Hence, after the truncation, one can get bosonic states with more than 2S bosons, and those states are unphysical.
In an attempt to tackle this issue, in ref. 237 they propose to expand the square root with Newton finite differences236 instead of a Taylor expansion. They show that the square root can be written as:
![]() | (48) |
![]() | (49) |
Recall that in the context of the HP operators, a magnetic ion of spin S can have at most 2S magnons. Then the results in (49) proof that truncating (48) for r = 2S and plugging it into the HP operators in (40) gives a set of spin operators that follow exactly the traditional commutation relations and are equal to the spin operators matrix-element-wise. Therefore, this would be an exact representation.
The presented advantages should lead to an improvement when carrying out spin-wave theory with this finite differences representation. However, to our best knowledge, there is no quantitative comparison with respect to the usual Taylor expansion approach.
There have been other approaches towards the obtention of an exact spin representation up to a certain number of bosons as done in ref. 238 by using a set of differential equations. The interested reader can find more information in ref. 238 and the references therein.
The Metropolis Monte Carlo algorithm33,242,243 (MMC) addresses this challenge by following a systematic criterion to explore various spin configurations of the system. First, assume that all spins behave as classical particles, a premise whose validity will be discussed later. The probability of finding the system in a microstate k with energy εk is given by Boltzmann statistics:
![]() | (50) |
![]() | (51) |
![]() | (52) |
When this ratio is greater than one, it follows that , or equivalently, εk < εj. The Metropolis algorithm thus attempts to find the system's ground state by systematically sweeping through spin configurations and evaluating which microstate is more probable by repeatedly calculating the ratio in (52). The full algorithm proceeds as follows:
1. Initialize the system in a specific microstate, e.g., a ferromagnetic configuration.
2. Select one spin of the system and rotate it.
3. Calculate the ratio (52) α. Choose a random number r in (0,1).
4. If n > r, update the microstate of the system.
5. Repeat steps 2–4 until the system stabilizes.
6. Calculate desired quantities: total energy, magnetization, etc.
Note that step 4 is designed such that some less probable configurations can still be accepted. The optimal value for this acceptance rate is 0.5.33,242 The acceptance rate is controlled by the algorithm that rotates the spins. A naive approach would involve simply performing a 180-degree flip, resulting in an Ising-like simulation. The Ising model implicitly assumes infinite single-ion anisotropy (SIA). Since anisotropy is responsible for long-range magnetic order in two dimensions, this approach would lead to an overestimation of the critical temperature.205 This can be confirmed with the examples shown in Table 3.
The specific algorithm used to rotate the spins significantly impacts the efficiency and accuracy of Monte Carlo simulations. Various methods have been developed to rotate spins continuously, aiming to achieve an optimal acceptance rate of approximately 0.5.247 For a detailed discussion on these approaches, see ref. 243, 247 and 248.
It is essential to note that this method relies on the assumption that spins behave as classical particles. This classical approximation becomes invalid at low temperatures, where quantum effects play a significant role.205,249 Consequently, the Metropolis Monte Carlo approach may yield inaccurate results for magnetization values at temperatures near absolute zero or far from the critical point. This limitation should be considered when calculating temperature- and magnetization-dependent properties.
As temperature increases, thermal fluctuations dominate, effectively reducing the significance of quantum effects. At sufficiently high temperatures, the classical approximation becomes appropriate, allowing the Metropolis algorithm to perform well in estimating critical temperatures and magnetization.14 The primary limitation of this method lies in its computational expense, which can become significant for large or complex systems.
One last important aspect to consider is the fact that MMC takes as an input a finite size system and the Mermin–Wagner theorem applies only to infinite size system. Therefore, the convergence of the MMC result as a function of the system size must always be taken into account to analyze finite size effects. This has been studied in ref. 250, where they showed that finite size effects can even stabilize 2D magnetism in absence of anisotropy.
![]() | (53) |
The results of applying this method are shown in Fig. 7 for 2D CrI3.205 The plot clearly shows that the magnetization given by the MMC underestimates the experimental magnetization up to the critical point, featuring an unphysical linear decrease of the magnetization for low temperatures. However, after the quantum rescaling, the magnetization overlaps with the experimental one even for low temperatures. Same trend is observed for bulk CrI3 in ref. 205 and for bulk metals Co, Fe, Ni, Gd in ref. 249.
![]() | ||
Fig. 7 Comparison of the experimental magnetization and the one obtained from a Metropolis Monte Carlo method with the non-Heisenberg Hamiltonian in ref. 205 before and after the temperature rescaling of eqn (53). The non-Heisenberg Hamiltonian is a Heisenberg Hamiltonian with SIA, exchange anisotropy on the z-axis and biquadratic exchange. The symmetric exchange and the anisotropic exchange were taken up to three nearest-neighbours where as the biquadratic exchange was considered only to nearest neighbours. Figure extracted from ref. 205 © 2020 Wiley-VCH GmbH. |
To summarize, when looking for values of the magnetization at low temperatures, the MMC method will not give accurate results. Fortunately, when some experimental data is available, the quantum rescaling helps to recover the whole magnetization vs. T plot trivially by using (53).
VAMPIRE features an extended manual describing all the available options. The developers have implemented an adaptative spin rotation scheme for the MMC that aims dynamically for a 0.5 acceptance rate.247 It also includes a tensorial interaction input that allows to insert manually the strength of the interactions. Hence, the user can insert manually as many interactions as needed.
SPIRIT stands out for featuring its own GUI and AiiDA120–122 plugin package,254 allowing an immediate insertion in automated workflows (more information in the corresponding section).
FIDIMAG MMC is contrained to considering only symmetric exchange interactions, DMI, cubic anisotropy and an external magnetic field.
For more details about the individual codes, the reader is referred to the corresponding references.
GAB(t,t′) := 〈〈A(t); B(t′)〉〉 := −iθ(t − t′)〈[A(t),B(t′)]〉 | (54) |
![]() | (55) |
Exploiting the properties of these functions and their Fourier transforms, the following equation can be obtained:
![]() | (56) |
Which is an algebraic equation and is therefore easier to solve than a differential equation. We can see how this equation relates one Green function with another one that involves a commutation with the Hamiltonian. We say in this context that the equation is relating a Green function with a higher order Green function.34 This relation could be extended indefinitely for higher order Green functions so in practice, the chain has to be stopped at a certain order so that a system of equations is obtained that yields the Green functions. The chain is stopped by factoring the Green function of a certain order in terms of those of less order. This process is called decoupling.34 The utility of these decouplings in our context lies on how they can isolate the value of the 〈Ŝz〉 which serves to calculate the magnetization with temperature of a ferromagnet thanks to the translation symmetry.
Consider the operators Si− and Si+ and let us denote Gij(ω) = 〈〈Si−;Si+〉〉. We are omitting the hat over the spin operators on purpose for the sake of simplicity.
![]() | (57) |
![]() | (58) |
Before proceeding on, we need to specify a magnetic atomistic Hamiltonian. We will follow the process and notation in ref. 34 for a Heisenberg Hamiltonian with single ion anisotropy along the z-axis:
![]() | (59) |
![]() | (60) |
The next term to compute is:
![]() | (61) |
With this equality we can write (57) as:
![]() | (62) |
In order to simplify further, we need to introduce a decoupling. The decoupling known as Random Phase Approximation (RPA) goes as follows:
![]() | (63) |
And the justification for this approach is purely heuristic.86 Note that implicitly, the RPA is setting Sz = 〈Sz〉, which is an assumption that is valid only at low or far from the critical point temperatures. Using the RPA (63) in (62), the Fourier expansions and the translational invariance of the ferromagnet (〈Szi〉 = 〈Szj〉 = 〈Sz〉), one can arrive to the system of equations:
![]() | (64) |
ωkRPA = B + 〈Sz〉(J0 − Jk) | (65) |
![]() | (66) |
So finally, the Green functions method with the RPA yields a system of eqn (65) and (66) to be solved self-consistently in order to obtain the total magnetization of the ferromagnet with temperature.
The Random Phase approximation can be extended for systems with arbitrary spin, yielding:225,232
![]() | (67) |
Another famous decoupling scheme is the Callen decoupling34,258
〈〈SziSl+;Sj−〉〉 ≃ 〈Szi〉〈〈Sl+;Sj−〉〉 − α〈Si−Sl+〉〈〈Si+;Sj−〉〉 | (68) |
![]() | (69) |
![]() | (70) |
Which has the same form as for the RPA but the Callen decoupling introduces an additionally term in ωk to be considered along the self-consistency process.
So in the end, using the Green's functions approach is about finding an appropriate decoupling scheme in order to arrive to a system of equations to be solved self-consistently to obtain the magnetization of the system. The decoupling schemes we have introduced here do not require using higher-order Green's functions but this is not a constrain and the philosophy can be extended to higher order Green's functions.
One of the main disadvantages of this approach is that every single atomistic magnetic Hamiltonian has to be studied separately, including both the derivation of the equation as well as their implementation into a code. Additionally, the decoupling process can be regarded as opaque in comparison with other methods that have a more clear physical meaning. Last but not least, every single decoupling scheme has a temperature regime in which it makes sense, but the extent on this validity interval is unknown and so far the only choice we have to test the method is via benchmarking and comparison with other methods.14,232
Despite the inconveniences, the Green's functions approach is widely used as well to compute critical temperatures of 2D magnetic materials as shown in recent literature.14,229,232,259
For 2D magnets, using eqn (47) with a dispersion relation as (42) without spin-wave gap results in no magnetic order at nonzero temperature. With nonzero spin-wave bandgap, the dependence can be obtained by solving (47), which gives the dependence for low temperature.
The magnetization in the vicinity of the critical point is predicted by renormalization group theory, which predicts for the Ising and Heisenberg Hamiltonians that the magnetization close to the critical point evolves as where β is the so called critical exponent.260 In fact, for these magnetic Hamiltonians, the magnetization close to the critical temperature is well described by the equation:
![]() | (71) |
It is reasonable to assume that a more general Hamiltonian would follow a similar behaviour to (71) in the vicinity of the critical point. However, the exact value of the exponent is unknown and in fact it will depend on the exact form of the Hamiltonian. However, we can use the value of the critical exponent as a metric of how Ising-like the material is i.e. the closer the critical exponent is to 1/8, the more Ising like it is.
Using the temperature rescaling in (53) in (71) leads to the so called Curie–Bloch equation:249
![]() | (72) |
This new model solves the problem of the difference in shape of the experimental vs. calculated M(T) curves for low temperatures. The magnetic Hamiltonians and the statistical distribution in the MMC are of classical nature and therefore, they do not precisely take into account quantum effects that are relevant at low temperatures. By carrying out this temperature rescaling, the shape of the M(T) curve is improved, replicating the shape of experimental curves.205
![]() | (73) |
![]() | ||
Fig. 8 Critical temperature for a square lattice as a function of rescaled SIA (A/J) calculated for S = 1, S = 3/2, and S = 2 with ferromagnetic exchange coupling using the RPA (In this context, RPA refers to nonlinear spin wave theory). The dashed line represents the critical temperature given by the Ising Hamiltonian, which sets the case for infinite SIA. Figure extracted from ref. 15 D. Torelli and T. Olsen, Calculating critical temperatures for ferromagnetic order in two-dimensional materials, 2D Materials, 2018, 6, 015028. Published 17 December 2018. DOI: 10.1088/2053-1583/aaf06d. © IOP Publishing. Reproduced with permission. All rights reserved. |
The first remarkable result is that nonlinear spin-wave theory dramatically fails for strong SIA ratio A/J. The authors attribute this failure to the appearance of strong magnon–magnon interactions that are not well described by the characteristic mean field treatment in nonlinear spin-wave theory. The failure is so bad that the predicted critical temperature even surpasses the Ising temperature which corresponds to the infinte SIA case.
However, nonlinear spin-wave theory agrees well with MMC at low SIA ratio. In fact, it does a better job for low SIA ratio since it predicts a 0 Curie temperature for no SIA (completely isotropic system in this case) where as MMC predicts a nonzero Tc, contradicting the Mermin–Wagner theorem.
The critical temperatures obtained by MMC are shown in Fig. 9 for three types of 2D lattices. The MMC results agree with the Ising limit (dashed lines) for high SIA ratio. The key result in ref. 15 is the expression of the critical temperature obtained by fitting the results of the MMC in Fig. 9:
![]() | (74) |
![]() | ||
Fig. 9 Critical temperature as a function of scaled SIA (A/J) calculated with classical Monte Carlo simulations for honeycomb, square, and hexagonal lattices with ferromagnetic exchange. The solid lines are obtained from the empirical fitting function (74). The Ising limit is indicated by dashed lines for the three lattices. Figure extracted from ref. 15 D. Torelli and T. Olsen, Calculating critical temperatures for ferromagnetic order in two-dimensional materials, 2D Materials, 2018, 6, 015028. Published 17 December 2018. DOI: 10.1088/2053-1583/aaf06d. © IOP Publishing. Reproduced with permission. All rights reserved. |
When B ≠ 0 and S ≠ 1/2, they found an analogous empirical relation:
![]() | (75) |
The magnetic Hamiltonian they proposed to model the ferromagnets is:
![]() | (76) |
![]() | (77) |
Which features only two-ion anisotropy (no SIA). The following notation has been used: and
is the out-of-plane anisotropic exchange strength. With this notation Jxxij = Jyyij = Jij(1 − Δij) and Jzzij = Jij(1 + Δij).
Tiwari et al. looked for 2D magnets in C2DB with positive first neighbours exchange constant Jij := J and anisotropic exchange strenght Δij := ΔNN. They found 34 2D magnets with J > 0 and ΔNN > 0. For those 34 magnets, they calculated the critical temperature using non-linear spin-wave theory, the Green's functions method and MMC. With the resulting Tc values, they fitted the results to the following model:
![]() | (78) |
Lattice | Parameter | MC | Green | RNSW |
---|---|---|---|---|
Honeycomb | α1 | 0.49 | 0.07 | 0.40 |
α2 | 0.14 | 0.37 | 0.62 | |
Hexagonal | α1 | 0.24 | 0.24 | 0.32 |
α2 | 0.045 | 0.14 | 0.21 | |
Square | α1 | 0.37 | 0.34 | 0.43 |
α2 | 0.08 | 0.24 | 0.36 |
The calculated values for Tc with each method and the resulting fit for a honeycomb lattice are plotted in Fig. 10. The same plot for an hexagonal and a square lattice can be found in the ESI of ref. 14 but they follow the same trends as Fig. 10, hence, we discuss only the results for the honeycomb lattice.
![]() | ||
Fig. 10 Comparison of the calculated critical temperatures with the Green's functions methods, the MMC and non-linear spin-wave theory (RNSW) as a function of the nearest neighbour anisotropic exchange strength ΔNN for a honeycomb lattice. The dashed lines represent the fit using eqn (78). The solid line for the MC shows 25th to 75th percentile of the calculated Curie temperature. The yellow shaded area shows the region where the MC and the Green's function methods have a difference of less than 10%. The horizontal tick show the Ising limit (1.52S2, with S = 3/2) and the quantum mean field (1/3[S(S + 1)NNN], NNN = 3 for honeycomb lattice). Figure extracted from ref. 14 without changes under Creative Commons Attribution 4.0 International License https://creativecommons.org/licenses/by/4.0/. |
Fig. 10 has some noteworthy results. Firstly, among the three methods, MMC is the only one that predicts a nonzero critical temperature in the absence on anisotropy, contradicting the Mermin–Wagner theorem. Secondly, nonlinear spin-wave theory gives always a lower bound of the critical temperature with respect to the other two methods. Additionally, the Green's functions method serves as an upper bound for the MMC for high anisotropic exchange strength. Last but not least, there is an interval of anisotropic exchange strength (the yellow shaded area in Fig. 10) in which the three methods differ in less than 10%. These results, combined with those in Fig. 8 suggest that spin-wave theory is not an appropriate method when treating high SIA or exchange anisotropy.
We started introducing the mean field theory expression for the critical temperature and the critical temperatures for the Ising model in 2D lattices. The first one does not account for dimensionality nor anisotropy effects and the second one assumes infinite single ion anisotropy. These facts make them both usually huge overestimations of the critical temperature of the 2D magnets as shown in Fig. 10. Moreover, these formulas are obtained when considering only isotropic exchange and z-component of the anisotropic exchange tensor respectively. Therefore, they do not take into account other effects such as DMI. Their main advantage is the fact that computing them is trivial and we can use their overestimation as an upper bound for the critical temperature. This upper bound can be used for example as an input for the maximum temperature for a Metropolis Monte Carlo or for a spin-wave theory algorithm.
Spin wave theory in both its linear and nonlinear flavours has the strong disadvantage that its explicit development has to be done for every different atomistic magnetic Hamiltonian. We are then forced to write the equations in advance before writing an automatic code. Moreover, the expressions can get long and complicated when including all the anisotropies we explained in this text. Additionally, the truncation of the Holstein–Primakoff expansion leads to neglecting the magnon–magnon interactions that can have importance at high temperatures. This is the reason why the magnetization obtained with spin-wave theory at high temperatures will not be reliable. Moreover, in ref. 15 and Fig. 8, it is shown how spin wave theory fails to give a sensible value of the critical temperature for high SIA/J situations. The authors attribute this result to magnon–magnon interactions introduced by SIA/J that are not correctly described within NLSW theory. However, for low SIA/J regimes, Fig. 8 and 10 show how spin-wave theory yields results compatible with MMC. In fact, it gives 0 Tc for zero anisotropy, something that MMC does not show. Additionally, there is an anisotropic exchange window shown in Fig. 10 in which MMC and NLSW theory differ no more than 10%. For all these reasons, we consider spin-wave theory suitable only for low SIA/J and exchange anisotropy cases, which in the end will give low Tc values.
For cases with higher SIA/J, MMC is more suitable as concluded from Fig. 8 and 9. We also want to highlight that treating the spins as classical particles is an assumption with validity only when quantum effects become unimportant because they are quenched by temperature. Consequently, we consider MMC the most appropriate method when looking for above room temperature 2D magnets. Its main downsides are its extremely high computational cost and its unphysical low magnetization obtained for below Tc temperatures as explained for Fig. 7. However, this unusual low magnetization is not a problem if we are only interested on obtaining Tc.
The Green's functions method presents similar problems to those of spin-wave theory. Firstly, the equations have to be derived in advance for every atomistic magnetic Hamiltonian, and these get more and more complex when the number of different interactions is increased. Secondly, the decoupling scheme dictates the validity of the approach. For example, the RPA decoupling sets Sz = 〈Sz〉, which is an assumption valid only at low temperatures when fluctuations are low. The Callen decoupling serves as an interpolation between the RPA and a high temperature regime34 and therefore, there is no other way to estimate its accuracy rather than doing benchmarks as done in ref. 14 with Fig. 10. This figure along with Fig. 8 discards spin-wave theory as a good method in cases with high SIA/J or exchange anisotropy and suggest that there is a regime of exchange anisotropy in which MMC and the Green's functions method give compatible results. The comparison of the methods as a function of SIA/J is done in ref. 232 and shown in Fig. 12 for a honeycomb lattice. This figure finally shows that for low SIA/J and zero anisotropic exchange, all methods give comparable results. It is worth noting that NLSW theory remains being a lower bound for the rest of the methods where as the RPA serves as an upper bound.
In conclusion, once the parameters of the chosen Hamiltonian are known, choosing the method to calculate the Curie temperature is not an arbitrary choice. The researcher has to choose the method depending on the value of the parameters, its estimation of the Curie temperature with the mean field approach and the availability of computational resources. Last but not least, we want to clarify that the methods we introduced above in section 5 can be extended to antiferromagnetic materials.
As a summary, in this review we have revisited the fundamentals of magnetic interactions in 2D materials from the approximated description of these systems using DFT to the picture of atomistic magnetic models. We have discussed the main critical points to consider in the simulation of magnetism in this kind of materials, which introduce important challenges driven by the small energies of both exchange and anisotropies and the complexity of the models that involve spin–orbit effects, correlation and van der Waals interactions. We present a guide to address the important challenges related to the modelling of 2D materials, mainly the selection of the Hubbard parameters, descriptions of the van der Waals interactions and approaches to properly model the magnetic ground state and compute exchange parameters, anisotropies and critical temperatures. This review also presents a comparison of the different DFT codes available, highlighting their relevant features for magnetism in DFT. In the end, the whole work compiles the necessary tools to model 2D magnetic materials from DFT to their critical temperature following the steps shown in Fig. 11.
![]() | ||
Fig. 12 Comparison of the critical temperatures obtained for a 2D honeycomb lattice as a function of the SIA/J ratio (A represents the SIA). The anisotropic exchange is set to 0. HP represents nonlinear spin-wave theory with the Holstein–Primakoff expansion. RPA and CD refer to the Green's functions method with the random phase approximation and the Callen decoupling respectively. RPA + CD refers to a combination of both and MC refers to the Metropolis Monte Carlo algorithm. Figure extracted from ref. 232: Varun Rajeev Pavizhakumari, Thorbjørn Skovhus and Thomas Olsen, Beyond the random phase approximation for calculating Curie temperatures in ferromagnets: application to Fe, Ni, Co and monolayer CrI3, Journal of Physics: Condensed Matter, 2025. Accepted manuscript online 6 January 2025. DOI: 10.1088/2053-1583/aaf06d. © IOP Publishing. Reproduced with permission. All rights reserved. |
We finally like to highlight the importance of all the content, tools and guidelines introduced in this section by comparing computed critical temperatures from DFT with the values obtained experimentally. Table 5 contains the critical temperature of several 2D magnets obtained from DFT as well as the DFT flavour used and the method to compute the transition temperature. The critical temperature is a good magnitude to illustrate the complexity of the modeling because it suffers from the propagation of all the errors that have appeared during the modeling process.
Monolayers | Code | DFT flavour (U or standard) | Tc/TN (K) | Tc/TN method | Experimental Tc/TN (K) | J1 (meV) | Ref. |
---|---|---|---|---|---|---|---|
CrI3 | VASP | PBE | 38.42 | Montecarlo fit (75) | 451 | 1.9 | 25 |
VASP | PBE+U (U = 2) | 46.06 | Montecarlo fit (75) | 1.92 | 25 | ||
VASP | PBE+U (U = 1.7) | 42.2 | Montecarlo | 1.06 | 224 | ||
GPAW | PBE | 28 | Montecarlo fit (75) | 0.97 | 16 | ||
GPAW | PBE23 | 41 | Green's Functions | 1.02523 | 14 | ||
GPAW | PBE23 | 31 | Montecarlo | 1.02523 | 14 | ||
GPAW | PBE23 | 23 | NLSW | 1.02523 | 14 | ||
Quantum Espresso | PBE+U (4.55) | 94.2 | Green's Functions | 3.197 | 124 | ||
VSe2-2H | Quantum ATK | PBE+U (U = 2) | 249.69 | Montecarlo | 418263 (multilayer) | 19.52 | 5 |
VASP | PBE+U (U = 2, J = 0.87) | 291 | Montecarlo | 38.8 | 264 | ||
VASP | PBE | 359 | Montecarlo | 18.93 | 265 | ||
VASP | PBE | 427.8 | Exchange-mean-field | 44.85 | 266 | ||
CrSBr | Quantum Espresso | PBE+U (U = 3) | 122 | NLSW | 146162 | 3.54 | 125 |
VASP | PBE+U (U = 3) | 160 | Montecarlo | 3.87 | 267 | ||
VASP | PBE+U (U = 4, J = 0.9) | 175 | Montecarlo | 3.65 | 50 | ||
VASP | PBE+U (U = 3) | 168 | NLSW | 5.68 | 268 | ||
CrBr3 | VASP | PBE | 23.53 | Montecarlo fit (75) | 27,269 34270 | 1.24 | 25 |
VASP | PBE+U (U = 2) | 27.52 | Montecarlo fit (75) | 2.56 | 25 | ||
VASP | PBE+U (U = 1.7) | 23.1 | Montecarlo | 0.67 | 224 | ||
Quantum Espresso | PBE+U (4.20) | 55.3 | Green's Functions | 2.682 | 124 | ||
VTe2-2H | Quantum ATK | PBE+U (U = 2) | 126.61 | Montecarlo | Not done yet to our best knowledge | 9.84 | 5 |
GPAW | PBE23 | 387.6 | Green's Functions | 21.9623 | 14 | ||
GPAW | PBE23 | 371.1 | Montecarlo | 21.9623 | 14 | ||
GPAW | PBE23 | 263.4 | NLSW | 21.9623 | 14 | ||
VASP | PBE+U (U = 2 J = 0.87) | 553 | Montecarlo | 44.3 | 264 | ||
VASP | PBE+U (U = 2) | 301 | Mean field theory | 29.25 | 271 | ||
CrCl3 | VASP | PBE | 13.49 | Montecarlo fit (75) | 12.95,272 17 (2 layers)269 | 0.89 | 25 |
VASP | PBE+U (U = 2) | 16.77 | Montecarlo fit (75) | 1.2 | 25 | ||
VASP | PBE+U (U = 1.7) | 12.1 | Montecarlo | 0.39 | 224 | ||
GPAW | PBE | 9.2 | Montecarlo fit (75) | 0.64 | 16 | ||
Fe3GaTe2 | VASP | LDA | 320 | Montecarlo | 240273 367 (few layers)274 | 57.18 | 165 |
QE | GGA | 320 | Montecarlo | 35.52 | 164 | ||
Fe3GeTe2 | Wien2k | LDA+U (U = 4) | 152 | Montecarlo | 1303 | 11.91 | 275 |
VASP | LDA | 154 | NLSW | 1.21 | 276 |
The normalized deviation of the computed critical temperature from the experimental value is shown in Fig. 13. The plot shows how some approaches have lead to a rather good estimation of the critical temperature (with error of less than 20% or even 10%) while others fail dramatically. Our take home message with this figure is that full control during the whole modeling process is essential to end up getting accurate results and to understand the limitations of the work done.
![]() | ||
Fig. 13 Comparison of calculated critical temperatures of 2D magnets with the experimental value. The data shown comes from Table 5. Since we are not aware of any experimental estimation of the critical temperature of VTe2, we take its value as the mean from those in Table 5. |
We expect this review to motivate future readers, to get familiar with the critical steps during the process of materials simulation using DFT and to be aware of the challenges and limitations of the state of the art.
All authors acknowledge support from PID2019-106684GB-I00 also funded by MCIN/AEI/10.13039/501100011033/FEDER, UE, as well as PID2022-138283NB-I00 funded by MCIU/AEI/10.13039/501100011033 and SGR and DLE funded by Generalitat de Catalunya. ICN2 is funded by the CERCA Programme from Generalitat de Catalunya, and is currently supported by the Severo Ochoa Centres of Excellence programme, Grant CEX2021-001214-S, both funded by MCIN/AEI/10.13039.501100011033.
This work has received funding from the European Union's Horizon Europe research and innovation programme – European Research Council Executive Agency under grant agreement no. 101078370 – AI4SPIN. Disclaimer: Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.
Footnote |
† And mainly those of the current approximate exchange correlation functionals. |
This journal is © The Royal Society of Chemistry 2025 |