# Is it possible to calculate the Curie temperature for magnetic systems?

Using DFT calculations we can know the electronic and magnetic property of the system at 0 K. I am wondering if there is any computational method to know the curie temperature of magnetic materials? Can I know the Curie temperature using Quantum Espresso?

Thank you!

This is an important question, and I will do my best to answer it in a way that adds to the nice suggestions made in earlier answers.

In short, the answer is yes, but only with the help of some other models and approaches. In this response, I will focus on the conventional (computational) approach to this question, using lattice models and Monte Carlo Methods.

This approach, generally speaking, involves running a series of spin-polarized DFT calculations for different magnetic configurations. The next step is to map the results of these calculations (total energies and/or electronic structure) to a lattice model. There a two main ways to do this, which I touch on below. The transition temperature (Curie temperature for ferromagnets) can then be identified by running Monte Carlo (MC) simulations on this lattice model, and observing at what temperature your order parameter (in this case total magnetization) goes to zero (paramagnetic transition), or heat capacity "spikes," etc... The VAMPIRE documentation for this is included here:

https://vampire.york.ac.uk/tutorials/curie-temperature-simulation/

In theory, you could perform MC runs with DFT energy evaluations. However, practically speaking the lattice model provides a cheap alternative to exploring the magnetic configuration space of your system.

I encourage you explore the magnetic exchange workflow recently implemented in atomate:

https://atomate.org/atomate.vasp.workflows.base.html#module-atomate.vasp.workflows.base.exchange

which calculates the parameters of a lattice model from DFT energies (using pymatgen functions linked below), and estimates the transition temperature using the MC method implemented in the VAMPIRE code suite (suggested in a previous answer).

The workflow calls the Heisenberg model mapper implemented in pymatgen: https://pymatgen.org/pymatgen.analysis.magnetism.heisenberg.html

Some more details

In the lattice picture, the energy of arrangements of different spins can be calculated in a relatively straightforward way using the effective Hamiltonian:

$$\mathcal{H} = - \sum_{\langle i,j \rangle} J_{ij} \delta_{ij}$$

where the sum is over all pairwise interactions between sites (the pair of sites $$i$$ and $$j$$ is indicated by $$\langle i,j \rangle$$). In the Ising picture, $$\delta_{ij} = \sigma_i \sigma_j$$ (where $$\sigma_k = \pm 1$$ indicates the up/down spin state at site $$k$$). In the classical Heisenberg picture, $$\delta_{ij} = \vec{s}_i \cdot \vec{s}_j$$ (the dot product between spins represented as vectors in 3D space).

One way to calculate $$J_{ij}$$ values from DFT using the KKR-Green's function method, another is by solving a matrix vector system (implemented in the atomate workflow). For a physical intuition behind these exchange interaction values, I suggest referring to a book in solid state physics/chemistry, such as the well known text by Ashcroft and Mermin.

Famous examples of these lattice models include Ising and Heisenberg (both "classical" and "quantum" formulations). However, there are many more, such as the Potts and XY models. Theorists have studied these models extensively for years. Despite the simplicity of these models compared to density functional theory approaches, these lattice models are able to capture the physics relevant to the order-to-disorder (ordered to paramagnetic) phase transition in magnets.

Common assumptions underlying these models include the approximation that each atom has an effective "spin" or magnetic moment on each atom. This assumption is reasonable for many systems with "localized" moments - more on some of the major limitations of this simplification here. Alternatively, there are other approaches that include "higher order" expansions to the free energy of your system using "cluster expansion" or "cluster multipole expansion," however, I will not focus on those approaches here.

I am currently working on extending this mapping to include the effect of atomic displacements in the model, in addition to magnetic interactions. These spin-lattice (magnon-phonon) coupling effects have been shown to also shift the phase transition temperature. I've run out of allowed references, however, if your curious, search for "compressible Heisenberg/Ising," "spin-lattice coupling Hamiltonian," or "Bean-Rodbell" model.

A helpful lesson that I've learned is that it's important to be aware of the assumptions of different lattice models, and understanding the relevant physics that we want to capture. As the famous saying goes:

"All models are wrong but some are useful."

• This is a fantastic answer! – Matt Horton Dec 11 '20 at 20:15
• @Guy Moore thanks for the answer. It is really informative and exact. – Ujjawal M. Dec 22 '20 at 7:05

I am wondering if there is any computational method to know the curie temperature of magnetic materials?

Yes, there is!. Simulations to determine Curie Temperature are usually done via Monte-Carlo methods, with the assumptions taken from the Heisenberg model. I haven't done any Curie Temperature Simulations using DFT but it is possible as this answer on ResearchGate suggests.

Can I know the Curie temperature using Quantum Espresso?

Rather than using Quantum ESPRESSO or any such DFT based codes, you're better off using Atomistic Spin model based codes like VAMPIRE

There is even a tutorial on Curie Temperature simulations as given in their website

VAMPIRE does, in fact, have issues when considering unit cells of different symmetry in one simulation if you're creating something similar to a bi-layer or a core-shell geometry.

You could check out a Curie Temperature calculation I've done for a core-shell system here

