Solid-State Thermochemistry

Solid-State Thermochemistry

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.

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. This should be considered required reading. The methods below are generally too expensive to do with ab initio methods and will require an MLIP.
Some rules of thumb:
  • For quick analyses of free energies of solids, QHA is usually the first method to try and can even be fast enough to use with DFT in many cases — but it struggles to capture anharmonicity. This method is best suited for solids with stiff bonds at moderate temperatures.
  • For a more rigorous treatment of solid-state free energies than QHA, SSCHA is typically a better way to go. It is still based on vibrational degrees of freedom via phonons but can capture anharmonicity.
  • If you want the highest quality free energies for solids or are dealing with liquids where phonons are ill-defined, use TI if you can afford it. Convergence may be tricky though.

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

SSCHA

One way to address the anharmonicity issues of the (Q)HA approximation is called SSCHA [ code ,  website ]. This method is nicely demonstrated with MLIPs in  this paper . It is well-suited for (T, P) phase diagrams of solids but should not be used for liquids.

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 .