Calculating enthalpy of a chemical reaction using DFT

Question

Let's say we have a decomposition reaction

$\tag{1}A \rightarrow B + C$

how would I go about calculating the Enthalpy of this chemical reaction using VASP. I am gonna describe what I am thinking below please let me know what you guys think. We are actually trying to see the thermodynamic stability of A. So we have the structure(cif file) for A in which we want to see the stability, B and C are compounds which have been experimentally synthesized so we will use the most common structures.

We start by doing geometrical optimization on all the compounds and then running an SCF calculation using the converged system(POSCAR). Now the ground state energy found from this will be used,

$\Delta H = H_{products} - H_{reactants}\tag{2}$

and so

$\Delta H = H_{B} + H_{C} - H_{A}\tag{3}$

where $H_i$ is the ground state energy obtained from the SCF calculation for i.

Now as VASP does calculation with $T = 0$ and,

$\tag{4}\Delta G = \Delta H - T \Delta S$ we get:

$\tag{5}\Delta H = \Delta G$ so our enthalpy can be used to predict the stability.

Let me know if I am going wrong somewhere or not. Thanks!

Answer 1

The main question here is whether the reaction is favourable or not, and for that we may need to go beyond the usual energies computed in a ground state DFT calculation, and consider the free energy. DFT's "total energy" is the same as the "internal energy", $E$, in thermodynamics, and you can compute this for each reactant and product, as you note. However, in order to know whether a reaction is favourable under different experimental conditions, we might also need the $PV$ and $TS$ terms in order to calculate the appropriate free energy, where $P$ is the pressure, $V$ is the volume, $T$ is the temperature and $S$ is the entropy.

Pressure (stress)

The $PV$ term is straightforward, since the stress tensor is computable using perturbation theory, directly from the DFT solution - the relevant perturbation being an infinitesimal strain of the cell. Calculating the stress tensor does usually require an additional term (called the Pulay stress) due to how the basis set changes when the cell is strained, but this is well-known for plane-wave basis sets and ultrasoft/PAW formalisms such as those used by ABINIT, CASTEP, Quantum Espresso and VASP. The PV term is sufficient to calculate the enthalpy, $H = E + PV$.

Note that if we have a particular environmental pressure, $P$, then we should optimise the geometry of our system under this applied pressure (positive or negative), for which we will have to do a full variable-cell optimisation; i.e. the optimisation should be performed by varying the cell vectors, in addition to the atomic coordinates. In the past, not all DFT software was capable of variable-cell optimisation, so people sometimes used an energy-volume (E-V) plot and fit a Birch-Murnaghan equation of state; this is not recommended, as it is extremely prone to error, and is not as accurate as a proper optimisation approach.

Entropy

The Helmholtz free energy, $F$, is defined as $F = E - TS$, so requires a calculation of $TS$; this is usually very difficult! $S$ contains contributions from the vibrations of the system around its equilibrium ("vibrational entropy") as well as contributions from all of the different ways the atoms could be arranged ("configurational entropy"). Vibrational entropy can be calculated in the quasi-harmonic approximation by computing the $\sim~3N$ vibrational normal modes of your $N$-atom system and thermally populating them; these modes are the "phonons" in a periodic material system, but in molecules you may need to consider other modes such as rotational modes. The configurational entropy is extremely difficult to compute, to the point where it is impractical in most cases. This is one of the reasons why $T = 0$K is so much easier to work with, because the entropic contribution is zero.

If we have both $PV$ and $TS$ terms then we can calculate the Gibbs free energy, $G = E + PV - TS$.

Conclusions

At low temperatures, using the DFT energy ($E$) or enthalpy ($H = E + PV$), as appropriate, can often give a good indication of how favourable a reaction will be, and what energy will be released (or required). At higher temperatures, the entropic contribution (via $TS$) becomes significant, and this is very difficult to model accurately. It is possible to calculate the vibrational contributions to the entropy in the quasiharmonic approximation, via a phonon calculation, which is relatively accurate for modest temperatures, but it is not usually practical to calculate the configurational entropy.