How to evaluate the heat of formation with computational chemistry?

Question

I'm sorry if this is a stupid question, but since I'm a beginner with computational chemistry, please, be patient.

Let's consider the reaction of formation of water:

$$\ce{1/2 O2(g) + H2(g) -> H2O(g)}$$

This reaction shoud give a $\Delta_\mathrm fH^\circ = -241.83~\mathrm{kJ/mol}$. The reaction is intentionally simple to let me run fast calculations and tests with single point energies.

First, I've performed a geometry optimization on $\ce{H2O}$, $\ce{O2}$ and $\ce{H2}$ (HF/6-31G*). Then, I've used M06/Def2-QZVPP for the single points:

• $E_\mathrm{SP}(\ce{H2O}) = -76.431519932791~E_\mathrm h$
• $E_\mathrm{SP}(\ce{O2}) = -150.262365177069~E_\mathrm h$
• $E_\mathrm{SP}(\ce{H2}) = -1.171819656537~E_\mathrm h$

In this case, doing $$E_\mathrm{SP}(\ce{H2O}) - \left(\frac12 E_\mathrm{SP}(\ce{O2}) + E_\mathrm{SP}(\ce{H2})\right) = \Delta E$$ I'm getting $\Delta E = -0.1285176877195084~E_\mathrm h$ that is equal to $-337.42312874~\mathrm{kJ/mol}$, almost $96.6~\mathrm{kJ/mol}$ higher than the reference of $-241.83~\mathrm{kJ/mol}$.
Is this the correct approach for a case like this? Or is this discrepancy normal?

Answer 1

Your approach is quite correct, but as Jan already pointed out, it is incomplete. What you calculated is the difference in the electronic energy of the reaction $$\ce{H2 (g) + 1/2 O2 (g,\,{}^1\Delta_{g}) -> H2O (g)}.$$

Let's turn this into some kind of a tutorial and a little exercise, that you can try to reproduce at home. For a more detailed description on how it works in Gaussian, you can read the paper by Joseph W. Ochterski. I will use Gaussian09 rev. D for the calculations involved.

The enthalpy of formation of water from the standard state of the elements involved, is basically the reaction enthalpy of the combustion of hydrogen and oxygen. You are in luck, because most of the quantum chemical programs shall provide you with the necessary data for reactions straight away. What you do need to do is calculate the thermal corrections to the reaction. Essentially this is done by calculating the Hessian, i.e. a frequency calculation. The second derivatives only have a interpretable meaning if they are performed at a stationary point. Hence they must be performed with the same method, that you used to optimise the geometry.
Jan already pointed out, that HF/6-31G(d) is probably not a well enough fit for geometries and frequencies. I'd always consider some level of theory with at least some coverage of correlation energy. I often prefer BP86 or M06L and the def2-SVP basis set as a first anchor point. The smaller the system, the more accurate the methods you can use (or afford).
Jan also pointed out, that especially for oxygen, the spin state is important. (I did a quick check and am almost 100% certain you calculated it in the singlet state; see above equation.)

For this exercise I have matched your suggestion with using M06/def2-QVZPP as it is still feasible. Keep in mind, that most of the times density functional approximations have to be calibrated. It's all about what you can or cannot afford. You should perform calculations to the best of your abilities.

Okay, let's dive into some mathematics. First, we consider the reaction equation $$\ce{H2 (g) + 1/2 O2 (g,\,{}^3\Sigma_{g}^{-}) -> H2O (g)}.$$

And we would like to know the enthalpy of formation for water at standard state (which we will use as $298.15~\mathrm{K}$ and $1~\mathrm{atm}$) $$\begin{align} \Delta_{f}H^\circ(\ce{H2O}) &= H^\circ(\ce{H2O}) - (H^\circ(\ce{H2}) +\frac12 H^\circ(\ce{O2}))\\ \Delta_{f}H^\circ(\ce{H2O}) &= [E_\mathrm{el}(\ce{H2O}) + H_\mathrm{corr}(\ce{H2O})]\\ &\qquad - \left([E_\mathrm{el}(\ce{H2}) + H_\mathrm{corr}(\ce{H2})] +\frac12 [E_\mathrm{el}(\ce{O2}) + H_\mathrm{corr}(\ce{O2})]\right)\\ \end{align}$$