(I'll answer with crystalline solids in mind.)

You could study the material's magnetic and non-magnetic phases and construct their Gibbs free energy (G) vs temperature (T) profiles.

$$$$G\: =\: H_{T=0}\: +\: H_{T>0}\: +\: \text{zero point energy}\: -\: T\cdot S, \tag{1}$$$$ Where H is enthalpy, and S is entropy (any kind).

Based on your system, you may have to alter equation $$1$$ and add/remove terms.

The following are certain ways to calculate each of the terms in equation $$1$$:

$$H_{T=0}$$: Density functional theory based geometry optimization calculations are my best guess. Luckily, you wanted this option as well.

$$H_{T>0}$$ and zero point energy: Lattice dynamics studies (vibrational frequency calculations, for molecules) should be sufficient. Codes like Phonopy, PHONON, and CRYSTAL, all work with DFT codes as force calculators. There is more to lattice dynamics than meets the eye: Quasiharmonic approximation and anharmonicity.

Finally, S: this is tricky. You need to figure out which kinds of entropy matter for your system in the temperature range of interest. Accordingly, you need to then calculate each of those entropies.

$$2^{nd}$$ method: I'm not too familiar with the CALPHAD (Calculation of phase diagram) method. As far as I understand, you create Gibbs free energy models from experimental or DFT data (or their hybrid) and compute G at various T. I just happened to ask this very question two days ago here: What information about a crystalline solid material goes into calculating phase diagrams in Thermo-Calc, and how?

For fluids, the $$P\cdot V$$ term may also have to be added on a case by case basis.

• CALPHAD method may not always results the correct curie temperature. It depends on the model we use to evaluate it ( i.e. Inden model, Inden-Jarl-Hillert model). This magnetic effect is implemented in the heat capacity term and if heat capacity at constant pressure (Cp) is plotted with respect to temperature, then it will be discontinuous at the curie temperature. – Niraja moharana Dec 14 '20 at 17:43
• @Nirajamoharana what are these models of? – Hitanshu Sachania Dec 14 '20 at 18:25
• Usually we express heat capacity (Cp) using Meyer-Kelly polynomial. But this model is not suitable to take care the magnetic effect. So there is a need of separate polynomial to implement the magnetic effect for Cp. Models by Inden, Jarl-Hillert are some of them. These models consist of material specific constants, a parameter called tau = T/(critical temperature for magnetic transition) and magnetic moment. – Niraja moharana Dec 15 '20 at 8:04
• Yes, you can call it as an iterative process. I didn't get that from your previous comment. Sorry for the delayed reply. @Hitanshu – Niraja moharana Dec 31 '20 at 17:00
• Nice to see a lot of discussion here. @Nirajamoharana have you seen the chat rooms: chat.stackexchange.com/… ? If this or any other discussion in comments gets too long, rather than starting a new room, you can use one of the existing rooms (I'm just letting people here know!). – Nike Dattani Apr 6 at 4:49

[Two codes are suggested]

The previous answers have given well reply, especially the one from Guy Moore. I wonder to recommend two more code.

First one: curie_calculator

https://github.com/Volkov-da/curie_calculator

Second one: mcsolver mcsolver: a user friendly and efficient tool implementing Monte Carlo simulations to estimate Curie/Neel temperature

Support multiple ocassions, e.g. standard ferromangetic/anti-ferromagnetic systems, DMI, Kiteav non-diagonal exchange interactions, dipole-dipole long-range couplings, with external fields.

https://github.com/golddoushi/mcsolver

You can either choose Windows version or python package.
Both of them need same input file.

How to define parameters?

1. Define the three lattice vectors of primitive cell.

2. Define all basic spins in primitive cell, note that the fractional coordinates are supposed. Ani represents the single-ion anisotropies in xyz directions (It is useless in Ising model, and only former two are used in XY model). As well as, note that the units of anisotropies are in Kelvin.

3. Define all exchange interactions (bonds). There are nine matrix elements for one J including Jxx, Jyy, Jzz, Jxy, Jxz, Jyz, Jyx, Jzx, Jzy, respectively. Each element discribes the coupling between two components of spins. For example, a basic bond term can be expressed as
$$S1\dot J\dot S2 = S_{1x}J_{xx}S_{2x} + S_{1y}J_{yy}S_{2y} + S_{1z}J_{zz}S_{2z} + S_{1x}J_{xy}S_{2y} + ...$$
For Ising model, since only one component is available for each spin, only the first element Jxx is used. As well as for XY model, only Jxx, Jyy, Jxy, Jyz are used.
In this step, you can click one of the bonds to review the actual linking in lattice on view pannel. Activated bond is depicted with bold and yellow line while others are green. You may drag left/right mouse Btn to rotate and expand/shrink the model shown in view pannel.
By the way, these parameters can be obtained via fitting experimental spin-wave spectra or DFT calculations. One of the general method for calculating anisotropic bonds and orbitals was introduced in dio: 10.1021/acs.jpclett.0c01911.

4. Define other parameters, including
the start and end temperatures, number of temperature interpolations (for the temperature scanning)
the start, end and number of samplings for external field (for the magnetic field scanning)
nthermal is the total steps to make system enter balanced states, nsweep is the total steps involved in mesearing, tau denotes the MC updates for each step
xAxis denote the physical quantity put in x-axis of right-hand Result viewer, it can be either T(for illustration of M-T curv) or H (for illustration of hysteresis loop).
model type, algorithm (only Metropolis and Wolff are supported now)
nFrame is the num. of output spin configurations, using for illustrating spin configurations in equilibrium or non-equilibrium states.

5. Set spin_i and spin_j and the lattice vector between them, for correlation mesearments. (if spin_i=spin_j and overLat=0 0 0, then you will get susceptibility for spin_i)

6. Set the core resources for parallel calc.

7. (Optional) Save current parameters into file.

8. Click startMC Btn to start.

9. Wait for the diagram update in right pannel. Afterwards, you can find a file result.txt in the root directory of mcsolver, there are many useful informations including the averaged spin (on spin_i and j defined in step 5), correlation between spin_i and j, internal energy, specific heat capacith, and Binder cumulant U4, etc. If you handle the sims with more than one cores then the results may not be ordered according to temperature, however, the correspondences in every line are ok.


Example: A details tutorial of mcsolver have been given in its README file. To repeat the CrI3 case in samples directory is a good choice.