552. Very briefly: Enthalpy of correction and different number of reactants and products in nwchem and gaussian

Those who are well-versed in the computational arts won't get much out of this post. On the other hand, happy amateurs like myself might find it a bit more useful as it clarifies something that's been bothering  me.

Short version: H_corr in Nwchem and Gaussian include PV.

As usual, I am not an expert in either computational or theoretical chemistry. I try to use it as a tool, and I try to use it as well as I can. But I am not an authority. Also, if you consult older posts on the blog you'll find plenty of examples of me misunderstanding basic computational concepts (with 550+ posts it's difficult to go back and erase all the embarrassing little gems).

---

The background
I had a bit of a fright the other day. I'm currently working on computationally characterising a system which undergoes a reaction that can lead to a large number of isomers. Only one of them is experimentally observed, and so it was of interest to see whether this is the thermodynamically favoured product as predict by reasonably cheap computational methods.

Because DFT calculations like these are based on gas phase reactions (even if you use a solvation model) the free energy that you get is based on the standard conditions for gas phase corrections i.e. 1 bar of partial pressure of each reactant.

If you want to calculate the free energy of reaction in solution you need to use concentrations instead. As you will see, you'll only have to worry about this for reactions that involve an unequal number of reactants and products. Normally your best results will be obtained by using isodesmic reaction schemes anyway, which is a great way of avoiding this. However, if you do have unequal numbers of reactants and products you /must/ correct for it when making solution phase predictions.

A gas at 1 bar of pressure is ca

V=n*R*T/P= 1*8.314*298.15/101350=0.024458 m^3= 24.46 litres

So for an example reaction like this:

A+B -> C

Using

G_SATP=G(gas phase) + R*T*ln Q
=G(gas phase)+ R*T*ln([C]/([A]*[B])
=G(gas phase)+ R*T*ln((1/24.46)/((1/24.46)*(1/24.46)))
=G(gas phase)+1.89 kcal/mol 

In other words: you can't ignore the standard state change when doing solution phase calculations. This is obviously of extra importance in pH calculations, which are notoriously tricky.

---

Enthalpy
So knowing the need for standard state corrections I was a bit paranoid about how to treat the reaction enthalpy and came across this document: http://chemistry.illinoisstate.edu/standard/che38037/handouts/380.37assign3.pdf

On page, equation 2 states that for Gaussian

E_e+H_corr=E_e+H_vib+H_rot+H_trans

(where E_e is the electronic energy) which is an expression that doesn't include Δ(PV). In that case

Δ(H)=Σ(H_products)-Σ(H_reactants)+Δn*R*T

where Δn is the difference in number of moles of products and reactants.

Consulting the excellent Gaussian thermochemistry whitepaper (http://www.gaussian.com/g_whitepap/thermo.htm) offers the following:

H_corr=E_tot+k_BT

and

The Gibbs free energy includes Δ(PV) , so when it's applied to calculating ΔG for a reaction, ΔNRT=ΔPV is already included. This means that ΔG will be computed correctly when the number of moles of gas changes during the course of a reaction.

[Note that H=E_e+H_corr=E_e+E_tot+k_BT]

At a first glance, k_BT isn't equivalent to RT, but in fact k_B=R/N_A -- in the words of Wikipedia: "[The gas constant R] is equivalent to the Boltzmann constant, but expressed in units of energy".

In other words, Δ(PV) is already accounted for in H_corr in gaussian.

Somewhat clearer: the Freq() page, http://www.gaussian.com/g_tech/g_ur/k_freq.htm, on the gaussian website now states

Sum of electronic and zero-point Energies=-527.492585  E0=Eelec+ZPE
Sum of electronic and thermal Energies= -527.489751     E= E0+ Evib+ Erot+Etrans
Sum of electronic and thermal Enthalpies=-527.488807  H=E+RT

We test that claim:
E=-527.489751 Hartree
RT=(1.987*298.15/1000)/627.503=9.4410e-04 Hartree
E+RT=-527.489751+9.4410e-04=-527.488807

It holds up. 

A similar example from an nwchem calculation:

Zero-Point correction to Energy = 120.416 kcal/mol ( 0.191895 au)
Thermal correction to Energy = 127.114 kcal/mol ( 0.202569 au)
Thermal correction to Enthalpy = 127.706 kcal/mol ( 0.203513 au)

(0.203513-0.202569)=9.4400e-04 Hartree(1.987*298.15/1000)/627.503

Nwchem also includes Δ(PV) in the thermal correction to enthalpy.