The data in the expanded form we can directly take from the output of the frequency calculation:

Oxygen  mo6qzvpp.freq.log/.com
 SCF Done:  E(UM06) =  -150.322374247     A.U. after    1 cycles
 Zero-point correction=                           0.003897 (Hartree/Particle)
 Thermal correction to Enthalpy=                  0.007203
 Thermal correction to Gibbs Free Energy=        -0.016049

Hydrogen mo6qzvpp.freq.log/.com
 SCF Done:  E(RM06) =  -1.17221845130     A.U. after    1 cycles
 Zero-point correction=                           0.009902 (Hartree/Particle)
 Thermal correction to Enthalpy=                  0.013206
 Thermal correction to Gibbs Free Energy=        -0.001590

Water mo6qzvpp.freq.log/.com
 SCF Done:  E(RM06) =  -76.4321360023     A.U. after    1 cycles
 Zero-point correction=                           0.021635 (Hartree/Particle)
 Thermal correction to Enthalpy=                  0.025415
 Thermal correction to Gibbs Free Energy=         0.004010

Therefore we will obtain $$\Delta_{f}H^\circ(\ce{H2O}) = -236.6~\mathrm{kJ/mol}$$ which I think is reasonably close to the experimental value.

---

Appendix The following input files for the computation of above values.

%chk=mo6qzvpp.chk
#p M06/def2QZVPP
opt
int(ultrafinegrid)

oxygen

0 3
O
O 1 1.2
 (blank)

---

%chk=mo6qzvpp.freq.chk
%oldchk=mo6qzvpp.chk
#p M06/def2QZVPP
freq
geom=check guess=read
int(ultrafinegrid)

oxygen freq

0 3
 (blank)

---

%chk=mo6qzvpp.chk
#p M06/def2QZVPP
opt
int(ultrafinegrid)

hydrogen

0 1
H
H 1 1.2
 (blank)

---

%chk=mo6qzvpp.freq.chk
%oldchk=mo6qzvpp.chk
#p M06/def2QZVPP
freq
geom=check guess=read
int(ultrafinegrid)

hydrogen freq

0 1
 (blank)

---

%chk=mo6qzvpp.chk
#p M06/def2QZVPP
opt
int(ultrafinegrid)

water

0 1
O
H 1 1.2
H 1 1.2 2 109.0
 (blank)

---

%chk=mo6qzvpp.freq.chk
%oldchk=mo6qzvpp.chk
#p M06/def2QZVPP
freq
geom=check guess=read
int(ultrafinegrid)

water freq

0 1
 (blank)

---

Have a lot of fun.

Answer 2

1. you need to perform a Hessian calculation and add the enthalpy correction to each energy
2. what spin state are you using for $\ce{O2}$?
3. HF/6-31G(d) is probably not good enough for geometries and frequencies.
4. Depending on what accuracy you're after, M06/Def2-QZVPP may not be good enough

Answer 3

The problem is that such a simple calculation has several associated errors (geometries, energies, vibrational frequencies and, therefore, thermodynamic functions) which render the result very inaccurate.

In the case of water, using the same M06 method you used, one obtains -212.9 kJ/mol with the cc-pVDZ basis set and -246.0 kJ/mol with the much larger aug-cc-pV6Z basis set. Using a much more accurate method, CCSD(T), the enthalpy of formation with the same previous two basis sets are -200.19 and -241.89 kJ/mol. The latter result is nearly identical to the experimental one.

Incidentally, notice that quoting enthalpies of formation with more than one decimal digit is only justified if you are using highly accurate methods. You will very rarely have accuracies of better than 1 kJ/mol and certainly not using DFT methods.