Challenges and strategies for first-principles simulations of two-dimensional magnetic phenomena

Abstract

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.

---

Jaime Garrido Aldea

Jaime Garrido Aldea received his bachelor's degrees in Physics and Mathematics from the University of Cantabria, and a master's degree in Quantum Science and Technology from the University of Barcelona. He has held a research internship at the Instituto de Física Fundamental (CSIC), and is currently pursuing a PhD at the Catalan Institute of Nanoscience and Nanotechnology (ICN2) in Spain. He works under the supervision of Jose Hugo Garcia Aguilar and Stephan Roche in the Theoretical and Computational Nanoscience Group. His research focuses on the computational study of materials using Density Functional Theory, with particular emphasis on two-dimensional magnetic materials.

Dorye L. Esteras

Dorye L. Esteras holds a PhD in Nanoscience & Nanotechnology and is a member of the Theoretical and Computational Nanoscience group led by Stephan Roche in the catalan institute of nanoscience (ICN2). With a background in physics and an early research career focused on the first-principles simulation of two-dimensional magnetic materials, he currently works on the development of automated and AI-accelerated workflows dedicated to the discovery of novel quantum materials for emerging technological applications.

Stephan Roche

Stephan Roche is an ICREA Research Professor and leads the Theoretical and Computational Nanoscience Group at ICN2. He has made significant contributions to the understanding of spin, charge, and thermal transport in disordered systems and two-dimensional materials, notably graphene. He pioneered large-scale quantum transport simulations enabling billion-atom modeling. Roche has been involved in the European Graphene Flagship for over a decade, serving as leader of the Spintronics work package and head of Division 1, “Enabling Science & Materials”. He also contributes to international research efforts on advanced materials within the European IAM-I initiative.

Jose H. Garcia

José Hugo Garcia Aguilar is a staff member at the Catalan Institute of Nanoscience and Nanotechnology (ICN2), specializing in quantum transport, spintronics, and 2D materials. With a PhD in Physics from the Federal University of Rio de Janeiro, he has developed scalable numerical methods for simulating quantum effects in disordered systems, some of which are now implemented in widely used simulation tools. His work has led to high-impact publications and collaborations across Europe. He is the principal investigator of the ERC-funded project AI4SPIN, which combines AI and quantum simulations to discover novel spintronic materials.

---

1. Introduction

The first experimental measurement of ferromagnetic order in a monolayer was carried out in 2017 for the out of plane ferromagnet CrI₃,¹ obtaining a Curie temperature (T_c) of 45 K. Shortly after, ferromagnetism in bilayer Cr₂Ge₂Te₆ was found,² with a T_c that can be tuned by introducing small out-of-plane magnetic fields. In 2018, ferromagnetic order was found for Fe₃GeTe₂ up to 130 K for a single monolayer.³

2D magnetic materials have remained so elusive all this time as a consequence of the Mermin–Wagner theorem,⁴ 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,⁹ energy harvesting systems and green energy applications^10,11 and nonvolatile magnetic memories.¹² Moreover, their 2D nature introduces unique opportunities, such as stacking and twisting, to tailor properties for specific applications.¹³

For practical device applications, achieving magnetic stability above room temperature remains a crucial challenge. Efforts are underway to identify 2D magnets with high intrinsic T_c.^5,14–19 While various techniques exist to enhance magnetic stability,²⁰ 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 materials²¹ and the development of extensive databases^22,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 T_c 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.³⁰

Beyond determining magnetic ground states, DFT enables the extraction of exchange parameters through total energy calculations³⁰ or related methods.³¹ These parameters can be integrated into complementary approaches like spin-wave theory,³² Metropolis Monte Carlo simulations,³³ or Green's function methods³⁴ 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.

2. Magnetic interactions in 2D magnets

It would be futile to attempt performing accurate simulations of the magnetic properties of two-dimensional magnets without first grasping their fundamental aspects. Fortunately, there are many books^32,35,36 and comprehensive reviews^13,37 that cover this topic in detail. In this section, we aim to highlight the essential concepts needed to understand the core challenges.

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:

Ĥ = H_kin + H_lat + H_e–e (1)

where the first term corresponds to the kinetic energy of the electrons, the second to their interaction with the atomic lattice, and the last term models the electron–electron interaction arising from the Coulomb field. The first two terms are single-particle energies and can be significantly simplified using Bloch's theorem,^38,39 which describes the system as a collection of non-interacting electrons with an effective mass and constrained momentum defined by the crystal structure. Therefore, the single-particle contribution favors delocalized electrons. The last contribution in the second quantization framework can be written aswhere c^†_i,σ and c_i,σ represent the creation and annihilation operators of an electron in the i-th orbital with a given spin σ = ±1, respectively. These orbitals can be any single-particle basis, but for the purposes of this discussion, it is useful to think of them as localized orbitals at each atom in the crystal. In this sense, the Coulomb field term, characterized by the matrix , contains the Coulomb interaction between pairs of electrons.

The Coulomb interaction can be further split into three contributions:

• The direct or Hartree term V^H_ij = V_i,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 V^Ex_ij = V_i,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.

2.1. Localized electrons and atomistic magnetic models

At the beginning of this section we remind that the operator associated with the exchange interaction can be expressed as:where V^Ex_i,j is the exchange integral between sites i and j, and c^†_i,σ and c_i,σ are the creation and annihilation operators for electrons at site i with spin (σ) = ±1. This form emphasizes the role of the Pauli exclusion principle and the antisymmetry of the wavefunction.

Using the number operator n̂_i,σ = ĉ^†_i,σĉ_i,σ, we can define the local spin density operators at site i:

These operators represent the z-component of the spin and the spin–flip processes, respectively. Defining 2Ŝ^x_i = Ŝ_i⁺ + Ŝ_i⁻ and 2iŜ^y_i = Ŝ_i⁺ − Ŝ_i⁻ we can rewrite the exchange contribution to the Hamiltonian as:

where Ŝ_i = (Ŝ^x_i, Ŝ^y_i, Ŝ^z_i) the spin vector operators at sites i.

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:³⁶

where J_ij is the isotropic exchange constant and S_i and S_j are the total magnetic moments around ions i and j. In practice, S_i and S_j are regarded as unit vectors so that the exchange constant can be expressed in energy units. The isotropic exchange is the dominant interaction in magnetic materials. The energy associated with the isotropic exchange interactions depends solely on the relative orientation of neighbouring spins and the existent overlapping between orbitals that form an interaction path.³⁰ In most cases, it is way higher than the rest of the interactions and it controls the parallel or antiparallel alignment of the spins. With our notation, J_ij > 0 favours parallel alignment of the spins (ferromagnetism) where as J_ij < 0 favours antiparallel alignment (antiferromagnetism).

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 model⁴² 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:

here, 〈i,j〉 stands for a summation that runs over the nearest neighbors of i.

• 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:

With an exchange constant J_ij 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):

where the vector denotes the direction that minimizes the total energy contribution for the magnetic moment , equivalently, the preferential alignment for the spin . is called the easy axis. The parameter K_i is the SIA energy for the spin . Single-ion anisotropy is originated from the interaction between the SOC and the crystal field. It tells us about the spin interacting with the environment and then it is a local property, that involves a magnetic atom and the surrounding coordination sphere formed commonly by the ligands.

• Anisotropic exchange:

where J^ani_ij is a 3 × 3 matrix. The anisotropic exchange tensor can be divided in two different contributions. The first one is the so called two-ion anisotropy,⁴¹ given by the terms k_i on the diagonal. When the two ion-anisotropy goes along bond directions, it is called Kitaev exchange. These terms introduce two-spins exchange anisotropy. The other contribution is given on the antisymmetric terms of the tensor. This contribution is the so-called Dzyaloshinski–Moriya interaction (DMI)^43–45 and it can be rewritten as:where is the so-called DMI vector. The DMI interaction is usually much weaker than the isotropic exchange interaction and therefore, it introduces a small canting of the spins with respect to the direction forced by the isotropic exchange.

These are not the only interactions frequently used in atomistic magnetic Hamiltonians. Some other less frequent but still worth mentioning are:⁴⁶

• Biquadratic exchange:

Which differs from the isotropic exchange in the square dependence on the dot product of the spins. Therefore, when B_ij > 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:

Similarly to the Biquadratic exchange, the four-spin interaction serves to break the degeneracy on the Heisenberg Hamiltonian.

• Dipole–dipole interactions:

where is the vector connecting sites i and j. Magnetic dipole–dipole interactions help to stabilize magnetic orders, specially in 2D vdW systems.^50–52 However, these interactions are usually negligible when calculating critical temperatures since their value is much smaller than the exchange interactions. This energy contribution is usually refered as shape anisotropy and it is SOC independent.

2.2. Itinerant electrons and stoner magnetism

While the Heisenberg and the atomistic models describe localized magnetic moments around the atoms, Stoner magnetism provides a framework for understanding magnetism in itinerant electron systems. In these systems, the magnetic moments arise not from localized spins but from the collective behavior of delocalized electrons in a metallic band structure. The origin of this magnetism lies in the interplay between the kinetic energy of electrons and their Coulomb interaction, which can be expressed within a mean-field framework.

As previously mentioned the Kinetic and lattice contribution can be combined into a single particle non-interactive Hamiltonian

where: t_ij is the hopping integral between sites i and j, describing the kinetic energy of the delocalized electrons. Moreover, in systems where correlations are negligible, one can approximate the Coulomb interaction as an onsite-termwhere U is the on-site interaction strength. Using the mean-field approximation, this term can be linearized as:

Here, 〈n̂_i,σ〉 represents the average occupation for spin σ. The difference in spin populations defines the magnetization:

