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.
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.
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. from matcalc import QHACalc
import matplotlib.pyplot as plt
logging.basicConfig(level=logging.INFO)
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)
thermal_expansion = results["thermal_expansion_coefficients"]
c_P = results["heat_capacity_P"]
results["qha"].plot_qha()
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 ( allow_shape_change=True ).
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.
To use QHA without matcalc for maximum flexibility, you can build upon the example below:
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. 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 . Some additional resources: