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?
Any suggestion will be helpful.
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."
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.
\begin{equation} G\: =\: H_{T=0}\: +\: H_{T>0}\: +\: \text{zero point energy}\: -\: T\cdot S, \tag{1} \end{equation} 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.
The previous answers are good responses, especially the one from Guy Moore. However, I'd like to recommend two more solutions: curie_calculator and mcsolver.
A user friendly and efficient tool implementing Monte Carlo simulations to estimate Curie/Neel temperature. Supports multiple occasions, e.g. standard ferromagnetic/anti-ferromagnetic systems, DMI, Kiteav non-diagonal exchange interactions, dipole-dipole long-range couplings, with external fields.
You can either choose the Windows version or python package. Both of them need same input file.
Define the three lattice vectors of primitive cell.
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.
Define all exchange interactions (bonds). There are nine matrix elements for one $J$ (i.e. $J_{xx}$, $J_{yy}$, $J_{zz}$, $J_{xy}$, etc.). Each element describes the coupling between two components of spins. For example, a basic bond term can be expressed as $$S_1\cdot J\cdot S_2 = 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 the Ising model, since only one component is available for each spin, only the first element $J_{xx}$ is used. As well as for the XY model, only $J_{xx}$, $J_{yy}$, $J_{xy}$, $J_{yz}$ are used. In this step, you can click one of the bonds to review the actual linking in lattice on the view panel. Activated bond is depicted with bold and yellow line while others are green. You may drag the left/right mouse button to rotate and expand/shrink the model shown in the view panel. 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 Zhang et. al.$^1$ 10.1021/acs.jpclett.0c01911.
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 statesnsweep is the total steps involved in measuringtau denotes the MC updates for each stepxAxis denotes the physical quantity put in x-axis of right-hand Result viewer, it can be either T(for illustration of M-T curve) or H (for illustration of hysteresis loop) model type, algorithm (only Metropolis and Wolff are supported now).nFrame is the number of output spin configurations, using for illustrating spin configurations in equilibrium or non-equilibrium states.Set spin_i and spin_j and the lattice vector between them, for correlation measurements. (if spin_i=spin_j and overLat=0 0 0, then you will get susceptibility for spin_i)
Set the core resources for parallel calculation.
(Optional) Save current parameters into file.
Click the startMC button to start.
Wait for the diagram update in the right panel. Afterward, you can find a file result.txt in the root directory of mcsolver, there are many useful information including the averaged spin (on spin_i and spin_j defined in step 5), correlation between spin_i and spin_j, internal energy, specific heat capacity, 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.
A detailed tutorial of mcsolver have been given in its README file. To repeat the $\ce{CrI3}$ case in samples directory is a good choice.
