EDIT: see here for a rough description of my take on a paper on reorganisational energies.
I'm reproducing papers to improve my understanding of what computational tools are available to 'wet' chemists. A good place to start for this is in the past -- typically, cutting edge methods are difficult to implement for someone who's not an expert in the field (Would you take a soldering iron to an NMR probe? We do it routinely, but most people wouldn't since it takes training.) and lacks the intuition to apply the right type of fudging. 'Old' papers (5-10 years) offer the advantage of having enough papers citing them that you get a crowd-sourced approach to evaluation of them, the methods described are often implemented directly in today's software packages and the computations are much, much cheaper with today's hardware -- often parts or even an entire paper can be reproduced in a day or two.
WARNING I may be doing something wrong WARNING
Here's my attempt at reproducing part of Baik and Friesner in J. Phys. Chem. A, 2002, 106, 7407-7412 - the exact approach is somewhat flavoured by the specific output that gaussian and nwchem return.
The basic principle is that if we know the ΔG of a redox reaction, then we can calculate the redox potential:
ΔG=-n*F*E,
where n is the number of electrons, F is Faraday's constant and E is the (unreferenced) redox potential. If we know ΔG in kcal/mol, then E=ΔG/(-n*F), where F is ca 23.06 kcal/(mol.V).
The main complication is that a normal dft calc will give us the electronic energy, ε0, while we ultimately need the thermodynamic paramaters, e.g. ΔH, ΔS, ΔG. Luckily, they are easy (although computationally costly) to get by doing a frequency calculation.
The short of it:
We here define G=ε0+Gcorr.
Gaussian will give you Gcorr directly, while e.g. nwchem makes you work for it -- but not too much -- so that
Gcorr=Hcorr-T*Srot+vib+trans
NOTE:In the paper (using Jaguar) the ZPE is explicitly included, but in nwchem and gaussian the Hcorr/Gcorr will already contain it. Approaches that I've seen online do not include ZPE and the computed data agrees better with experimental data if it's excluded. If I had to put money on anything, I'd say: make sure your software package includes ZPE in the Hcorr reported, and don't add it explicitly.
ΔGin gas phase is obtained by taking the difference in G of the product and reactant. All that remains is the solvation energy.
Throwing in solvation using a continuum solvent model is fairly easy, but can also lead to long computational times, so we will cheat -- the original paper is from 2002, when some of the calcs took them 40 days. So we will do a single-point/energy calculation to get the solvation.
(I think a more modern approach -- given more powerful hardware -- would be to do everything with a solvent model included for small molecules. For e.g. polyoxometalates, however, this can still be costly when using COSMO, but using PCM is almost a requirement when using g09 since anions converge very slowly. Bottomline: decide on the approach and be prepared to justify it.)
Anyway:
Gsolvation≅ (ε0(using solvation model)+Gcorr(gas phase))-(ε0(in gas phase)+Gcorr(gas phase))=ε0(using solvation model)-ε0(in gas phase)=εsolvation.
When I only write ε0 I mean in gas phase. Here's for a reaction where A goes to B through a redox process:
ΔG≅(ε0+εsolvation+Gcorr)B-(ε0+εsolvation+Gcorr)A=ΔGin gas phase+ΔGsolvation
How:
You need to do the following steps:
- Draw benzophenone
- Optimise the structure in the gas phase i.e. without a solvent model
- Do a frequency calc of the optimised structure in the gas phase
- Do an energy (single point) calculation of the gas-phase optimised structure using a solvent model (e.g. PCM or COSMO). That is, you use the structure you optimised in the gas phase, then apply a solvent model and do a single-point energy calculation.
- Draw the benzophenone anion
- Optimise the structure in the gas phase i.e. without a solvent model
- Do a frequency calc of the optimised structure in the gas phase
- Do an energy (single point) calculation of the gas-phase optimised structure using a solvent model
2 will give you ε0. 3 will give you ε0 again, and Gcorr. 4 will give you ε0(using solvation model).
Steps 5-8 will do the same, but for the reduced species.
Gaussian inputs:
#p opt rb3lyp/6-31g**and
benzophenone, neutral, gas phase
0 1
...xyz coordinates of starting guess...
#p freq rb3lyp/6-31g**and
benzophenone, neutral, frequency
0 1
...xyz coordinates of optimised structure...
#p rb3lyp/6-31g** scrf=(pcm,solvent=acetonitrile)and
benzophenone, neutral, solvation
0 1
...xyz coordinates of optimised structure...
#p opt rb3lyp/6-31g**
benzophenone, anion, gas phase
-1 2
...xyz coordinates of starting guess...
etc. Note the multiplicity -- one unpaired electron gives a multiplicity of 2 and a charge of -1. You can of course use the .chk file as a restart point, or do both opt and freq in a single calculation.
Nwchem input:
title "Benzophenone, neutral, gas phase"and
start benzophenone_gas_opt
charge 0
geometry
...xyz coordinates of starting guess...
end
dft
xc b3lyp
mult 1
grid fine
end
task dft optimize
title "benzophenone, neutral, gas phase, frequency"and
start benzophenone_gas_freq
charge 0
geometry
...xyz coordinates of optimised structure...
end
dft
xc b3lyp
mult 1
grid fine
end
task dft freq
title "benzophenone, neutral, COSMO"and
start benzophenone_solv_energy
charge 0
geometry
...xyz coordinates of optimised structure...
end
cosmo
dielectric 37.5
end
dft
xc b3lyp
mult 1
grid fine
end
task dft energy
title "Benzophenone, anion, gas phase"etc.
start benzophenone_anion_gas_opt
charge -1
geometry
...xyz coordinates of starting guess...
end
dft
xc b3lyp
mult 2
grid fine
end
task dft optimize
You can of course throw two task directives in one and the same calculcation etc. You can also re-use movecs files. But I'm trying to make this as simple as possible to follow.
Gaussian output:
Neutral, gas phase opt:
SCF Done: E(UB3LYP) = -576.648098680 A.U. after 16 cycles
Convg = 0.3611D-08 -V/T = 2.0097
Neutral, gas phase freq:
Zero-point correction= 0.191763 (Hartree/Particle)
Thermal correction to Energy= 0.202482
Thermal correction to Enthalpy= 0.203427
Thermal correction to Gibbs Free Energy= 0.154218
Sum of electronic and zero-point Energies= -576.456336
Sum of electronic and thermal Energies= -576.445616
Sum of electronic and thermal Enthalpies= -576.444672
Sum of electronic and thermal Free Energies= -576.493881
neutral, solvation=acetonitrile, energy:
SCF Done: E(UB3LYP) = -576.654962012 A.U. after 13 cycles
Convg = 0.2187D-08 -V/T = 2.0097
Anion gas phase opt:
SCF Done: E(UB3LYP) = -576.656219870 A.U. after 20 cyclesAnion gas frequency:
Convg = 0.6297D-08 -V/T = 2.0095
Zero-point correction= 0.187970 (Hartree/Particle)
Thermal correction to Energy= 0.198790
Thermal correction to Enthalpy= 0.199734
Thermal correction to Gibbs Free Energy= 0.150074
Sum of electronic and zero-point Energies= -576.468250
Sum of electronic and thermal Energies= -576.457430
Sum of electronic and thermal Enthalpies= -576.456486
Sum of electronic and thermal Free Energies= -576.506146
Anion, solvation, energy:
SCF Done: E(UB3LYP) = -576.732948458 A.U. after 16 cycles
Convg = 0.5837D-08 -V/T = 2.0096
Putting it all together:
ε0,anion,gasphase=-576.656219870
Ganion,gasphase=-576.656219870+0.150074=-576.506145870
εsolvation,anion=-576.732948458-(-576.656219870)=-0.076728588
ε0,anion,gasphase=-576.656219870
Ganion,gasphase=-576.656219870+0.150074=-576.506145870
εsolvation,anion=-576.732948458-(-576.656219870)=-0.076728588
ε0,neutral,gasphase=-576.648098680
Gneutral,gasphase= -576.648098680+0.154218=-576.493880680
εsolvation,neutral=-576.654962012-(-576.648098680)=-0.006863332
Gneutral,gasphase= -576.648098680+0.154218=-576.493880680
εsolvation,neutral=-576.654962012-(-576.648098680)=-0.006863332
ΔG=(-576.506145870-0.076728588)-(-576.493880680-0.006863332)=-0.082130446 Hartree=-51.537594039014 kcal/mol=ΔG/(-1*F)= -51.537594039014/(-1*23.06)=2.23493469379939288811 V. Referencing it against the SCE at 4.1888 V, we get -1.95386530620060711189 V vs SCE
Looking at the paper, they got -2.426+4.1888=1.7628 V absolute, corresponding to -40.6502 kcal/mol=-0.06478026 hartree. Remember that this was done using G09 and not Jaguar, which was used in the paper. The experimental value is -1.88 V. The maximum obtainable accuracy of DFT (according to the paper) is about 2 kcal/mol.
Nwchem output:
Neutral, gas
Total DFT energy = -576.648088887542neutral, gas, freq
One electron energy = -2313.472465671473
Coulomb energy = 1046.450004818160
Exchange-Corr. energy = -82.832342635322
Nuclear repulsion energy = 773.206714601093
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)
Total Entropy = 103.010 cal/mol-K
- Translational = 41.486 cal/mol-K (mol. weight = 182.0732)
- Rotational = 31.544 cal/mol-K (symmetry # = 1)
- Vibrational = 29.980 cal/mol-K
neutral, solvated, energy
COSMO solvation results
-----------------------
gas phase energy = -576.6480876954
sol phase energy = -576.6603289676
(electrostatic) solvation energy = 0.0122412722 ( 7.68 kcal/mol)
Total DFT energy = -576.656222314875
One electron energy = -2320.751987024161
Coulomb energy = 1058.312412372760
Exchange-Corr. energy = -83.108633373896
Nuclear repulsion energy = 768.891985710422
Anion gas, freq
Zero-Point correction to Energy = 118.177 kcal/mol ( 0.188327 au)
Thermal correction to Energy = 124.894 kcal/mol ( 0.199031 au)
Thermal correction to Enthalpy = 125.486 kcal/mol ( 0.199975 au)
Total Entropy = 102.152 cal/mol-K
- Translational = 41.486 cal/mol-K (mol. weight = 182.0732)
- Rotational = 31.577 cal/mol-K (symmetry # = 1)
- Vibrational = 29.090 cal/mol-K
Anion solvated, energy
gas phase energy = -576.6562158410
sol phase energy = -576.7436282813
(electrostatic) solvation energy = 0.0874124403 ( 54.85 kcal/mol)
Which works out to
Gcorr=(127.706-298.15*103.010/1000)/627.509=0.15456920697551748261
ε0,neutral,gasphase=-576.648088887542
Gneutral,gasphase=-576.648088887542+0.15456920697551748261+0.191895=-576.49351968056648251739+0.191895
εsolvation,neutral=0.0122412722
εsolvation,neutral=0.0122412722
Gcorr=(125.486-298.15* 102.152/1000)=0.15141494338035305421
ε0,anion,gasphase= -576.656222314875
Ganion,gasphase= -576.656222314875+0.15141494338035305421+0.188327=-576.50480737149464694579+0.188327
εsolvation,anion=0.0874124403
ΔG=(-576.50480737149464694579-0.0874124403)-(-576.49351968056648251739-0.0122412722)=-0.08645885902816442840 Hartree=-54.25371216990443230085 kcal/mol => 2.35271952167842250687 V (absolute) <=> -1.836 V vs SCE. This is very close to the experimental value of -1.88 V but not at all close to the reported value in the paper of -2.426 V vs SCE.
Discussion:
Whether I'm doing something wrong or whether the differences in solvation models and implementation in differernt software packages is to blame, I don't know. 15 kcal/mol is a lot.
Ganion,gasphase= -576.656222314875+0.15141494338035305421+0.188327=-576.50480737149464694579+0.188327
εsolvation,anion=0.0874124403
ΔG=(-576.50480737149464694579-0.0874124403)-(-576.49351968056648251739-0.0122412722)=-0.08645885902816442840 Hartree=-54.25371216990443230085 kcal/mol => 2.35271952167842250687 V (absolute) <=> -1.836 V vs SCE. This is very close to the experimental value of -1.88 V but not at all close to the reported value in the paper of -2.426 V vs SCE.
Discussion:
Whether I'm doing something wrong or whether the differences in solvation models and implementation in differernt software packages is to blame, I don't know. 15 kcal/mol is a lot.
I'll post data for larger basis sets in a week or so to see whether the nwchem/cosmo value was a fluke or 'real'.
[g09 using aug-cc-pwdCVDZ:
single point E, gas phase using lanl2dz structure:
neutral : -576.718523917
anion: -576.746041821
ΔG=(-576.746041821+0.150074-0.076728588)-(-576.718523917+0.154218-0.006863332)= -0.101527160 Hartree =-1.426 V relative to SCE, vs -1.88 experimental and -1.953 with LANL2DZ. That didn't help...]
I should've looked this up ages ago haha...very helpful!
ReplyDeleteVery nice and clear instruction. Halim UWO.CA
ReplyDeletebase on the thermodynamic cycle you need to calculate solvation free energies for both reduced and oxidized forms.
ReplyDeleteThanks for the feedback.
DeleteThat's right. And as far as I can see that's what I did in the post -- let me know where it is missing and I'll fix the post.
Everything is fine... my bad... i missed out something...
DeleteAbout gaussian, do i have to add ZPE to the thermal correction to free energy?
This is a good source of information for thermochemistry in gaussian: http://www.gaussian.com/g_whitepap/thermo/thermo.pdf
DeleteIn that document, the ZPE is only used for the atomization energy calculation, but not the enthalpy/ free energy calcs.
Can you please help me. What changes in the calculation do I have to make if instead of calculating the reduction potential, I will calculate the oxidation potential. I understand that the Delta G will be that corresponds to the Free energy of oxidation. I am confused with regards to the substitution of this value to the E=-DeltaG/nF... I tried some calculations and the results were questionable as to its value.
ReplyDeleteAt its simplest the oxidation reaction is simply the reduction reaction reversed.
DeleteOx -> Red
Red -> Ox
Same goes for the calculation of Ox vs Red potentials -- just change the sign of DG.
As for questionable values: assuming that you are aware of DG giving the ABSOLUTE, not relative, potential, it may simply be due to the inadequacy of the solvation model, or poorly optimised structures, or the real reaction proceeding in a more complex way than what you have modelled.
Can you give me a simple example of what you are trying to do? E.g. Benzene -> Benzene^- vs Benzene -> Benzene^+?
The system that i've been working has an experimental value of 0.48V oxidation potential (vs Ag/AgCl). Hence, the reaction would be
Delete[MnL] -> [MnL]^+ +e-.
With regards to the changing of the sign fo DG, does is it mean that the equation will now be:
E=DG/nF instead of E=-DG/nF? and DG is the free energy of ionization which has always a positive value?
Please try to see if what I did is correct,
Delete-the reported value is 0.49 V oxidation potential.
-the reaction i presume shud be : [ML] ---> [ML]^+ + e^-
-and so I shud calculate DG as:
DG = G([ML]^+) - G([MN])
-then use the DG above to this equation
E=DG/nF <--- is this correct? the sign of DG here is + ?
No, if you write the reaction like this it should be -DG -- i.e. what I meant was that the free energy difference for the oxidation potential is the same as DG for the reduction potential, but with a sign change.
DeleteRemember that DG gives you the absolute potential, whereas the expt value is a relative value. You'll need to find the absolute potential for Ag/AgCl and correct for that.
Thanks a lot, very useful information.
ReplyDeleteIn Gaussian:
ReplyDelete1. We look for Gibbs's free energy but when you calculate solvation energy you don't consider thermochemistry calculations. Why? Is it the same?
2. Why do not consider molecule relaxation in solvent? Is not necessary?
3. How to include thermal contributions in the IP or EA calculations? Are they included in the Gcorr?
I really want to say that your Blog is awesome!!! It's so useful!
Best Regards
In Gaussian:
ReplyDelete1. We look for Gibbs's free energy but when you calculate solvation energy you don't consider thermochemistry calculations. Why? Is it the same?
2. Why do not consider molecule relaxation in solvent? Is not necessary?
3. How to include thermal contributions in the IP or EA calculations? Are they included in the Gcorr?
I really want to say that your Blog is awesome!!! It's so useful!
Best Regards
The answer to 2 is probably the easiest and most important: when using COSMO (and PCM /I think/) you should never use it during geometry optimisation. The method was only designed to provide solvation energies for _gas phase_ structures -- you shouldn't use COSMO/PCM during geometry optimisation as it wasn't designed for this.
DeleteOn the other hand, a lot of people DO use COSMO/PCM in that way, but it just isn't 'right' (note that that says nothing about the accuracy about either the energies obtained or the faithfulness of the structure).
Also, the more I learn (and I've got a long, long way to go) the less I trust implicit solvation models and the more hesitant I become in using data obtained with them to draw any conclusions. If the gas phase energy trends and the COSMO/PCM trends diverge at all you'll need to spend some time hammering out exactly why that is. And given the complexity of typical computational targets today you're likely better off not playing with implicit solvation models if solvation is important.
Because here's the final issue: implicit solvation models like PCM and COSMO are often all but useless for anything other than neutral molecules which don't take part in any special interactions with the solvent (e.g. pi stacking, hydrogen bonding).
If solvent interaction is key, then you'll need to include explicit solvent molecules.
Either way, the answer to question 2 is: at least for COSMO you should only apply it to the gas phase geometries.
Answer to 3: yes, Gcorr=H-T*S, so that's the most immediate way of correcting for temperature.
DeleteNote that that assumes that H and S don't change with temperature -- I think they are calculated assuming T=298.15 K in nwchem and presumably something similar in G09, but I'm just speculating here.
I'm not entirely sure what you mean by question 1, but I think I answered it as part of the answer to question 2 -- the implicit solvation model just gives an energy correction term for the gas phase structure.
DeleteIn the first question I try to say that we use the SCF energy to calculate solvation free energy.
DeleteStrictly, we must to use the thermochemical function using:
DeltaG (solvation)= Gcorr(solv) - Gcorr(gas)
I don't know if this relation is correct. What do you think? Is correct relate SCF energy with Free energy in a direct way (without considering thermal correction to free energy)?
I'm really grateful for your answers
Sorry for the delayed reply.
DeleteWhereas the solvation correction tends to be referred to as the free energy of solvation, I think a more appropriate (practical) way of looking at it is as a correction to the electronic energy i.e. you can use the solvation correction even if you're just doing an electronic calculation.
So whereas I think that what you say makes sense, in practice the solvation is used as a correction to the electronic energy.
Amazing work. Is this code can calculate reduction potential of solvent?
ReplyDelete