18 February 2014

553. Briefly: Drawing molecules on Debian Wheezy: ISISDraw under wine (very briefly)

I'm cheating here: this isn't a program for linux. Instead we're running a windows program under wine, which is not ideal. On the other hand, it's so easy to do that it's worth exploring as an option to 'pure' linux offerings.

I'll be posting about native linux offerings later.

Either way, I've kept a copy of ISISDraw 2.5 around since the early 2000s. Luckily, other people are mirroring the installation file for it as well. See e.g. http://fc.smcdsb.on.ca/~rszerminski/addl_files_for_web/DRAW24.EXE for ISISDraw 2.4

 You'll need to install wine:i386 and for that you'll need to enable multiarch if you're on amd64:
sudo dpkg --add-architecture
sudo apt-get update
sudo apt-get install wine:i386
wget http://fc.smcdsb.on.ca/~rszerminski/addl_files_for_web/DRAW24.EXE
wine DRAW24.EXE

As always, if you have problems with e.g. weird fonts, use winetricks to install at least corefonts, or ideally allfonts:
wget http://winetricks.org/winetricks
sudo mv winetricks /usr/bin/winetricks
chmod +x /usr/bin/winetricks 
winetricks allfonts

Then run ISISDraw:
A blank slate
 ISISDraw has a lot of options, including templates and journal settings.
Chem Inspector
 There's a bug under wine: to get the menus for the buttons on the left side to show up, click on the button on the left panel, and then move the mouse left while still holding the button down.
To get the menus to work, click and while still holding the left button down move right
 ISISDraw is not going to blow you away, but it's free and it works under wine.

15 February 2014

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: Hcorr 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
Ee+Hcorr=Ee+Hvib+Hrot+Htrans

(where Ee is the electronic energy) which is an expression that doesn't include Δ(PV). In that case
Δ(H)=Σ(Hproducts)-Σ(Hreactants)+Δ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:
Hcorr=Etot+kBT
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=Ee+Hcorr=Ee+Etot+kBT]

At a first glance, kBT isn't equivalent to RT, but in fact kB=R/NA -- 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 Hcorr 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.

05 February 2014

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";
-}
+#}
--