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.

551. Very Briefly: Getting ECCE to work with Gaussian 09 (G09) part 3: Energy gradients (EGRADVEC)

This is a really simple fix to make ECCE return the energy gradient vectors for geometry optimisation jobs in G09:

Edit apps/scripts/parsers/gaussian-03.egradvec and comment out lines 38 (if () {... ) and 70 (}):

diff --git a/ecce-v7.0/scripts/parsers/gaussian-03.egradvec b/ecce-v7.0/scripts/parsers/gaussian-03.egradvec
index a2f5633..4a4b691 100755
--- a/ecce-v7.0/scripts/parsers/gaussian-03.egradvec
+++ b/ecce-v7.0/scripts/parsers/gaussian-03.egradvec
@@ -35,7 +35,7 @@ $| = 1;
 #Don't parse gradients for runtypes with geometry optimization
 #since we can't match gradients with correct geometry.
 
-if (!($runtype =~ /Geo/i)) {
+#if (!($runtype =~ /Geo/i)) {
   $natom = 0;
   while () {
     if (/-----/) { last; }
@@ -67,4 +67,4 @@ if (!($runtype =~ /Geo/i)) {
   }
   print "units:\nHartree/Bohr\n";
   print "END\n";
-}
+#}
-- 

550. Briefly: Getting ECCE to work with Gaussian 09 (G09) part 2: Geometry trace when using NoSymm

I've made a couple of posts about how to modify ECCE to work better in general, but in particular with Gaussian 09  (G09).

See here for some of the previous posts:
Adding new basis sets, part 1: http://verahill.blogspot.com.au/2013/06/455-adding-nwchem-basis-sets-to-ecce.html
Adding new basis sets, part 2:  http://verahill.blogspot.com.au/2013/06/456-adding-nwchem-basis-sets-to-ecce.html
Adding new exchange/correlation functionals:
http://verahill.blogspot.com.au/2013/12/533-adding-new-exchangecorrelation.html
Adding new options to ECCE
http://verahill.blogspot.com.au/2013/12/534-adding-new-options-to-ecce-adding.html
Getting ECCE to read frequency calc output from G09 calcs:
http://verahill.blogspot.com.au/2013/12/536-briefly-getting-ecce-to-work-with.html


Either way, I've noticed that when I'm using the NoSymm option during optimisation ECCE fails to import the geometries from the geometry steps.

No geometry trace option

What is even worse is that if you use an optimisation job like that as the basis for another job by selecting Duplicate Setup with Last Geometry ECCE will in fact use the starting geometry -- not the optimised one.

Luckily the fix is simple -- edit apps/parsers/gaussian-03.desc and add the following lines:

[GEOMTRACE]
Script=gaussian-03.geomtrace
Begin=Z-Matrix orientation\:
Skip=5
Frequency=all
End=--------------------------------------
[END]

You can now 'Reconnect job monitoring' also to your old jobs if you still have access to the run directory. Either way, the view will now be a sweet one:

Et voila: A geometry trace