164. A rough approach to calculating redox potentials using dft

I'm not a computational or theoretical chemist, but I'd like to expand my skills in the area. Ergo, if you do post, don't ask me questions about theory -- I'm only at the learning stage myself. In addition, this is my interpretation of the paper I'm following.

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, ε₀, 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=ε₀+G_corr.
Gaussian will give you G_corr directly, while e.g. nwchem makes you work for it -- but not too much -- so that
G_corr=H_corr-T*S_{rot+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.

ΔG_{in 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:
G_solvation≅ (ε₀(using solvation model)+Gcorr(gas phase))-(ε₀(in gas phase)+Gcorr(gas phase))=ε₀(using solvation model)-ε₀(in gas phase)=ε_solvation.

When I only write ε₀ I mean in gas phase. Here's for a reaction where A goes to B through a redox process:
ΔG≅(ε₀+ε_solvation+G_corr)_B-(ε₀+ε_solvation+G_corr)_A=ΔG_{in gas phase}+ΔG_solvation

How:
You need to do the following steps:

1. Draw benzophenone
2. Optimise the structure in the gas phase i.e. without a solvent model
3. Do a frequency calc of the optimised structure in the gas phase
4. 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.


5. Draw the benzophenone anion
6. Optimise the structure in the gas phase i.e. without a solvent model
7. Do a frequency calc of the optimised structure in the gas phase
8. Do an energy (single point) calculation of the gas-phase optimised structure using a solvent model

2 will give you ε₀. 3 will give you ε₀ again, and G_corr. 4 will give you ε₀(using solvation model).
Steps 5-8 will do the same, but for the reduced species.

---

Gaussian inputs:

#p opt rb3lyp/6-31g**

benzophenone, neutral, gas phase

0 1
...xyz coordinates of starting guess...

and

#p freq rb3lyp/6-31g**

benzophenone, neutral, frequency

0 1
...xyz coordinates of optimised structure...

and

#p rb3lyp/6-31g** scrf=(pcm,solvent=acetonitrile)

benzophenone, neutral, solvation

0 1
...xyz coordinates of optimised structure...

and

#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"
start benzophenone_gas_opt
charge 0
geometry
       ...xyz coordinates of starting guess...
end
dft
       xc b3lyp
       mult 1
       grid fine
end
task dft optimize

and

title "benzophenone, neutral, gas phase, frequency"
start benzophenone_gas_freq
charge 0
geometry
       ...xyz coordinates of optimised structure...
end
dft
       xc b3lyp
       mult 1
       grid fine
end
task dft freq

and

title "benzophenone, neutral, COSMO"
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

and

title "Benzophenone, anion, gas phase"
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

etc.

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 cycles
             Convg  =    0.6297D-08             -V/T =  2.0095

Anion gas frequency:

 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
G_{anion,gasphase}=-576.656219870+0.150074=-576.506145870
ε_{solvation,anion}=-576.732948458-(-576.656219870)=-0.076728588

ε_{0,neutral,gasphase}=-576.648098680
G_{neutral,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.648088887542
      One electron energy =    -2313.472465671473
           Coulomb energy =     1046.450004818160
    Exchange-Corr. energy =      -82.832342635322
 Nuclear repulsion energy =      773.206714601093

neutral, gas, freq

 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)

Anion gas

       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

G_{neutral,gasphase}=-576.648088887542+0.15456920697551748261+0.191895=-576.49351968056648251739+0.191895
ε_{solvation,neutral}=0.0122412722

Gcorr=(125.486-298.15* 102.152/1000)=0.15141494338035305421

ε_{0,anion,gasphase}= -576.656222314875
G_{anion,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. 

However, not including solvation at all yields values which are completely off: -3.85 V vs SCE for gaussian and 3.88 for nwchem. The solvation free energy correction is the largest factor in yielding a 'correct' (as in order of magnitude) result -- and if I'd be forced to venture a guess, solvation models/implementations (and their empirical input) both differ a lot between modelling packages and have improved a lot in the intervening 10 year.

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...]