M = 〈n̂_i,+〉 − 〈n̂_i,−〉. (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(E_F) > 1, (19)

where D(E_F) is the density of states at the Fermi energy. When this condition is satisfied, the system favors a spin-polarized ground state, resulting in spontaneous magnetization.

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.

2.3. Magnetic order and the role of exchange interaction

The exchange interaction is the fundamental mechanism governing the emergence and type of magnetic order in materials. Depending on the nature of the electronic system and the interplay of additional interactions, different forms of magnetic order can arise. This section provides an overview of three main categories of magnetic order: ferromagnetism, antiferromagnetism, and non-collinear magnetism, discussing their origins within localized and itinerant frameworks and their relation with the exchange interaction.

• 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 (J_ij < 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, J_ij < 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.⁵³

• 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 γ-iron^54–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.²⁹

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.

2.4. Special features in 2D magnetism

So far, our discussion about magnetism has been independent of the dimensionality of the material. However, the dimensionality of the 2D magnetic materials makes magnetism different from the 3D counterparts. One of the most notable differences arises as a consequence of the Mermin–Wagner theorem.⁴ This theorem states that for a 2D material with infinite size and modelled with an isotropic Heisenberg model with finite range interactions, long range magnetic order is not possible at nonvanishing temperature. This theorem is sometimes formulated saying that long range magnetic order is not possible for a 2D material with infinite system size and with any continuous symmetry.^13,58 However, this theorem does not contemplate what happens when anisotropy is present.

What we know as of today is that anisotropy can stabilize magnetic order at nonvanishing temperature as shown for CrI₃ experimentally¹ 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.⁶ Monolayer Fe₃GeTe₂ presents a similar feature with a ferromagnetic phase driven by strong out of plane anisotropy up to 130 K.³ 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.³⁷ 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.

3. Introduction to DFT

The hydrogen atom was the first important benchmark to illustrate the success of quantum mechanics and the Schrödinger equation. Regretfully, more complex systems such as heavier atoms or solids were still intractable at the time, mainly because of the many body nature these systems. As Dirac said, the fundamental laws necessary for the mathematical treatment of a large part of physics and the whole of chemistry are thus completely known, and the difficulty lies only in the fact that application of these laws leads to equations that are too complex to be solved. In the particular case of a solid, solving the Schrödinger equation for a wavefunction Ψ(r₁…,r_N) with N spatial coordinates (3N scalar variables) is completely out of reach, even computationally. This is why the birth of Density Functional^{26,27,40,59,60} (DFT) is one of the most important milestones in condensed matter physics since it provides an exact theory to describe many-body systems with interacting particles.²⁷

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 n₀(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 n₀(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:²⁷

where T[n] is the kinetic energy contribution, E_int[n] is the contribution of the electron–electron interaction, E_II is the repulsion between cores and the term of the integral corresponds to the external potential contribution (core-electron interaction).

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 n₀(r). As a result, knowing the ground-state electronic density n₀(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.²⁶ 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 n₀(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:

where T_s[n] represents the kinetic energy of the particles of the auxiliary system, E_xc[n] represents the contribution of the many-body effects written as a functional of the density and E_Hartree[n]:

Is the classical interaction energy of the charge density with itself. Since (21) is just (20) rewritten, by making them equal and solving for E_xc[n]:

E_xc[n] = T[n] − T_s[n] + E_int[n] − E_Hartree[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 gas^26,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.⁶⁸

The practical development of the Kohn–Sham approach leads to a set of equations:

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 configurations⁶⁹ 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 tutorial⁷⁸ and will not be elaborated here, we consider it critically important to address the delocalization problem inherent to all approximations of the exchange–correlation functional^79–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

3.1. Improving the localization by using DFT+U

It was shown back in 1982 that the exact exchange–correlation functional follows the constraint that the total energy behaves in a linear piecewise manner as a function of the total number of electrons.⁸⁸ It has been covered already in previous reviews^79–82 how the violation of this constraint makes approximate-exchange correlation functionals delocalize electrons exceedingly, leading to poor description of the electronic structure for certain systems, huge underestimations of the bandgap or even the contradictory prediction of a metal instead of an insulator. Some typical examples of this failure are metal oxides such as NiO,^84,89,90 FeO and MnO.^91–93 These examples have all in common that the partially filled strongly localized d or f shells play an essential role in the electronic structure. The approximate exchange–correlation functionals lead then to a particularly bad description of these l-shells. In order to correct this prominent excess of delocalization, DFT+U was invented.⁹⁴ Inspired by the Hubbard model,⁹⁵ DFT+U aims to favor localization of the electrons by introducing an energy penalty when the orbitals have fractional occupation. This way, DFT+U favours cases in which the orbitals are either with one electron or empty, avoiding the cases in between. This is done by rewritting the exchange–correlation functional as:^71,86where E_DFT[ρ(r)] is the total energy given by the used exchange–correlation functional, is the Hubbard-like correction and E_dc[n^Iσ] is the double counting term that aims to remove the contribution of the corrected orbitals from E_DFT[ρ(r)]. By construction, the computational cost of the DFT+U approach is practically the same as a usual DFT calculation. I is an index running over the atoms with correlated/localized electrons, and the indices m and m′ run over the localized states of atom I. Most of the times, m and m′ run over states with the same angular momentum quantum number l i.e. they belong to the same l-shell. are the occupation numbers of the localized states so that adds an energy correction depending on how populated they are. The details on how is calculated are code-specific, but a common way is projecting the Kohn–Sham eigenstates {ψ^σ_ki} on a basis set of localized functions {φ^I_m}.^71,86

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,⁸⁶ being implemented in most of the DFT packages (see Table 1).

Table 1 List of popular Density Functional Theory (DFT) codes with capabilities and features. These features have been compiled after searching through each code's manual

DFT code	License	Basis set	sc-DFT	GBT	vdW corrections	DFT+U	U obtention	DFT+U+V	V obtention
VASP^150,151	Commercial	Plane-waves (pseudopotentials)	Yes^166	Yes	DFT-D2^167	Yes	Yes^100	No	No
					DFT-D3^168,169
					DFT-D4^170
					Method of libMBD^171
					Tkatchenko–Scheffler^172–174
					MBD@rSC^175,176
					MBD@rSC/FI^177,178
SIESTA^179	GPL (academic)	Localized atomic orbitals (pseudopotentials)	Yes^134 (but not public)	Yes (but not tested)	DRSLL^180	Yes^71	No	No	No
					LMKLL^181
					KBM^182
					C09^183
					BH^184
					VV^185
FLEUR^152,153	GPL (academic)	FP-LAPW (all electron)	Yes	Yes^186	DFT-D3^168,169	Yes^187	No	No	No
					DRSLL^180
Quantum ESPRESSO^114,118	GPL (academic)	Plane-waves (pseudopotentials)	Yes	No	DFT-D2^167,188	Yes	Yes^100,105	Yes^87	Yes^100,105
					DFT-D3^168
					Tkatchenko–Scheffler^172
					MBD^176
					XDM^189,190
OPENMX^154,155	GPL (academic)	Localized atomic orbitals (pseudopotentials)	Yes	Yes^157,158	DFT-D2^167	Yes^191,192	No	No	No
					DFT-D3^168,169
GPAW^146,147	GPL (academic)	Plane-waves (pseudopotentials)	Yes	Yes	DRSLL^180,193	Yes	No	No	No
					LMKLL^181
					BH^184
					KBM^182
					C09^183
ABINIT^194	GPL (academic)	Plane-waves (pseudopotentials)	Yes^195	No	Yes^196	Yes^197	Yes^100,108	No	No
					DRSLL^180
					DFT-D2^167
					DFT-D3^168
CASTEP^198	Commercial	Plane-waves (pseudopotentials)	Yes	Yes	Tkatchenko–Scheffler^172	Yes	Yes	No	No
					MBD@rSC^175,176
					DFT-D2^167
					DFT-D3^168
					OBS^199
					XDM^200
ELK^201	GPL (academic)	FP-LAPW (all electron)	Yes	Yes	No	Yes	No	No	No
WIEN2K^202	Academic (paid)	FP-LAPW (all electron)	No	No	DRSLL^180	Yes	No	No	No
					LMKLL^181
					SCAN+rVV10^203

---

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:

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 {U^I}, 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 U^I 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.⁸⁶ 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.⁸⁶

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 approximation^106–108 (cRPA) and some machine learning approaches^109,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.⁸⁸ 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 HP¹⁰⁵ available as part of the DFT code QUANTUMESPRESSO.¹¹⁴ This process is shown in Fig. 3. The structure is initially relaxed for an initial value of U = U_in, which can be zero. Then, the ground state is obtained and a value of U = U_out is calculated for that ground state. Afterwards, a new ground state is calculated with the obtained U_in = U_out. For this new ground state, a new value of U_out is calculated. The process is then repeated until the difference between the input U and the output U is less than a threshold.

Fig. 3 Scheme of the algorithm to obtain the U parameter in a self-consistent procedure.

The DFT+U approach serves to describe more localized regimes. Following a similar idea, DFT+U+V was born:⁸⁷ 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

where the star over the sum denotes that the sum is taken over all the neighbours up to a given shell. The value V^IJ is the V Hubbard parameter between shells I and J. The Hubbard V favors states with localization in the neighbouring atoms, positively affecting hybridization that is usually suppressed by the exceed of localization introduced by the Hubbard U in the atomic centers.¹⁰⁴

The DFT+U+V approach was used for NiO, Si and GaAs in its original publication.⁸⁷ 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 used^{110,115–117} and the obtention of V is fully automated thanks to the code HP¹⁰⁵ which is part of QUANTUM ESPRESSO^114,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 Hubbard¹²³ 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 system¹²⁴ 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 effects⁴⁷ 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.

3.2. Improving the localization by using hybrid functionals

We already discussed when introducing DFT+U how common functionals tend to delocalize charge and suffer from the self-interaction error. A different approach to compensate the self-interaction error and to improve the performance of the functionals is by using the so called hybrid functionals.⁴¹ Hybrid functionals treat the exchange–correlation energy following:

E^hybr_xc = a₀E^exact_x + (1 − a₀)E^DFT_x + E^DFT_c, 0 ≤ a₀ ≤ 1. (28)

where E^exact_x is the exact exchange, E^DFT_x is the exchange energy coming from a traditional DFT functional (LDA, GGA…) and E^DFT_c is the correlation energy which comes from the DFT exchange–correlation functional. a₀ is a parameter to be tuned to balance between the exact exchange contribution and the usual DFT exchange. The exact exchange energy is taken from a Hartree–Fock-like equation:where φ_iσ are the Kohn–Sham eigenstates labelled by band index and spin. Therefore, hybrid functionals act over all electrons where as in DFT+U, the correction was applied only on those electrons considered as correlated or localized (usually d and f electrons). Moreover, DFT+U assumes an orbital independent value of the interaction energy among electrons with the same l-number. This assumption comes from the idea that localized states should retain atomic character and therefore, spherical symmetry.⁸⁶ However this assumption can break down in cases where there is a mechanism capable of removing the energy degeneracy of those states such as situations with strong crystal field or strong SOC.^86,127 On the other hand, hybrid functionals treat those electrons with the same l individually.

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 a₀ 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.⁸⁶

It is important to remark that a higher amount of exact exchange (increasing the a₀ parameter) does not necessarily imply a better result. This is because the exchange contribution from DFT, E^DFT_x 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 a₀ = 1 can be critical.

The mixing scheme of hybrid functionals can be extended even further by letting the mixing coefficient be spatially dependent:

where ε^exact_x,σ(r) is the exact exchange–correlation energy density:ε^DFT_x,σ(r) is the DFT exchange–correlation energy density and g_σ(r) is the so called local mixing function (LMF), which follows 0 ≤ g_σ(r) ≤ 1.

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.

3.3. van der Waals corrections in 2D magnetic materials

One of the issues of semilocal exchange–correlation functionals is their inability to account for long range correlation interactions such as London dispersion,^130,131 resulting in an overrepulsive potential energy between atoms as a function of the interatomic distance.¹³⁰ This interaction is one of the van der Waals type interaction and it can be of special relevance in condensed systems.

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’:¹³¹

E_DFT+vdW = E_DFT + E_disp (32)

where the term on the left side is the total energy, E_DFT is the total energy obtained from standard DFT and E_disp is the energy correction term to account for the long range corrections. This energy contribution is calculated as:where the subindex AB means that the quantity refers to the atoms A and B. This way, the energy is the sum of potentials proportional to r_AB⁻ⁿ with dispersion coefficient C_n,AB and damping function f_damp,n(r_AB). The damping functions appear because every energy contribution applies only for long range. The damping function then makes sure that the contribution vanishes at short distances, avoiding potential overbinding effects.¹³⁰ Eqn (33) is obtained from the second order perturbation theory of the correlation energy after writing the Coulomb potential in a multipole expansion.¹³⁰ This multipole expansion is the reason for the appearance of the sum , being the term with n = 6 the one corresponding to the dipole–dipole interaction.

The details of the calculations of the terms C_n,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.

3.4. Obtaining the magnetic order of the ground state

Modelling a magnet within DFT is about finding its ground-state electronic density, but the electronic density and the spin order are not decoupled. The spin order influences the electronic density and structure. Therefore, a crucial step of the modelling process is focusing on finding the ground state magnetic order of the material under study.

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.¹³²

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 DFT¹³⁴ (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 optimization^111,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:¹⁴⁰

where q is the spin-spiral wavevector and k labels the state as done in the traditional BT. This result allows to simulate spin-spirals within the primitive cell. Its biggest caveat is that its derivation assumes that all directions are equivalent and therefore, the rotation axis can be set to the z-axis in such a way the spins behave as:where θ is the cone angle (angle between the precession axis and the spin direction). The problem of this assumption is that important effects such as SOC can destroy the equivalency of directions and therefore the GBT is no longer a fully valid approach. At least, the effect of SOC can be incorporated perturbatively.^28,142,143

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 spirals²⁸ 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 VASP^150,151 and FLEUR.^152,153 The localized atomic orbitals DFT code OPENMX^154–156 has also a GBT implementation.^157,158

3.5. Magnetic ground states in 2D magnetic materials

2D magnetic materials can show a variety of magnetic orders in their ground state (see Fig. 5). In ref. 29, 192 magnetic materials from C2DB^22,23 were studied with the GBT. 50 of them were found to be FM, 21 AFM, 34 commensurate non-collinear spin-spirals, 36 incommensurate spin-spirals and 15 chiral spin-spirals. This results exemplifies the many possibilities that can appear when trying to find the magnetic order of a 2D magnet.

Fig. 5 Classification of some important 2D materials in the monolayer limit.

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 as^19,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 CrX₃ (X = I, Cl, Br) or CrSBr monolayers. Other materials are well known by their intralayer-AFM such as the MPS₃ (M: Mn, Fe, Co, Ni). Fe, Co and Ni compounds are examples of the zigzag AFM introduced in Fig. 4, by the other side MnPS₃ presents an example of a Néel AFM. Also these materials provide a good example of different AFM (FePS₃) and FM (CoPS₃, NiPS₃ and MnPS₃) 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 CrCl₃ 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 CrI₃ is an out of plane ferromagnet up to 45 K.¹ 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.⁶ 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 results¹⁶³ only when shape anisotropy is considered.¹⁶⁰

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 (MX₂) such as the Ni or Co compounds.

In Fig. 6 we show one of the results of ref. 29 for CoI₂. 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 CoI₂. 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 Fe₃GeTe₂ and Fe₃GaTe₂. 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.

3.6. What DFT codes can I use?

There are many DFT codes available and not all of them offer the same relevant features in the context of magnetism. That is why we summarize the availability of some of the features we discuss in this review in Table 1 for some of the most known and widespread DFT packages. These features have been compiled by inspection of the manuals of the corresponding codes.

4. Obtaining exchange parameters

Once the magnetic ground state is known, the parameters of an atomistic magnetic Hamiltonian can be obtained from DFT. In this review, we will describe in simple terms the usual approaches to do so.

4.1. Energy mapping method

The energy mapping method is the most popular technique to calculate exchange interactions.^{6,14,16,17,30,204–207} The fundamental approximation behind this method consists in dividing the contributions to the total energy into two components:where E₀ is just a reference energy which is assumed to be independent of the spin configuration, χ_i denotes a spin configuration and corresponds to the magnetic spin-dependent Hamiltonian for that specific spin configuration. It is obvious that this assumption is quite strong: a change in the spin configuration can also change the energy contributions of the bands, for example, and hence, the decoupling in (36) is very naive. Nevertheless, the assumption works well for insulators and semiconductors.

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 collinear³⁰ and noncollinear configurations.⁴⁷

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).³² 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:

where S is the value of the spin (only one magnetic ion per unit cell was considered) J* is the exchange constant considering the FM and AFM configurations and β is a factor depending on the type of latticem, β ≥ 0. Hence, it is clear that the bigger the value of S, the smaller the correction is. Nevertheless, the exact expression of the parameter will change between different cases.

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.

4.2. Fitting the spin-spiral spectra

Another way to calculate the exchange parameters is by fitting the spin-spiral energy spectra obtained after using the GBT. This approach is perfectly illustrated in ref. 144, 145 and 209. By inserting eqn (35) in the magnetic Hamiltonian, an expression of the energy as a function of q can be obtained. Then, the spin-spiral spectra can be fitted with the parameters of the atomistic magnetic Hamiltonian.

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,³⁷ 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.

4.3. Approaches based on the magnetic force theorem and the LKAG approach

The approaches that can circumvent the problems of the previous methods are those based on the so called LKAG (Liechtenstein, Katsnelson, Antropov and Gubanov) approach and the magnetic force theorem.^210,211

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.³¹ When the perturbation is small enough, the energy variation can be calculated using the so called magnetic force theorem:³¹

where is the density of states of the single-particle states (in our case, the Kohn–Sham states), E_F is the Fermi energy and δZ is the total number of electrons, which is equal to zero for changes in the spin configuration. The intuitive idea of the approach is to write the left-hand side of the equation in terms of the magnetic Hamiltonian parameters when considering two infinitesimal spin variations at sites i and j i.e. δE_ij. The right hand side of the equation can be expanded in terms of Green function's^212–214 of a Hamiltonian obtained in the last iteration of the DFT calculation. The spin rotations can be introduced in this Hamiltonian by rotating some of its parts²¹⁵ and then, expressions for the exchange parameters can be deduced.²¹⁴ This rotation can be done easily only if the basis sets used in the DFT calculation is of localized character. This is why practical implementations usually require LCAO calculations or a Wannierization although a plane-waves implementation of the magnetic force theorem exists.²¹⁶ The LKAG to this day has been used to compute symmetric exchange,^31,217 anisotropic exchange,²¹⁴ DMI interaction,^213,218 biquadratic exchange²¹⁹ … The formalism behind this approach is too complex to be introduced here but the interested reader can find more information in extensive reviews.^210,211

4.3.1. Available codes implementing LKAG methods. Probably, the most famous package to calculate exchange constants is TB2J.²¹⁴ This package takes Wannier90²²⁰ functions or LCAO DFT results to construct a tight-binding Hamiltonian in order to calculate interactions such as isotropic exchange, anisotropic exchange and DMI. TB2J can use directly the results of Siesta and OPENMX, while for plane waves or FLAPW codes, localized functions have to be constructed via Wannier90.

Jx is another code based on the same LKAG and magnetic force theorem approach.²²¹ 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,²²² capable of calculating the symmetric exchange from Siesta calculations. Nojij is being extended to a package named Grogu²²³ which is capable of calculating symmetric exchange and SIA. Unfortunately, to our best knowledge, Grogu is not yet publicly available.

4.4. Exchange mechanisms driving magnetic order in 2D magnets

As we have explained before, the long range magnetic order in 2D magnets with infinite system size has to be stabilized by anisotropy. The stabilization can be done mainly by two different means: single ion anisotropy or anisotropic exchange. For example, CrI₃ is an out of plane ferromagnet, but its magnetic order is stabilized by two-ion anisotropic exchange along the z direction since the Cr atoms SIA is very small.^6,25 The same behaviour has been obtained^25,224 for CrBr₃. Contrary to these cases, the in-plane ferromagnetism of monolayer CrSBr is stabilized by SIA since its anisotropic exchange results negligible from DFT calculations.⁹⁹ However, SIA is not the only interaction driving its magnetic order since dipolar interactions have been proved to affect it.⁵¹

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”.

5. Calculation of the Curie temperature

Several methods exist for calculating the critical temperature of a material based on a known magnetic Hamiltonian, each with varying levels of accuracy, applicable regimes, and distinct advantages and limitations. A comprehensive understanding of these methods is valuable for interpreting results effectively and, where appropriate, for combining techniques to improve accuracy or achieve rapid estimations of the critical temperature. In this section, we will introduce the most famous ones for ferromagnetic materials (although most of them can be extended for antiferromagnets as well), highlighting their limitations and how they can be use in synergy with others.

5.1. Using mean field theory

Mean Field Theory (MFT) for ferromagnetism^32,225 introduces the concept of an effective, uniform magnetic field acting on all ions in a lattice, referred to as the molecular field. This molecular field is assumed to be directly proportional to the magnetization of the system, expressed as , where λ is a proportionality constant. In MFT, this molecular field affects all ions uniformly, allowing the ferromagnetic system to be treated as a paramagnet subject to an internal magnetic field equivalent to the molecular field.^36,225 The primary function of the molecular field is to align all spins parallel to it, effectively substituting the ferromagnet's exchange coupling. For the model to be self-consistent, λ must be positive.

In a lattice of identical magnetic ions, the critical temperature with MFT is as follows:²²⁶

where represents the sum of exchange constants between nearest neighbors. Some sources, such as,^15,32 include an additional factor, S(S + 1), to account for the spin of the ions. This factor is omitted here since, in our convention, the exchange constants are expressed in energy units, and spins are treated as unit vectors.

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.

5.2. The formulas for the 2D lattices with the Ising model

Consider the Ising model with no anisotropy and only nearest-neighbor exchange interactions. This model depends solely on the exchange parameter J, assumed constant across all nearest neighbors, and on the number of nearest neighbors, which is determined by the lattice. Consequently, the critical temperature in this model is a function only of the lattice type and the value of J. The relationship between critical temperature, lattice structure, and J has been investigated both numerically²²⁷ and analytically.²²⁸ Table 2 provides critical temperatures for common lattice types.

Table 2 Relation between the critical temperature and the J exchange constant with the Ising model for several lattices

Lattice
Square	2.2692
Hexagonal	1.5186
Honeycomb	3.6410

---

The Mermin–Wagner theorem⁴ 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 CrI₃ was found to exceed experimental values significantly – more than twice the experimental value for monolayer CrI₃ and nearly four times that of the bulk material. Similarly, ref. 230 reported a critical temperature of 161 K for CrI₃ using the Ising model, far above the experimentally observed 45 K.¹

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.

5.3. Methods based on calculating the magnetization with temperature

It is well established that ferromagnets exhibit zero magnetization at and above the critical temperature. Consequently, the critical temperature can be estimated by analyzing the behavior of the system's magnetization as a function of temperature. In this section, we outline various theoretical approaches and methodologies for calculating the magnetization–temperature relationship, spanning multiple levels of theory. In the following section, we will discuss methods for fitting the resulting data to extract precise estimates of the critical temperature.

5.3.1. Using spin-wave theory. Spin-wave theory studies magnons, quasiparticles that describe the collective motion of the spins in a wave-like behaviour. Magnons can appear as thermal excitations, and they are capable of breaking the long range magnetic order in magnetic materials. Hence, studying their population in the material can lead to the value of the critical temperature. One important remark is that in this section, the spins will be considered as quantum operators instead of unit vectors.

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

where the subindex ai represents the magnetic ion a in the unit cell i and b^†_ai (b_ai) are the bosonic creation (annihilation) operators. This representation is the most common in the literature although there are some others that serve the same purpose such as the Dyson–Maleev,^233–235 the finite differences²³⁶ representation²³⁷ or the differential equations approach in ref. 238. This representation can also be generalized to antiferromagnets.³²

The physical intuition of these operators is that the spin state |S,m_S = S 〉 is mapped to the vacuum bosonic state. Then, the states with lower m_S are obtained as excitations of the vacuum bosonic state i.e. acting with b^†_i creates a boson and decreases m_S by one. As a consequence of the (2S + 1) possible values of m_S, 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⁻] = 2S^z and [S^zS^±] = ±S^±. The second one, to keep the number of bosons to be 2S at most. Consider the state |S,m_S = −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:

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 dominant^6,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 manner^6,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 theorem²⁴⁰ 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) ≈ Δ₀ + ρk² (42)

In this context, Δ₀ denotes the spin-wave gap, while ρ represents the spin-wave stiffness. The spin-wave gap Δ₀ 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.⁶

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;N_magnons) (43)

But since magnons are spin 1 bosons, the obey Bose–Einstein statistics and the total number follows the Bose–Einstein distribution:^15,232

where the sum runs over all the k points in the first Brillouin zone. For a small enough spacing between k states, we can approximate the first Brillouin zone by a continuum. By doing so, the number of magnons is:^32,225where the prefactor is to account for the density of k states i.e. every k state occupies an area of in the reciprocal space. So in the end, the total number of magnons depends simultaneously on the energy dispersion. This forces the need of solving E(k) and (44) or (45) in a self-consistent manner when nonlinear spin-wave theory is used.

Since magnons are spin 1 bosons, we can calculate the total magnetization (assuming a single magnetic ion per unit cell) as:

where M₀ is the saturation magnetization of the system. This expression must be consequently modified when more than one magnetic atom is considered in the unit cell. The analogous expression for a continuous approximation is:

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.

5.3.1.1. Alternative representations to the Holstein–Primakoff representation. In the non-linear spin-wave theory approach, the ladder operators were expanded using the Taylor expansion of the square root (41). The expansion was truncated, keeping terms with up to four bosonic operators that were approximated afterwards in a Hartree–Fock like approach. This truncation leads to neglecting magnon–magnon interactions.

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 differences²³⁶ instead of a Taylor expansion. They show that the square root can be written as:

where n̂ = b^†b is the number operator, Δf(x) = f(x + 1) − f(x) and Δ^k = (Δ)^k. Consider an r order expansion of (48) and let us call S_r⁺ and S_r⁻ the resulting expansions of the ladder operators. The authors found that the resulting commutation relations are:where n̂^(k) = n̂(n̂ − 1)…(n̂ − k + 1). Note that n̂^(r+1) vanishes unless there are at least r + 1 bosons.²³⁷ Hence, the first commutation relation in (49) yields the correct result when 2S ≤ r and with that condition, the operators given by an order r truncation respect the original spin commutation relations. An order r truncation of the Taylor expansion yields and hence, the commutation relation does not improve when r increases, contrarily to the finited differences expansion.²³⁷

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.

5.3.2. Using Metropolis Montecarlo (MMC). Consider any of the magnetic Hamiltonians discussed in this review. For any lattice, the number of possible collinear spin microstates for the system is 2^N, where N is the number of magnetic ions in the lattice. Therefore, if we aim to find a configuration that minimizes the energy, examining each configuration individually becomes infeasible.

The Metropolis Monte Carlo algorithm^33,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:

where the sum runs over all microstates with energy ε. Here, is the partition function of the system:where the index i runs over all possible microstates of the system. Thus, computing the partition function and is impractical for real systems. However, we can still determine whether one configuration is more probable than another by calculating the ratio:

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.²⁰⁵ This can be confirmed with the examples shown in Table 3.

Table 3 Critical temperatures of chromium based 2D ferromagnets obtained with Ising based Monte Carlo vs. experiment

Material	TcIsing/K	Tc (experimental)/K
CrI3	107^244	45^245
CrBr3	86^244	27,^245 34^246
CrCl3	66^244	16^245 (bilayer)

---

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.²⁴⁷ 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.¹⁴ 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.

5.3.2.1. Empirical rescaling of the temperature. As explained previously, the MMC method yields accurate values for the magnetization for high and close to T_c temperatures. The values it provides for the magnetization at low temperatures are not accurate and this can turn into an obstacle when calculating quantities that depend on the magnetization or its derivatives. Evans et al. managed to solve this issue with a simple rescaling of the temperature of the MMC.²⁴⁹ Since MMC considers spins as classical particles, the real temperature we would need to obtain the same results (making quantum effects negligible) in a laboratory should be higher. The quantitative relation they suggest is:where T_MMC is the temperature of the MMC and T_experiment represents the real temperature that we would need in an experiment. The parameter α is just the effective rescaling factor that can be obtained from experimental measurements of the magnetization with temperature. Eqn (53) is called the quantum rescaling of the temperature in the MMC.

The results of applying this method are shown in Fig. 7 for 2D CrI₃.²⁰⁵ 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 CrI₃ 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).

5.3.2.2. Available codes implementing Monte Carlo spin dinamics. There are already available codes implementing MMC for spin lattices. Some examples are VAMPIRE,^{33,249,251,252} SPIRIT^253,254 and FIDIMAG.²⁵⁵ The code UPPASD²⁵⁶ incorporates the thermal relaxation by Langevin Dynamics²⁵⁷ instead of MMC, but this approach leads to analogous results for the critical temperature.

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.²⁴⁷ 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 AiiDA^120–122 plugin package,²⁵⁴ 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.

5.4. The Green's functions method

The magnetization with temperature can be calculated by employing Green's functions. Let us define the retarded Green Function of Operators A(t) and B(t′) in the Heisenberg picture:^34,225

G_AB(t,t′) := 〈〈A(t); B(t′)〉〉 := −iθ(t − t′)〈[A(t),B(t′)]〉 (54)

where θ(t) is the step function. And we are writing:

Exploiting the properties of these functions and their Fourier transforms, the following equation can be obtained:

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.³⁴ 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.³⁴ 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 S_i⁻ and S_i⁺ and let us denote G_ij(ω) = 〈〈S_i⁻;S_i⁺〉〉. We are omitting the hat over the spin operators on purpose for the sake of simplicity.

where we have used:

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:

The next term to compute is:

With this equality we can write (57) as:

In order to simplify further, we need to introduce a decoupling. The decoupling known as Random Phase Approximation (RPA) goes as follows:

And the justification for this approach is purely heuristic.⁸⁶ Note that implicitly, the RPA is setting S^z = 〈S^z〉, 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 (〈S^z_i〉 = 〈S^z_j〉 = 〈S^z〉), one can arrive to the system of equations:

ω_k^RPA = B + 〈S^z〉(J₀ − J_k) (65)

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

Another famous decoupling scheme is the Callen decoupling^34,258

〈〈S^z_iS_l⁺;S_j⁻〉〉 ≃ 〈S^z_i〉〈〈S_l⁺;S_j⁻〉〉 − α〈S_i⁻S_l⁺〉〈〈S_i⁺;S_j⁻〉〉 (68)

where the parameter α is commonly set to α = 〈S^z〉/2S². The Callen decoupling leads to a system of equations similar to that for the RPA:

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}

5.4.1. Models to fit the magnetization. We already know that from the critical temperature and above, ferromagnets present zero magnetization. Therefore, the critical temperature can be estimated by observing the evolution of the magnetization of the system as a function of the temperature. For a fully systematic calculation of the Curie temperature, it is essential to have a model to fit the magnetization curve. The dependence of the magnetization at low temperatures can be obtained with linear spin-wave theory. Using eqn (47). For example, for a 3D material with a dispersion for low leads to the famous Bloch equation^32,225 and including more terms in the expansion of results in the appearance of terms proportional to T^5/2 and T^7/2.³²

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.²⁶⁰ In fact, for these magnetic Hamiltonians, the magnetization close to the critical temperature is well described by the equation:

where M_s is the saturation magnetization, the value of the critical exponent depends on the dimensionality of the system and the dimensionality of the spin lattice. For a 2D material, the exponent is predicted to be β = 1/2 for an Ising Hamiltonian.^13,261 There is no critical exponent for 2D materials and a Heisenberg Hamiltonian since the Mermin–Wagner theorem forbids the existence of long range magnetic order in this situation (for an infinite system size). However, there are already known critical exponents for the 3D counterpart.

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:²⁴⁹

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.²⁰⁵

5.5. Critical temperatures fitted from Montecarlo results

MMC is computationally expensive and its cost increases with the number of atoms and the number of interactions considered. However, its physical input is just the magnetic Hamiltonian and therefore, its results should depend only on its parameters. Two different systems with the same magnetic Hamiltonian should give the same result. This motivates the goal of having a closed formula relating the critical temperature with the Hamiltonian parameters, even if it is empirical. This was the goal of Torelli and Olsen.¹⁵ Using the following Hamiltonian that includes both SIA and exchange anisotropy for a single magnetic ion per unit cell:they carried out nonlinear spin-wave theory and MMC simulations for different ratios of the Hamiltonian parameters and different 2D lattices. Then, they fitted the results to an empirical model. The results for the case B = 0 for a square lattice are shown in Fig. 8.

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 T_c, 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:

where N_nn is the number of nearest neighbours, γ = 0.033 and T^Isingratio_c is the value of for the Ising model and the corresponding lattice (see Table 2). It is important to remark that the empirical formula (74) applies for systems with the given lattices and the Hamiltonian in (73) with B = 0 and with a single magnetic ion per unit cell. Still, the formula serves as a cheap first approximation of the critical temperature and can be combined to obtain a rough estimation of the necessary upper bound of an MMC simulation.

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:

5.6. Critical temperatures fitted from different methods

At this point, many ways of calculating the Curie point of 2D ferromagnets have been introduced and one may wonder how compatible the results of the different methods can be. In the each method's section we have discussed their limitations in a qualitative manner, but a more ambitious scientist might still want to have a quantitative relation or a way to compare the output Curie point of the different methods. We will discuss the work of Tiwari et al.¹⁴ for this purpose. In this work, they calculated the Curie point both with computational methods and with a closed formula they propose for the 2D magnets featured in the Computational 2D Materials Database^22,23 (C2DB). The three methods they used were: non-linear spin-wave theory, the Green's functions method and MMC.

The magnetic Hamiltonian they proposed to model the ferromagnets is:

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 J^xx_ij = J^yy_ij = J_ij(1 − Δ_ij) and J^zz_ij = J_ij(1 + Δ_ij).

Tiwari et al. looked for 2D magnets in C2DB with positive first neighbours exchange constant J_ij := 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 T_c values, they fitted the results to the following model:

where θ = 1 for the Green's functions method and spin-wave theory and θ = 0 for the MMC. α₁ and α₂ are parameters to fit as a function of J,S and Δ_nn for each method and type of lattice. The results they obtained for these two parameters are shown in Table 4 for the Honeycomb, hexagonal and square lattices. The authors state that this formula is only valid when Δ_nn ≤ 0.2 since it is based on the theoretical approach in ref. 262.

Table 4 Resulting parameters for eqn (78) across different lattices and methods. Data extracted from ref. 14

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 T_c 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.52S², with S = 3/2) and the quantum mean field (1/3[S(S + 1)N_NN], N_NN = 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.

5.7. Brief review of the methods and conclusion

We would like to close this section with a brief discussion about the limitations and the use cases of the methods explained above.

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 T_c 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 T_c 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 T_c temperatures as explained for Fig. 7. However, this unusual low magnetization is not a problem if we are only interested on obtaining T_c.

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 S^z = 〈S^z〉, 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 regime³⁴ 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.

6. Summary and outlook

Despite remarkable advances, simulating 2D magnetic materials from first principles remains a formidable challenge due to the extreme sensitivity of magnetic interactions to subtle structural, electronic, and computational parameters. The reduced dimensionality of these systems introduces unique theoretical and practical difficulties, particularly in capturing the intricate anisotropy mechanisms, and often delicate magnetic ordering. Unlike their 3D counterparts, 2D magnets rely critically on anisotropic interactions – such as single-ion anisotropy and different types of anisotropic interactions – to stabilize long-range magnetic order against thermal fluctuations, in line with the Mermin–Wagner theorem. However, these energy contributions typically lie in the sub-meV regime, pushing the limits of accuracy for standard DFT methods. Furthermore, the layered nature of many 2D materials introduces complex indirect exchange pathways mediated by ligands or chalcogen atoms, requiring careful electronic structure modeling. While approaches like DFT+U and hybrid functionals have improved our ability to capture localized electron behavior, the absence of universal, ab initio protocols for determining magnetic ground states (especially noncollinear or incommensurate orders) remains a major bottleneck. Addressing these challenges will require more systematic treatments of magnetic anisotropies, automate the exploration of large configuration spaces (e.g., via generalized Bloch theorem or spin-constrained DFT), and expand beyond traditional collinear models. Future progress will hinge on combining accurate, systematically parameter-free methods (e.g., self-consistent DFT+U/V, hybrid functionals, and generalized Bloch theorem-based approaches) with machine learning-driven optimization and high-throughput screening will be crucial to identify new candidate materials. Coupling these advances with high-throughput workflows and experimental feedback will be key to uncovering novel 2D magnets with robust and tunable properties, particularly those that can operate above room temperature for next-generation spintronic and quantum devices.

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. 11 Sketch of the worfklow to follow when modelling 2D magnetic materials from DFT to the calculation of their critical temperature. All the steps in this workflow are covered throughout this review paper.

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 CrI₃, 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.

Table 5 Compilation of calculated critical temperatures from DFT results

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)	45^1	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	PBE^23	41	Green's Functions		1.025^23	14
	GPAW	PBE^23	31	Montecarlo		1.025^23	14
	GPAW	PBE^23	23	NLSW		1.025^23	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	418^263 (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	146^162	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 34^270	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	PBE^23	387.6	Green's Functions		21.96^23	14
	GPAW	PBE^23	371.1	Montecarlo		21.96^23	14
	GPAW	PBE^23	263.4	NLSW		21.96^23	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	240^273 367 (few layers)^274	57.18	165
	QE	GGA	320	Montecarlo		35.52	164
Fe3GeTe2	Wien2k	LDA+U (U = 4)	152	Montecarlo	130^3	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 VTe₂, 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.

Author contributions

Jose-Hugo García was responsible of the conceptualization, funding acquisition, project administration, supervision, validation, writing, reviewing and editing. S. Roche participates in the project administration and editing. Stephan Roche improved the manuscript and contributed to the broad scope discussion. Dorye. L Esteras and Jaime Garrido Aldea have contributed with the investigation, the original draft preparation and have reviewed and edited the original draft.

Conflicts of interest

The authors declare that there are no conflicts of interest that could have influenced the content, analysis, or conclusions of this work.

Data availability

No new primary research data, software, or code were generated or analyzed in this review. All data supporting the findings discussed in this article are derived from previously published studies, which are cited appropriately in the text and reference section. Additionally, images included in this review have been reproduced from other sources after obtaining the necessary copyright permissions.

Acknowledgements

Jaime Garrido Aldea acknowledges financial support from MICIU grant no. PRE2022-103460, funded by CEX2021-001214-S-20-8/MCIU/AEI/10.13039/501100011033 and by FSE+. DLE acknowledges financial support from MICIU and the European Union for the next-generation grant PRTR-C17.I1.

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.

References

1. B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. Jarillo-Herrero and X. Xu, Nature, 2017, 546, 270–273.
2. C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, Z. Q. Qiu, R. J. Cava, S. G. Louie, J. Xia and X. Zhang, Nature, 2017, 546, 265–269.
3. Z. Fei, B. Huang, P. Malinowski, W. Wang, T. Song, J. Sanchez, W. Yao, D. Xiao, X. Zhu, A. F. May, W. Wu, D. H. Cobden, J.-H. Chu and X. Xu, Nat. Mater., 2018, 17, 778–782.
4. N. D. Mermin and H. Wagner, Phys. Rev. Lett., 1966, 17, 1133–1136.
5. M. Jafari, W. Rudziński, J. Barnaś and A. Dyrdał, Sci. Rep., 2023, 13, 20947.
6. J. L. Lado and J. Fernández-Rossier, 2D Mater., 2017, 4, 035002.
7. S. N. Kajale, J. Hanna, K. Jang and D. Sarkar, Nano Res., 2024, 17, 743–762.
8. E. Elahi, G. Dastgeer, G. Nazir, S. Nisar, M. Bashir, H. Akhter Qureshi, D. kee Kim, J. Aziz, M. Aslam, K. Hussain, M. A. Assiri and M. Imran, Comput. Mater. Sci., 2022, 213, 111670.
9. Y. Khan, S. M. Obaidulla, M. R. Habib, A. Gayen, T. Liang, X. Wang and M. Xu, Nano Today, 2020, 34, 100902.
10. T. N. Kobernik and A. I. Kartsev, J. Phys. Chem. Lett., 2024, 15, 12151–12155.
11. C. Chowde Gowda, A. Kartsev, N. Tiwari, A. A. Safronov, P. Pandey, A. K. Roy, P. M. Ajayan, D. S. Galvão and C. S. Tiwary, J. Mater. Chem. C, 2024, 12, 18691–18703.
12. X. Huang, L. Zhang, L. Tong, Z. Li, Z. Peng, R. Lin, W. Shi, K.-H. Xue, H. Dai, H. Cheng, D. de Camargo Branco, J. Xu, J. Han, G. J. Cheng, X. Miao and L. Ye, Nat. Commun., 2023, 14, 2190.
13. Q. H. Wang, A. Bedoya-Pinto, M. Blei, A. H. Dismukes, A. Hamo, S. Jenkins, M. Koperski, Y. Liu, Q.-C. Sun, E. J. Telford, H. H. Kim, M. Augustin, U. Vool, J.-X. Yin, L. H. Li, A. Falin, C. R. Dean, F. Casanova, R. F. L. Evans, M. Chshiev, A. Mishchenko, C. Petrovic, R. He, L. Zhao, A. W. Tsen, B. D. Gerardot, M. Brotons-Gisbert, Z. Guguchia, X. Roy, S. Tongay, Z. Wang, M. Z. Hasan, J. Wrachtrup, A. Yacoby, A. Fert, S. Parkin, K. S. Novoselov, P. Dai, L. Balicas and E. J. G. Santos, ACS Nano, 2022, 16, 6960–7079.
14. S. Tiwari, J. Vanherck, M. L. Van de Put, W. G. Vandenberghe and B. Sorée, Phys. Rev. Res., 2021, 3, 043024.
15. D. Torelli and T. Olsen, 2D Mater., 2018, 6, 015028.
16. D. Torelli, H. Moustafa, K. W. Jacobsen and T. Olsen, npj Comput. Mater., 2020, 6, 158.
17. D. Torelli, K. S. Thygesen and T. Olsen, 2D Mater., 2019, 6, 045018.
18. D. Torelli, K. S. Thygesen and T. Olsen, 2D Mater., 2019, 6, 045018.
19. A. Kabiraj, M. Kumar and S. Mahapatra, npj Comput. Mater., 2020, 6, 35.
20. K. Wang, K. Ren, Y. Hou, Y. Cheng and G. Zhang, J. Appl. Phys., 2023, 133, 110902.
21. F. Giustino, J. H. Lee, F. Trier, M. Bibes, S. M. Winter, R. Valentí, Y.-W. Son, L. Taillefer, C. Heil, A. I. Figueroa, B. Plaçais, Q. Wu, O. V. Yazyev, E. P. A. M. Bakkers, J. Nygård, P. Forn-Díaz, S. D. Franceschi, J. W. McIver, L. E. F. F. Torres, T. Low, A. Kumar, R. Galceran, S. O. Valenzuela, M. V. Costache, A. Manchon, E.-A. Kim, G. R. Schleder, A. Fazzio and S. Roche, J. Phys.: Mater., 2021, 3, 042006.
22. S. Haastrup, M. Strange, M. Pandey, T. Deilmann, P. S. Schmidt, N. F. Hinsche, M. N. Gjerding, D. Torelli, P. M. Larsen, A. C. Riis-Jensen, J. Gath, K. W. Jacobsen, J. J. Mortensen, T. Olsen and K. S. Thygesen, 2D Mater., 2018, 5, 042002.
23. M. N. Gjerding, A. Taghizadeh, A. Rasmussen, S. Ali, F. Bertoldo, T. Deilmann, N. R. Knøsgaard, M. Kruse, A. H. Larsen, S. Manti, T. G. Pedersen, U. Petralanda, T. Skovhus, M. K. Svendsen, J. J. Mortensen, T. Olsen and K. S. Thygesen, 2D Mater., 2021, 8, 044002.
24. H. Zhang, Electron. Struct., 2021, 3, 033001.
25. D. Wines, K. Choudhary and F. Tavazza, J. Phys. Chem. C, 2023, 127, 1176–1188.
26. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
27. R. M. Martin, Electronic Structure: Basic Theory and Practical Methods, Cambridge University Press, 2004.
28. J. Sødequist and T. Olsen, 2D Mater., 2023, 10, 035016.
29. J. Sødequist and T. Olsen, npj Comput. Mater., 2024, 10, 170.
30. H. Xiang, C. Lee, H.-J. Koo, X. Gong and M.-H. Whangbo, Dalton Trans., 2013, 42, 823–853.
31. A. Liechtenstein, M. Katsnelson, V. Antropov and V. Gubanov, J. Magn. Magn. Mater., 1987, 67, 65–74.
32. K. Yosida, Theory of magnetism, Springer, 1996.
33. R. F. L. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. A. Ellis and R. W. Chantrell, J. Phys.: Condens. Matter, 2014, 26, 103202.
34. P. Fröbrich and P. Kuntz, Phys. Rep., 2006, 432, 223–304.
35. R. M. White, Quantum Theory of Magnetism, Springer, Berlin, Heidelberg, 3rd edn, 2007.
36. S. Blundell, Magnetism in Condensed Matter, Oxford University Press, 1st edn, 2001.
37. X. Jiang, Q. Liu, J. Xing, N. Liu, Y. Guo, Z. Liu and J. Zhao, Appl. Phys. Rev., 2021, 8, 031305.
38. F. Bloch, Z. Phys., 1929, 52, 555–600.
39. O. Madelung, Introduction to Solid-State Theory, Springer, Berlin, Heidelberg, 1978, vol. 2.
40. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
41. A. V. Arbuznikov, J. Struct. Chem., 2007, 48, S1–S31.
42. E. Ising, Z. Phys., 1925, 31, 253–258.
43. I. Dzyaloshinskii, et al., Sov. Phys. JETP, 1957, 5, 1259–1272.
44. T. Moriya, Phys. Rev., 1960, 120, 91.
45. T. Moriya, Phys. Rev. Lett., 1960, 4, 228.
46. S. Blügel and G. Bihlmayer, in Magnetism of Low-dimensional Systems: Theory, John Wiley & Sons, Ltd, 2007.
47. A. I. Kartsev, K. V. Obraztsov and P. V. Lega, J. Commun. Technol. Electron., 2023, 68, 1169–1190.
48. A. Kartsev, M. Augustin, R. F. L. Evans, K. S. Novoselov and E. J. G. Santos, npj Comput. Mater., 2020, 6, 150.
49. A. Kartsev, S. Malkovsky and A. Chibisov, Nanomaterials, 2021, 11, 150.
50. K. Yang, G. Wang, L. Liu, D. Lu and H. Wu, Phys. Rev. B, 2021, 104, 144416.
51. C. Boix-Constant, S. Mañas-Valero, A. M. Ruiz, A. Rybakov, K. A. Konieczny, S. Pillet, J. J. Baldoví and E. Coronado, Adv. Mater., 2022, 34, 2204940.
52. P. Bruno, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 6015–6021.
53. N. A. Spaldin, in Ferrimagnetism, Cambridge University Press, 2010, pp. 113–129.
54. M. Uhl, L. Sandratskii and J. Kübler, J. Magn. Magn. Mater., 1992, 103, 314–324.
55. K. Knöpfle, L. M. Sandratskii and J. Kübler, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 5564–5569.
56. E. Sjöstedt and L. Nordström, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 66, 014447.
57. S. Mankovsky, G. H. Fecher and H. Ebert, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 144401.
58. Y. Li, B. Yang, S. Xu, B. Huang and W. Duan, ACS Appl. Electron. Mater., 2022, 4, 3278–3302.
59. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864.
60. A. D. Becke, J. Chem. Phys., 2014, 140, 18A301.
61. L. H. Thomas, Mathematical proceedings of the Cambridge philosophical society, 1927, pp. 542–548.
62. E. Fermi, Rend. Accad. Naz. Lincei, 1927, 6, 5.
63. P. A. Dirac, Mathematical proceedings of the Cambridge philosophical society, 1930, pp. 376–385.
64. W. Yang, Phys. Rev. A, 1986, 34, 4575.
65. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864.
66. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
67. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
68. R. Van Noorden, B. Maher and R. Nuzzo, Nature, 2014, 514, 550–553.
69. J. Kubler, K. H. Hock, J. Sticht and A. R. Williams, J. Phys. F: Met. Phys., 1988, 18, 469.
70. R. Cuadrado and J. I. Cerdá, J. Phys.: Condens. Matter, 2012, 24, 086005.
71. F. Gómez-Ortiz, N. Carral-Sainz, J. Sifuna, V. Monteseguro, R. Cuadrado, P. García-Fernández and J. Junquera, Comput. Phys. Commun., 2023, 286, 108684.
72. C. L. Fu and K. M. Ho, Phys. Rev. B: Condens. Matter Mater. Phys., 1983, 28, 5480–5486.
73. M. Methfessel and A. T. Paxton, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 40, 3616–3621.
74. N. Marzari, D. Vanderbilt, A. De Vita and M. C. Payne, Phys. Rev. Lett., 1999, 82, 3296–3299.
75. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
76. A. Tkatchenko and M. Scheffler, Phys. Rev. Lett., 2009, 102, 073005.
77. H. Rydberg, M. Dion, N. Jacobson, E. Schröder, P. Hyldgaard, S. I. Simak, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2003, 91, 126402.
78. F. Nogueira, A. Castro and M. A. L. Marques, in A Tutorial on Density Functional Theory, ed. C. Fiolhais, F. Nogueira and M. A. L. Marques, Springer Berlin Heidelberg, Berlin, Heidelberg, 2003, pp. 218–256.
79. P. Mori-Sánchez, A. J. Cohen and W. Yang, J. Chem. Phys., 2006, 125, 201102.
80. P. Mori-Sánchez, A. J. Cohen and W. Yang, Phys. Rev. Lett., 2008, 100, 146401.
81. P. Mori-Sánchez, A. J. Cohen and W. Yang, Phys. Rev. Lett., 2009, 102, 066403.
82. P. Mori-Sánchez and A. J. Cohen, Phys. Chem. Chem. Phys., 2014, 16, 14378–14387.
83. V. I. Anisimov, J. Zaanen and O. K. Andersen, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 44, 943–954.
84. V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyżyk and G. A. Sawatzky, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 16929–16934.
85. I. V. Solovyev, P. H. Dederichs and V. I. Anisimov, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 16861–16871.
86. B. Himmetoglu, A. Floris, S. de Gironcoli and M. Cococcioni, Int. J. Quantum Chem., 2014, 114, 14–49.
87. V. L. Campo and M. Cococcioni, J. Phys.: Condens. Matter, 2010, 22, 055602.
88. J. P. Perdew, R. G. Parr, M. Levy and J. L. Balduz, Phys. Rev. Lett., 1982, 49, 1691–1694.
89. M. D. Towler, N. L. Allan, N. M. Harrison, V. R. Saunders, W. C. Mackrodt and E. Aprà, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 5041–5054.
90. O. Bengone, M. Alouani, P. Blöchl and J. Hugel, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 16392–16401.
91. I. I. Mazin and V. I. Anisimov, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 55, 12822–12825.
92. Z. Fang, K. Terakura, H. Sawada, T. Miyazaki and I. Solovyev, Phys. Rev. Lett., 1998, 81, 1027–1030.
93. Z. Fang, I. V. Solovyev, H. Sawada and K. Terakura, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 762–774.
94. A. I. Liechtenstein, V. I. Anisimov and J. Zaanen, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 52, R5467–R5470.
95. J. Hubbard and B. H. Flowers, Proc. R. Soc. London, Ser. A, 1963, 276, 238–257.
96. A. I. Liechtenstein, V. I. Anisimov and J. Zaanen, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 52, R5467–R5470.
97. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
98. F. Menescardi and D. Ceresoli, Appl. Sci., 2021, 11, 12.
99. A. N. Rudenko, M. Rösner and M. I. Katsnelson, npj Comput. Mater., 2023, 9, 83.
100. M. Cococcioni and S. de Gironcoli, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 035105.
101. N. J. Mosey and E. A. Carter, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 155123.
102. N. J. Mosey, P. Liao and E. A. Carter, J. Chem. Phys., 2008, 129, 014103.
103. I. Timrov, N. Marzari and M. Cococcioni, Phys. Rev. B, 2018, 98, 085127.
104. I. Timrov, N. Marzari and M. Cococcioni, Phys. Rev. B, 2021, 103, 045141.
105. I. Timrov, N. Marzari and M. Cococcioni, Comput. Phys. Commun., 2022, 279, 108455.
106. F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann and A. I. Lichtenstein, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 70, 195104.
107. F. Aryasetiawan, K. Karlsson, O. Jepsen and U. Schönberger, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 74, 125106.
108. B. Amadon, T. Applencourt and F. Bruneval, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 125110.
109. M. Yu, S. Yang, C. Wu and N. Marom, npj Comput. Mater., 2020, 6, 180.
110. W. Yu, Z. Zhang, X. Wan, H. Guo, Q. Gui, Y. Peng, Y. Li, W. Fu, D. Lu, Y. Ye and Y. Guo, J. Chem. Theory Comput., 2023, 19, 6425–6433.
111. B. Shahriari, K. Swersky, Z. Wang, R. P. Adams and N. de Freitas, Proc. IEEE, 2016, 104, 148–175.
112. X. Wang, Y. Jin, S. Schmitt and M. Olhofer, ACM Comput. Surv., 2023, 55, 1–36.
113. S. Baroni, S. de Gironcoli, A. Dal Corso and P. Giannozzi, Rev. Mod. Phys., 2001, 73, 515–562.
114. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
115. M. Cococcioni and N. Marzari, Phys. Rev. Mater., 2019, 3, 033801.
116. I. Timrov, F. Aquilante, M. Cococcioni and N. Marzari, PRX Energy, 2022, 1, 033003.
117. N. Naveas, R. Pulido, C. Marini, J. Hernández-Montelongo and M. M. Silván, iScience, 2023, 26, 106033.
118. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. D. Jr, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, J. Phys.: Condens. Matter, 2017, 29, 465901.
119. P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru and S. Baroni, J. Chem. Phys., 2020, 152, 154105.
120. G. Pizzi, A. Cepellotti, R. Sabatini, N. Marzari and B. Kozinsky, Comput. Mater. Sci., 2016, 111, 218–230.
121. S. P. Huber, S. Zoupanos, M. Uhrin, L. Talirz, L. Kahle, R. Häuselmann, D. Gresch, T. Müller, A. V. Yakutovich, C. W. Andersen, F. F. Ramirez, C. S. Adorf, F. Gargiulo, S. Kumbhar, E. Passaro, C. Johnston, A. Merkys, A. Cepellotti, N. Mounet, N. Marzari, B. Kozinsky and G. Pizzi, Sci. Data, 2020, 7, 300.
122. M. Uhrin, S. P. Huber, J. Yu, N. Marzari and G. Pizzi, Comput. Mater. Sci., 2021, 187, 110086.
123. E. Macke, I. Timrov, N. Marzari and L. C. Ciacchi, J. Chem. Theory Comput., 2024, 20, 4824–4843.
124. D. L. Esteras and J. J. Baldoví, Mater. Today Electron., 2023, 6, 100072.
125. D. L. Esteras, A. Rybakov, A. M. Ruiz and J. J. Baldoví, Nano Lett., 2022, 22, 8771–8778.
126. Y. Wang, N. Luo, J. Zeng, L.-M. Tang and K.-Q. Chen, Phys. Rev. B, 2023, 108, 054401.
127. P. Verma and D. G. Truhlar, Theor. Chem. Acc., 2016, 135, 182.
128. X.-B. Feng and N. M. Harrison, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 035114.
129. I. de P. R. Moreira, F. Illas and R. L. Martin, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 155102.
130. S. Grimme, A. Hansen, J. G. Brandenburg and C. Bannwarth, Chem. Rev., 2016, 116, 5105–5154.
131. M. Kim, W. J. Kim, E. K. Lee, S. Lebègue and H. Kim, Int. J. Quantum Chem., 2016, 116, 598–607.
132. L. I. Kushchuk, D. K. Veretimus, P. V. Lega, A. Y. Antonenkova and A. I. Kartsev, J. Surf. Invest.: X-Ray, Synchrotron Neutron Tech., 2024, 18, 859–863.
133. B. L. Chittari, Y. Park, D. Lee, M. Han, A. H. MacDonald, E. Hwang and J. Jung, Phys. Rev. B, 2016, 94, 184428.
134. R. Cuadrado, M. Pruneda, A. García and P. Ordejón, J. Phys.: Mater., 2018, 1, 015010.
135. X. Jiang, Q. Liu, J. Xing, N. Liu, Y. Guo, Z. Liu and J. Zhao, Appl. Phys. Rev., 2021, 8, 031305.
136. F. Zheng and P. Zhang, Comput. Phys. Commun., 2021, 259, 107659.
137. M. K. Horton, J. H. Montoya, M. Liu and K. A. Persson, npj Comput. Mater., 2019, 5, 64.
138. L. M. Sandratskii, Phys. Status Solidi B, 1986, 136, 167–180.
139. L. M. Sandratskii, J. Phys.: Condens. Matter, 1991, 3, 8565.
140. P. Kurz, F. Förster, L. Nordström, G. Bihlmayer and S. Blügel, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 024415.
141. G. Bihlmayer, in Density-functional Theory of Magnetism, John Wiley & Sons, Ltd, 2007.
142. T. Olsen, Phys. Rev. B, 2016, 94, 235106.
143. L. M. Sandratskii, Phys. Rev. B, 2017, 96, 024450.
144. B. Zimmermann, G. Bihlmayer, M. Böttcher, M. Bouhassoune, S. Lounis, J. Sinova, S. Heinze, S. Blügel and B. Dupé, Phys. Rev. B, 2019, 99, 214426.
145. M. Dupé, S. Heinze, J. Sinova and B. Dupé, Phys. Rev. B, 2018, 98, 224415.
146. J. J. Mortensen, L. B. Hansen and K. W. Jacobsen, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 035109.
147. J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen, M. Dułak, L. Ferrighi, J. Gavnholt, C. Glinsvad, V. Haikola, H. A. Hansen, H. H. Kristoffersen, M. Kuisma, A. H. Larsen, L. Lehtovaara, M. Ljungberg, O. Lopez-Acevedo, P. G. Moses, J. Ojanen, T. Olsen, V. Petzold, N. A. Romero, J. Stausholm-Møller, M. Strange, G. A. Tritsaris, M. Vanin, M. Walter, B. Hammer, H. Häkkinen, G. K. H. Madsen, R. M. Nieminen, J. K. Nørskov, M. Puska, T. T. Rantala, J. Schiøtz, K. S. Thygesen and K. W. Jacobsen, J. Phys.: Condens. Matter, 2010, 22, 253202.
148. J. J. Mortensen, A. H. Larsen, M. Kuisma, A. V. Ivanov, A. Taghizadeh, A. Peterson, A. Haldar, A. O. Dohn, C. Schäfer, E. Ö. Jónsson, E. D. Hermes, F. A. Nilsson, G. Kastlunger, G. Levi, H. Jónsson, H. Häkkinen, J. Fojt, J. Kangsabanik, J. Sødequist, J. Lehtomäki, J. Heske, J. Enkovaara, K. T. Winther, M. Dulak, M. M. Melander, M. Ovesen, M. Louhivuori, M. Walter, M. Gjerding, O. Lopez-Acevedo, P. Erhart, R. Warmbier, R. Würdemann, S. Kaappa, S. Latini, T. M. Boland, T. Bligaard, T. Skovhus, T. Susi, T. Maxson, T. Rossi, X. Chen, Y. L. A. Schmerwitz, J. Schiøtz, T. Olsen, K. W. Jacobsen and K. S. Thygesen, J. Chem. Phys., 2024, 160, 092503.
149. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys.: Condens. Matter, 2017, 29, 273002.
150. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
151. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
152. The FLEUR project, https://www.flapw.de/.
153. D. Wortmann, G. Michalicek, N. Baadji, M. Betzinger, G. Bihlmayer, J. Bröder, T. Burnus, J. Enkovaara, F. Freimuth, C. Friedrich, C.-R. Gerhorst, S. Granberg Cauchi, U. Grytsiuk, A. Hanke, J.-P. Hanke, M. Heide, S. Heinze, R. Hilgers, H. Janssen, D. A. Klüppelberg, R. Kovacik, P. Kurz, M. Lezaic, G. K. H. Madsen, Y. Mokrousov, A. Neukirchen, M. Redies, S. Rost, M. Schlipf, A. Schindlmayr, M. Winkelmann and S. Blügel, FLEUR, Zenodo, 2023,  DOI:10.5281/zenodo.7576163.
154. T. Ozaki, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 67, 155108.
155. T. Ozaki and H. Kino, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 195113.
156. T. Ozaki and H. Kino, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 72, 045121.
157. T. B. Prayitno and F. Ishii, J. Phys. Soc. Jpn., 2018, 87, 114709.
158. T. B. Prayitno and F. Ishii, J. Phys. Soc. Jpn., 2019, 88, 054701.
159. Y. Zhu, X. Kong, T. D. Rhone and H. Guo, Phys. Rev. Mater., 2018, 2, 081001.
160. K. Yang, G. Wang, L. Liu, D. Lu and H. Wu, Phys. Rev. B, 2021, 104, 144416.
161. A. M. Ruiz, D. L. Esteras, A. Rybakov and J. J. Baldoví, Dalton Trans., 2022, 51, 16816–16823.
162. K. Lee, A. H. Dismukes, E. J. Telford, R. A. Wiscons, J. Wang, X. Xu, C. Nuckolls, C. R. Dean, X. Roy and X. Zhu, Nano Lett., 2021, 21, 3511–3517.
163. K. Lee, A. H. Dismukes, E. J. Telford, R. A. Wiscons, J. Wang, X. Xu, C. Nuckolls, C. R. Dean, X. Roy and X. Zhu, Nano Lett., 2021, 21, 3511–3517.
164. A. M. Ruiz, D. L. Esteras, D. López-Alcalá and J. J. Baldoví, Nano Lett., 2024, 24, 7886–7894.
165. X. Li, M. Zhu, Y. Wang, F. Zheng, J. Dong, Y. Zhou, L. You and J. Zhang, Appl. Phys. Lett., 2023, 122, 082404.
166. P.-W. Ma and S. L. Dudarev, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 054420.
167. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
168. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
169. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
170. E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer, S. Spicher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2019, 150, 154122.
171. J. Hermann, M. Stöhr, S. Góger, S. Chaudhuri, B. Aradi, R. J. Maurer and A. Tkatchenko, J. Chem. Phys., 2023, 159, 174802.
172. A. Tkatchenko and M. Scheffler, Phys. Rev. Lett., 2009, 102, 073005.
173. T. Bučko, S. Lebègue, J. Hafner and J. G. Ángyán, J. Chem. Theory Comput., 2013, 9, 4293–4299.
174. T. Bučko, S. Lebègue, J. G. Ángyán and J. Hafner, J. Chem. Phys., 2014, 141, 034114.
175. A. Tkatchenko, R. A. DiStasio, R. Car and M. Scheffler, Phys. Rev. Lett., 2012, 108, 236402.
176. A. Ambrosetti, A. M. Reilly, J. DiStasio, A. Robert and A. Tkatchenko, J. Chem. Phys., 2014, 140, 18A508.
177. T. Gould and T. Bučko, J. Chem. Theory Comput., 2016, 12, 3603–3613.
178. T. Gould, S. Lebègue, J. G. Ángyán and T. Bučko, J. Chem. Theory Comput., 2016, 12, 5920–5930.
179. J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Junquera, P. Ordejon and D. Sanchez-Portal, J. Phys.: Condens. Matter, 2002, 14, 2745.
180. M. Dion, H. Rydberg, E. Schröder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401.
181. K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist and D. C. Langreth, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 081101.
182. J. Klimeš, D. R. Bowler and A. Michaelides, J. Phys.: Condens. Matter, 2009, 22, 022201.
183. V. R. Cooper, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 161104.
184. K. Berland and P. Hyldgaard, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 035412.
185. O. A. Vydrov and T. Van Voorhis, J. Chem. Phys., 2010, 133, 244103.
186. M. Heide, G. Bihlmayer and S. Blügel, Phys. B, 2009, 404, 2678–2683.
187. A. B. Shick, A. I. Liechtenstein and W. E. Pickett, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 10763–10769.
188. V. Barone, M. Casarin, D. Forrer, M. Pavone, M. Sambi and A. Vittadini, J. Comput. Chem., 2009, 30, 934–939.
189. A. D. Becke and E. R. Johnson, J. Chem. Phys., 2007, 127, 154108.
190. A. Otero-de-la Roza and E. R. Johnson, J. Chem. Phys., 2012, 136, 174109.
191. M. J. Han, T. Ozaki and J. Yu, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 045110.
192. S. Ryee and M. J. Han, J. Phys.: Condens. Matter, 2018, 30, 275802.
193. M. Dion, H. Rydberg, E. Schröder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2005, 95, 109902.
194. X. Gonze, B. Amadon, G. Antonius, F. Arnardi, L. Baguet, J.-M. Beuken, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, N. Brouwer, F. Bruneval, G. Brunin, T. Cavignac, J.-B. Charraud, W. Chen, M. Côté, S. Cottenier, J. Denier, G. Geneste, P. Ghosez, M. Giantomassi, Y. Gillet, O. Gingras, D. R. Hamann, G. Hautier, X. He, N. Helbig, N. Holzwarth, Y. Jia, F. Jollet, W. Lafargue-Dit-Hauret, K. Lejaeghere, M. A. L. Marques, A. Martin, C. Martins, H. P. C. Miranda, F. Naccarato, K. Persson, G. Petretto, V. Planes, Y. Pouillon, S. Prokhorenko, F. Ricci, G.-M. Rignanese, A. H. Romero, M. M. Schmitt, M. Torrent, M. J. van Setten, B. V. Troeye, M. J. Verstraete, G. Zérah and J. W. Zwanziger, Comput. Phys. Commun., 2020, 248, 107042.
195. X. Gonze, B. Seddon, J. A. Elliott, C. Tantardini and A. V. Shapeev, J. Chem. Theory Comput., 2022, 18, 6099–6110.
196. B. Van Troeye, M. Torrent and X. Gonze, Phys. Rev. B, 2016, 93, 144304.
197. B. Amadon, F. Jollet and M. Torrent, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 155104.
198. S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. I. J. Probert, K. Refson and M. C. Payne, Z. Kristallogr. - Cryst. Mater., 2005, 220, 567–570.
199. F. Ortmann, F. Bechstedt and W. G. Schmidt, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 205101.
200. A. D. Becke and E. R. Johnson, J. Chem. Phys., 2007, 127, 124108.
201. E. D. Team, ELK: Embedded Learning Library, 2024, https://elk.sourceforge.io/, accessed: 2024-10-29.
202. P. Blaha, K. Schwarz, F. Tran, R. Laskowski, G. K. H. Madsen and L. D. Marks, J. Chem. Phys., 2020, 152, 074101.
203. H. Peng, Z.-H. Yang, J. P. Perdew and J. Sun, Phys. Rev. X, 2016, 6, 041005.
204. X. Lu, R. Fei and L. Yang, Phys. Rev. B, 2019, 100, 205409.
205. D. A. Wahab, M. Augustin, S. M. Valero, W. Kuang, S. Jenkins, E. Coronado, I. V. Grigorieva, I. J. Vera-Marun, E. Navarro-Moratalla, F. L. Evans, K. S. Novoselov and E. J. G. Santos, Adv. Mater., 2020, 33, 29.
206. T. Olsen, Phys. Rev. B, 2017, 96, 125143.
207. D. Ködderitzsch, W. Hergert, W. M. Temmerman, Z. Szotek, A. Ernst and H. Winter, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 66, 064434.
208. D. Torelli and T. Olsen, J. Phys.: Condens. Matter, 2020, 32, 335802.
209. B. Xu, S. Meyer, M. J. Verstraete, L. Bellaiche and B. Dupé, Phys. Rev. B, 2021, 103, 214423.
210. I. V. Solovyev, J. Phys.: Condens. Matter, 2024, 36, 223001.
211. A. Szilva, Y. Kvashnin, E. A. Stepanov, L. Nordström, O. Eriksson, A. I. Lichtenstein and M. I. Katsnelson, Rev. Mod. Phys., 2023, 95, 035004.
212. M. I. Katsnelson and A. I. Lichtenstein, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 61, 8906–8912.
213. V. V. Mazurenko and V. I. Anisimov, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 184434.
214. X. He, N. Helbig, M. J. Verstraete and E. Bousquet, Comput. Phys. Commun., 2021, 264, 107938.
215. D. M. Korotin, V. V. Mazurenko, V. I. Anisimov and S. V. Streltsov, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 224405.
216. F. L. Durhuus, T. Skovhus and T. Olsen, J. Phys.: Condens. Matter, 2023, 35, 105802.
217. V. Antropov, M. Katsnelson and A. Liechtenstein, Phys. B, 1997, 237–238, 336–340.
218. M. I. Katsnelson, Y. O. Kvashnin, V. V. Mazurenko and A. I. Lichtenstein, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 100403.
219. A. Szilva, M. Costa, A. Bergman, L. Szunyogh, L. Nordström and O. Eriksson, Phys. Rev. Lett., 2013, 111, 127204.
220. G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi and J. R. Yates, J. Phys.: Condens. Matter, 2020, 32, 165902.
221. H. Yoon, T. J. Kim, J.-H. Sim and M. J. Han, Comput. Phys. Commun., 2020, 247, 106927.
222. L. Oroszlány, J. Ferrer, A. Deák, L. Udvardi and L. Szunyogh, Phys. Rev. B, 2019, 99, 224412.
223. G. Martínez-Carracedo, L. Oroszlány, A. García-Fuente, B. Nyári, L. Udvardi, L. Szunyogh and J. Ferrer, Phys. Rev. B, 2023, 108, 214418.
224. X. Lu, R. Fei and L. Yang, Phys. Rev. B, 2019, 100, 205409.
225. N. Majlis, The Quantum Theory of Magnetism, World Scientific, 2nd edn, 2007.
226. M. Pajda, J. Kudrnovský, I. Turek, V. Drchal and P. Bruno, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 174402.
227. K. Malarz, M. Zborek and B. Wrobel, Curie temperature for an Ising model on Archimedean lattices, 2005.
228. A. Codello, J. Phys. A: Math. Theor., 2010, 43, 385002.
229. J. Vanherck, C. Bacaksiz, B. Sorée, M. V. Miloševič and W. Magnus, Appl. Phys. Lett., 2020, 117, 052401.
230. Y. Zhu, X. Kong, T. D. Rhone and H. Guo, Phys. Rev. Mater., 2018, 2, 081001.
231. T. Holstein and H. Primakoff, Phys. Rev., 1940, 58, 1098–1113.
232. V. Rajeev Pavizhakumari, T. Skovhus and T. Olsen, J. Phys.:Condens. Matter, 2025, 37, 115806.
233. F. J. Dyson, Phys. Rev., 1956, 102, 1217–1230.
234. S. V. Maleev, Sov. Phys. JETP, 1958, 6, 776.
235. D. Stanek, O. P. Sushkov and G. S. Uhrig, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 064505.
236. C. Jordan, Calculus of Finite Differences, American Elsevier Publishing Company, Inc., New York, 2nd edn, 1965.
237. J. König and A. Hucht, SciPost Phys., 2021, 10, 007.
238. M. Vogl, P. Laurell, H. Zhang, S. Okamoto and G. A. Fiete, Phys. Rev. Res., 2020, 2, 043243.
239. V. V. Mkhitaryan and L. Ke, Phys. Rev. B, 2021, 104, 064435.
240. W. P. Lima, F. R. V. Araújo, D. R. da Costa, S. H. R. Sena and J. M. Pereira, Braz. J. Phys., 2022, 52, 42.
241. T. Olsen, MRS Commun., 2019, 9, 1142–1150.
242. W. Krauth, Statistical Mechanics Algorithms and Computations, Oxford University Press, 1st edn, 2006.
243. D. Hinzke and U. Nowak, Comput. Phys. Commun., 1999, 121–122, 334–337.
244. J. Liu, Q. Sun, Y. Kawazoe and P. Jena, Phys. Chem. Chem. Phys., 2016, 18, 8777–8784.
245. H. H. Kim, B. Yang, S. Li, S. Jiang, C. Jin, Z. Tao, G. Nichols, F. Sfigakis, S. Zhong, C. Li, S. Tian, D. G. Cory, G.-X. Miao, J. Shan, K. F. Mak, H. Lei, K. Sun, L. Zhao and A. W. Tsen, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 11131–11136.
246. Z. Zhang, J. Shang, C. Jiang, A. Rasmita, W. Gao and T. Yu, Nano Lett., 2019, 19, 3138–3142.
247. J. D. Alzate-Cardona, D. Sabogal-Suárez, R. F. L. Evans and E. Restrepo-Parra, J. Phys.: Condens. Matter, 2019, 31, 095802.
248. Y. Zhang, B. Wang, Y. Guo, Q. Li and J. Wang, Comput. Mater. Sci., 2021, 197, 110638.
249. R. F. L. Evans, U. Atxitia and R. W. Chantrell, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 144425.
250. S. Jenkins, L. Rózsa, U. Atxitia, R. F. L. Evans, K. S. Novoselov and E. J. G. Santos, Nat. Commun., 2022, 13, 6917.
251. P. Asselin, R. F. L. Evans, J. Barker, R. W. Chantrell, R. Yanes, O. Chubykalo-Fesenko, D. Hinzke and U. Nowak, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 054415.
252. VAMPIRE software package, version 6, https://vampire.york.ac.uk/.
253. G. P. Müller, M. Hoffmann, C. Dißelkamp, D. Schürhoff, S. Mavros, M. Sallermann, N. S. Kiselev, H. Jónsson and S. Blügel, Phys. Rev. B, 2019, 99, 224414.
254. P. Rüßmann, J. Ribas Sobreviela, M. Sallermann, M. Hoffmann, F. Rhiem and S. Blügel, Front. Mater., 2022, 9, 1–14.
255. M.-A. Bisotti, D. Cortés-Ortuño, R. Pepper, W. Wang, M. Beg, T. Kluyver and H. Fangohr, J. Open Res. Softw., 2018, 6, 22.
256. B. Skubic, J. Hellsvik, L. Nordström and O. Eriksson, J. Phys.: Condens. Matter, 2008, 20, 315203.
257. W. T. Coffey and Y. P. Kalmykov, The Langevin Equation, World Scientific, 3rd edn, 2012.
258. H. B. Callen, Phys. Rev., 1963, 130, 890–898.
259. J. Vanherck, B. Sorée and W. Magnus, J. Phys.: Condens. Matter, 2018, 30, 275801.
260. A. Pelissetto and E. Vicari, Phys. Rep., 2002, 368, 549–727.
261. S. R. A. Salinas, Introduction to statistical physics, Springer New York, New York, 2001.
262. M. Bander and D. L. Mills, Phys. Rev. B, 1988, 38, 12015–12018.
263. X. Wang, D. Li, Z. Li, C. Wu, C.-M. Che, G. Chen and X. Cui, ACS Nano, 2021, 15, 16236–16241.
264. H.-R. Fuh, C.-R. Chang, Y.-K. Wang, R. F. L. Evans, R. W. Chantrell and H.-T. Jeng, Sci. Rep., 2016, 6, 32625.
265. H. Sheng, H. Long, G. Zou, D. Bai, J. Zhang and J. Wang, J. Mater. Sci., 2021, 56, 15844–15858.
266. M. Abdollahi and M. B. Tagani, Phys. Rev. B, 2023, 108, 024427.
267. Y. Guo, Y. Zhang, S. Yuan, B. Wang and J. Wang, Nanoscale, 2018, 10, 18036–18042.
268. H. Wang, J. Qi and X. Qian, Appl. Phys. Lett., 2020, 117, 083102.
269. H. H. Kim, B. Yang, S. Li, S. Jiang, C. Jin, Z. Tao, G. Nichols, F. Sfigakis, S. Zhong, C. Li, S. Tian, D. G. Cory, G.-X. Miao, J. Shan, K. F. Mak, H. Lei, K. Sun, L. Zhao and A. W. Tsen, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 11131–11136.
270. Z. Zhang, J. Shang, C. Jiang, A. Rasmita, W. Gao and T. Yu, Nano Lett., 2019, 19, 3138–3142.
271. W. Chen, J. M. Zhang, Y. Z. Nie, Q. L. Xia and G. H. Guo, J. Magn. Magn. Mater., 2020, 508, 166878.
272. A. Bedoya-Pinto, J.-R. Ji, A. K. Pandeya, P. Gargiani, M. Valvidares, P. Sessi, J. M. Taylor, F. Radu, K. Chang and S. S. P. Parkin, Science, 2021, 374, 616–620.
273. M. Wang, B. Lei, K. Zhu, Y. Deng, M. Tian, Z. Xiang, T. Wu and X. Chen, npj 2D Mater. Appl., 2024, 8, 22.
274. G. Zhang, F. Guo, H. Wu, X. Wen, L. Yang, W. Jin, W. Zhang and H. Chang, Nat. Commun., 2022, 13, 5067.
275. Z.-X. Shen, X. Bo, K. Cao, X. Wan and L. He, Phys. Rev. B, 2021, 103, 085102.
276. X. Chen, Z.-Z. Lin and L.-R. Cheng, Chin. Phys. B, 2021, 30, 047502.

---

Footnote

† And mainly those of the current approximate exchange correlation functionals.

---