Thermal Corrections

Thermal Corrections

Overview

With DFT or an interatomic potential, the computed energies are 0 K energies and therefore are not Gibbs free energies. Many times, what we want are Gibbs free energies at a given temperature and pressure. This document goes over how to do this.

Ideal Gas Thermochemistry

If you are calculating the thermochemistry of an ideal gas, you're in luck — this is straightforward. Most DFT codes will automatically compute the thermochemical corrections assuming an ideal gas model when you carry out a frequency calculation. If you need to do this yourself (e.g. for use with an MLIP), ASE has an IdealGasThermo  class  that will get the job done. All you need is the optimized structure and the corresponding vibrational energies (i.e. a unit conversion from the frequencies). The theory behind this is described in  ⚙️⁠⁠Ideal Gas Thermochemistry⁠ . There is also the  GoodVibes  code that includes some fancy features but is not yet compatible with MLIPs.

Solid-State Thermochemistry

Calculating the thermochemistry of a solid is much trickier. If you are doing such calculations, first read  this review  and the references therein as needed.

Quasi-Harmonic Approximation (QHA)

If the solid has only vibrational degrees of freedom, the most routine approach is to use the quasi-harmonic approximation (QHA) as outlined in the  phonopy code . QHA is not just for Gibbs free energies — it will also produce thermal expansion coefficients and heat capacities for you.
An example of how to do a QHA calculation using the QHACalc class in  matcalc  is included below. In practice, this is too expensive to do with DFT in most cases, and MLIPs can prove critical.
Until matcalc >0.4.5 is released, you must use the main branch of matcalc via pip install git+https://github.com/materialyzeai/matcalc.git
import logging
from matcalc import QHACalc
import matplotlib.pyplot as plt
from ase.io import read

# Set logging level
logging.basicConfig(level=logging.INFO)

# Read Atoms object
atoms = read("...")

# Instantiate calculator as usual
calc = ...

# Run QHA from 0 K to 600.0 K at 1 bar
results = QHACalc(calc, pressure=1e-4, t_max=600.0, fmax=1e-6, optimizer="FIRE", max_steps=10000, scale_factors=(0.97, 0.98, 0.99, 1.0, 1.01, 1.02, 1.03), fix_imaginary_attempts=1).calc(atoms)

# Fetch some results
thermal_expansion = results["thermal_expansion_coefficients"]
c_P = results["heat_capacity_P"]
results["qha"].plot_qha()
plt.savefig("qha.png")
There are some important features in the above code:
    The pressure is in units of GPa.
    The structure must be converged to a very tight force tolerance before doing a phonon calculation, to prevent spurious imaginary modes.
    The FIRE optimizer should be used with QHACalc if you are allowing the unit cell shape to change, which is the default.
    The QHA calculation is physically meaningless if there are imaginary modes. There is a fix_imaginary_attempts keyword argument that will automatically rattle the atoms and rerun the relaxation and phonon calculation for a specific volume if imaginary modes are detected. By setting it to 1, it means this correction will be attempted once. You may need to change this — read the logs. You may not be able to resolve all of the imaginary modes, but there should be relatively few and of low magnitude (again, check the logs). You should read  this paper  for an understanding of the impact that imaginary modes can have on the thermochemistry.
Note that QHA is not always perfect. If your system has rotational degrees of freedom, configurational entropy, or is metallic (i.e. electronic degrees of freedom), these will be neglected by QHA. Additionally, significant anharmonicity will not be properly treated without more advanced methods.

MD-Based Anharmonic Methods

To capture the anharmonicity missing from the QHA approximation, you can run an MD simulation and get the vibrational density of states from the MD trajectory itself. This is implemented in the DynaPhoPy code [ paper ,  code ]. The accuracy is expected to be improved compared to the QHA approximation.
Note that this method is also not perfect. In particular, it does not work well for floppy molecules as described several papers by Joachim Sauer, such as  this one . It should also not be used for (T, P) phase diagrams.

Other Notable Methods

  • SSCHA [ code ,  website ]
  • tdep [ code ]

Thermodynamic Integration (TI)

A separate approach to get the Gibbs free energy is through thermodynamic integration (TI). The  calphy code  has some automated routines for TI for use with LAMMPS, which many MLIP codes are interfaced with. Also, you may want to check out  this paper , which incorporates TI into the workflow code pyiron. Yet another code that supports TI is  dpti .
For calculating (T, P) phase diagram boundaries in solids, you want to find the free energy of each phase and report the lowest energy phase at each point. The simplest approach is to use the QHA. It is relatively straightforward, in that you would calculate G(T, P) for each phase under consideration at a given (T, P) pair. The phase with the lowest G(T, P) at a given set of conditions is the ground-state phase (an example for phases of ZrO2 can be found  here ). However, QHA is certainly not perfect as described in the  aforementioned review . TI resolves many of the problems with QHA and can be used to generate (T, P) phase boundaries, as described in the original  calphy paper .