Calculation of Redox Potentials¶
In this tutorial, we use ORCA together with the ORCA Python interface (OPI) to calculate the redox potentials of anisole and dimethoxybenzene in acetonitrile (MeCN) using density functional theory (DFT) with implicit solvation. The redox potentials are obtained from the Gibbs free energy difference between the neutral and oxidized species in solution, following a standard thermodynamic cycle. The calculations required to obtain these free energy differences include geometry optimizations and frequency calculations for both the oxidized and neutral species, as well as high-level single-point calculations using a hybrid DFT functional on the optimized structures. With this approach, four calculations are required for each redox couple: geometry optimization and single-point energy evaluation for both the oxidized and neutral states. Setting these calculations up, running them, and parsing the output can be simplified with OPI as shown below.
Step 1: Installing 3rd-party Dependencies¶
import sys
!{sys.executable} -m pip install --upgrade py3Dmol
Step 1: Import Dependencies¶
We begin by importing all required Python modules for this tutorial. These include:
OPI: for setting up, running, and parsing ORCA quantum chemistry calculations.
py3Dmol: for interactive 3D visualization of molecular structures directly in the notebook.
Note: We additionally import modules for visualization/plotting like
py3Dmol. For this, it might be necessary to installpy3Dmolinto your OPIvenv(e.g., by activating the.venvand usinguv pip install py3Dmol).
# > Import pathlib for directory handling
from pathlib import Path
import shutil
# > OPI imports for performing ORCA calculations and reading output
from opi.core import Calculator
from opi.input.structures.structure import Structure
from opi.input.simple_keywords import Dft, Task, SolvationModel, Solvent, BasisSet
from opi.input.simple_keywords import SimpleKeyword
from opi.input.blocks import Block, BlockScf
from opi.output.core import Output
from opi.utils.units import AU_TO_EV
# > Import py3Dmol for 3D molecular visualization
import py3Dmol
Step 2: Define Working Directory¶
All actual calculations will be performed in a subfolder redox.
# > Calculation is performed in `fock_matrix_diagonalization`
working_dir = Path("redox")
# > The `working_dir`is automatically (re-)created
shutil.rmtree(working_dir, ignore_errors=True)
working_dir.mkdir()
Step 3: Setting up the Structures¶
We define the input structures as xyz-files:
# > define cartesian coordinates as python strings
xyz_data_dimethoxybenzene="""
20
C -1.90530 2.31127 0.12469
C -3.12410 1.64361 0.01546
O 0.48884 -0.42420 0.16516
H -2.60202 5.51597 0.19386
H -3.42166 4.33208 -0.83191
C -3.15240 0.24805 -0.04572
H -4.06815 2.17619 -0.02342
C -1.96561 -0.48783 0.00114
H -4.10509 -0.26920 -0.13016
C -0.74187 0.17066 0.11006
H -2.03239 -1.56917 -0.04930
C -0.72242 1.56733 0.17175
O -1.73499 3.66676 0.19535
H -3.57809 4.26171 0.96693
C -2.90794 4.46750 0.12570
H 0.23029 2.08457 0.25726
C 0.52053 -1.84452 0.10448
H 1.56753 -2.15857 0.15754
H -0.00347 -2.28803 0.95775
H 0.11176 -2.20889 -0.84395
"""
xyz_data_anisole ="""
16
C -1.78955 2.19474 0.12638
C -2.99242 1.49535 0.01839
H 0.35762 -0.45367 0.15060
H -2.56504 5.37952 0.21701
H -3.33324 4.19503 -0.84863
C -2.98962 0.09748 -0.03760
H -3.94916 2.00459 -0.02068
C -1.78685 -0.60529 0.00944
H -3.92944 -0.44263 -0.11848
C -0.58360 0.08793 0.11419
H -1.78789 -1.69118 -0.03394
C -0.58764 1.48247 0.17211
O -1.65102 3.55370 0.19785
H -3.52457 4.08734 0.94531
C -2.84294 4.32528 0.12189
H 0.35314 2.02072 0.25508
"""
# > Visualize the input structure - dimethoxybenzene
view = py3Dmol.view(width=400, height=400)
view.addModel(xyz_data_dimethoxybenzene, 'xyz')
view.setStyle({}, {'stick': {}, 'sphere': {'scale': 0.3}})
view.zoomTo()
view.show()
# > Read the structure into OPI structure
structure_dimethoxybenzene = Structure.from_xyz_block(xyz_data_dimethoxybenzene)
# > Visualize the input structure - anisole
view = py3Dmol.view(width=400, height=400)
view.addModel(xyz_data_anisole, 'xyz')
view.setStyle({}, {'stick': {}, 'sphere': {'scale': 0.3}})
view.zoomTo()
view.show()
# > Read the structure into OPI structure
structure_anisole = Structure.from_xyz_block(xyz_data_anisole)
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
Step 4: Defining and Running the ORCA Calculations¶
First, we define a general function run_calc for handing simple keywords to an opi calculator, running the calculation, and retrieving the output.
def run_calc(basename : str, working_dir: Path, sk_list: list[SimpleKeyword], structure: Structure, block_list: list[Block] = [], ncores: int = 4) -> Output:
"""Perform an ORCA calculation and get the output"""
print(f"Running the calculation {basename} ... ", end="")
# > Set up a Calculator object, the basename for the ORCA calculation is also set
calc = Calculator(basename=basename, working_dir=working_dir)
# > Assign structure to calculator
calc.structure = structure
# > Use simple keywords in calculator
calc.input.add_simple_keywords(*sk_list)
if block_list:
for block in block_list:
calc.input.add_blocks(block)
# > Define number of CPUs for the calcualtion
calc.input.ncores = ncores # > CPUs for this ORCA run (default: 1)
# > perform the calculation
calc.write_and_run()
print("Done")
# > get the output
output = calc.get_output()
# > Check for normal termination, SCF, and if requested geometry optimization
if not output.terminated_normally():
raise RuntimeError(f"ORCA calculation failed! Path: {working_dir}/{basename}")
if not output.scf_converged():
raise RuntimeError(f"ORCA SCF did not converge! Path: {working_dir}/{basename}!")
if Task.OPT in sk_list and not output.geometry_optimization_converged():
raise RuntimeError(f"ORCA geometry optimization did not converge! Path: {working_dir}/{basename}!")
# > Parse the output
output.parse()
# > return the output
return output
Next, we use the predefined function run_calc to define the function run_redox for the DFT calculations of both species of a redox couple.
Here, we also define a level of theory:
Geometry optimization and frequency calculations are performed with the composite r²SCAN-3c method together with the implicit solvation model SMD.
On the optimized structures high-level single-points are performed with the hybrid DFT method ωB97X-V/def2-TZVP, also with the SMD model.
def run_redox(basename : str, working_dir: Path, structure: Structure, solvent: Solvent, ncores: int = 4) -> tuple[Output, Output, Output, Output]:
"""Run a geometry optimization and a single-point calculation for both charge states of a redox couple and return the OPI outputs for all calculations"""
# Define a level of theory for geometry optimization (r²SCAN-3c)
opt_sk_list = [
Dft.R2SCAN_3C, # > r²SCAN-3c method (Comes with a predefined basis set)
Task.OPT, # > Perform the geometry optimization
Task.FREQ, # > After geometry optimization perform a frequency calculation
SolvationModel.SMD(solvent) # SMD solvation model
# > Note that the order of OPT and FREQ does not matter! FREQ is always performed after OPT!
]
# Defin a level of theory for single-point calculation (wB97X-V/def2-TZVP)
sp_sk_list = [
Dft.WB97X_V, # > wB97X-V
BasisSet.DEF2_TZVP,
SolvationModel.SMD(solvent), # SMD solvation model
Task.SP
]
# SCF stability performance block for single-point calculations
block_list = [ BlockScf(stabperform=True) ]
# Optimize the structures
# > Neutral species
opt_neutral = run_calc(basename=f"{basename}_neutral_opt", working_dir=working_dir, sk_list=opt_sk_list, structure=structure, ncores=ncores)
# > Charged species
# > Adjust the charge
structure.charge += 1
# > Adjust the electron-spin multiplicity
structure.set_ls_multiplicity() # equal to structure.multiplicity = 2
opt_cation = run_calc(basename=f"{basename}_cation_opt", working_dir=working_dir, sk_list=opt_sk_list, structure=structure, ncores=ncores)
# Perform high-level single-point calculations
# > Neutral
sp_neutral = run_calc(basename=f"{basename}_neutral_sp", working_dir=working_dir, sk_list=sp_sk_list, structure=opt_neutral.get_structure(), block_list=block_list, ncores=ncores)
# > Charged
sp_cation = run_calc(basename=f"{basename}_cation_sp", working_dir=working_dir, sk_list=sp_sk_list, structure=opt_cation.get_structure(), block_list=block_list, ncores=ncores)
# > return outputs
return opt_cation, sp_cation, opt_neutral, sp_neutral
The function run_redox is then applied to dimethoxybenzene and anisole and the actual calculations are performed:
dimethoxybenzene_opt_cation, dimethoxybenzene_sp_cation, dimethoxybenzene_opt_neutral, dimethoxybenzene_sp_neutral = run_redox("dimethoxybenzene", working_dir=working_dir, structure=structure_dimethoxybenzene, solvent=Solvent.ACETONITRILE)
anisole_opt_cation, anisole_sp_cation, anisole_opt_neutral, anisole_sp_neutral = run_redox("anisole", working_dir=working_dir, structure=structure_anisole, solvent=Solvent.ACETONITRILE)
Running the calculation dimethoxybenzene_neutral_opt ... Done
Running the calculation dimethoxybenzene_cation_opt ... Done
Running the calculation dimethoxybenzene_neutral_sp ... Done
Running the calculation dimethoxybenzene_cation_sp ... Done
Running the calculation anisole_neutral_opt ... Done
Running the calculation anisole_cation_opt ... Done
Running the calculation anisole_neutral_sp ... Done
Running the calculation anisole_cation_sp ... Done
Step 5: Process the Results¶
After running the DFT calculations, the resulting free energies are obtained from the outputs. Electronic energies are taken from the single-point calculations using the function get_final_energy, while thermostatistical corrections are obtained from the frequency calculations in the geometry optimization outputs via the function get_final_energy_delta. These contributions are combined to yield a final Gibbs free energy for each species.
The calculated absolute free energies are then inserted into the Nernst equation to obtain standard redox potentials:
Here, \(\Delta G^\circ_{\mathrm{R/O}}\) denotes the Gibbs free energy difference of the neutral from the oxidized state, i.e., \(\Delta G^\circ_{\mathrm{R/O}} = G^\circ_{\mathrm{neutral}} - G^\circ_{\mathrm{oxidized}}\). Note that all redox processes considered here involve one-electron transfer, so \(n = 1\). The Faraday constant \(F\) is absorbed into the unit conversion factor used to convert the free energies from atomic units into volts.
To report meaningful redox potentials, a reference electrode (\(E^\circ_{\mathrm{ABS}}(\mathrm{REF})\)) is required. Typically, there are two approaches:
Take a pre-computed shift from the literature, e.g., 4.422 V for SCE in MeCN (Synlett 2016, 27, 714-723)
Or use an already known experimental redox potential, calculate its absolute redox potential, and use this as a reference.
Below, we first apply an absolute reference shift and then determine the redox potential of anisole relative to dimethoxybenzene.
# Absolute free energy difference (oxidation) in eV:
# ΔG_ox = G(cation) - G(neutral)
# For a one-electron couple, E_red(abs) = ΔG_ox (in V) because E_red = -ΔG_red and ΔG_red = -ΔG_ox.
SCE_ABS_MECN = 4.422 # V, Synlett 2016, 27, 714-723
# --- 1,3-dimethoxybenzene ---
G_ox = dimethoxybenzene_sp_cation.get_final_energy() + dimethoxybenzene_opt_cation.get_free_energy_delta()
G_red = dimethoxybenzene_sp_neutral.get_final_energy() + dimethoxybenzene_opt_neutral.get_free_energy_delta()
delta_g_ox_ev_dim = (G_ox - G_red) * AU_TO_EV # oxidation free energy in eV
E_red_abs_dim = delta_g_ox_ev_dim # for 1e-, because E = -ΔG_red and ΔG_red = -ΔG_ox
E_red_vs_SCE_dim = E_red_abs_dim - SCE_ABS_MECN # substract literature shift
print(f"The calculated reduction potential (vs SCE, literature shift) of 1,3-dimethoxybenzene is {E_red_vs_SCE_dim:.2f} V")
# --- anisole ---
G_ox = anisole_sp_cation.get_final_energy() + anisole_opt_cation.get_free_energy_delta()
G_red = anisole_sp_neutral.get_final_energy() + anisole_opt_neutral.get_free_energy_delta()
delta_g_ox_ev_anis = (G_ox - G_red) * AU_TO_EV # oxidation free energy in eV
E_red_abs_anis = delta_g_ox_ev_anis # for 1e-, because E = -ΔG_red and ΔG_red = -ΔG_ox
E_red_vs_SCE_anis = E_red_abs_anis - SCE_ABS_MECN # substract literature shift
print(f"The calculated reduction potential (vs SCE, literature shift) of anisole is {E_red_vs_SCE_anis:.2f} V")
# --- anisole vs dimethoxybenzene using experimental anchor ---
E_exp_dim_vs_SCE = 1.549 # V, J. Org. Chem. 2014, 79, 9297-9304
E_calc_anis_vs_SCE = (E_red_abs_anis - E_red_abs_dim) + E_exp_dim_vs_SCE
print(f"The calculated reduction potential (vs SCE, referenced to exp. dimethoxybenzene) for anisole is {E_calc_anis_vs_SCE:.2f} V")
# --- Experimental references ---
print("Experimental E°(anisole) vs SCE: 1.77 V in MeCN (J. Org. Chem. 2014, 79, 9297-9304)")
print("Experimental E°(1,3-dimethoxybenzene) vs SCE: 1.55 V in MeCN (J. Org. Chem. 2014, 79, 9297-9304)")
The calculated reduction potential (vs SCE, literature shift) of 1,3-dimethoxybenzene is 1.38 V
The calculated reduction potential (vs SCE, literature shift) of anisole is 1.45 V
The calculated reduction potential (vs SCE, referenced to exp. dimethoxybenzene) for anisole is 1.62 V
Experimental E°(anisole) vs SCE: 1.77 V (J. Org. Chem. 2014, 79, 9297–9304)
Experimental E°(1,3-dimethoxybenzene) vs SCE: 1.55 V (J. Org. Chem. 2014, 79, 9297–9304)
One can observe that the agreement with experiment is improved when using a calculated relative redox potential anchored to an experimentally known reference compound. This approach benefits from error cancellation between similar molecules, as systematic errors arising from the electronic structure method and solvation model cancel when taking free energy differences. Such error compensation is most effective when the reference compound is similar to the target system.
Step 7: Visualization¶
We can visualize molecular orbitals (MOs) or electron densities. For example, we can plot the HOMOs of 1,3-dimethoxybenzene and anisole to gain insight into where electron density is removed upon oxidation:
# > to keep the size of this notebook small the resolution is not too large
def view_homo(output: Output):
homo_id = output.get_homo().index
print(output.basename)
print(f"The id of the HOMO is {homo_id}")
cube_output = output.plot_mo(homo_id,resolution=40)
cube_data = cube_output.cube
view = py3Dmol.view(width=500, height=500)
view.addModel(output.get_structure().to_xyz_block(), "xyz")
view.setStyle({'stick': {'radius': 0.1}, 'sphere': {'scale': 0.2}})
view.addVolumetricData(cube_data, "cube", {"isoval": 0.05, "color": "blue", "opacity": 0.8})
view.addVolumetricData(cube_data, "cube", {"isoval": -0.05, "color": "red", "opacity": 0.8})
view.zoomTo()
view.show()
view_homo(anisole_sp_neutral)
view_homo(dimethoxybenzene_sp_neutral)
anisole_neutral_sp
The id of the HOMO is 28
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
dimethoxybenzene_neutral_sp
The id of the HOMO is 36
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
One can see that the HOMOs of both 1,3-dimethoxybenzene and anisole are delocalized over the aromatic \(\pi\)-system and show significant contributions from the methoxy substituents. This indicates that upon oxidation, electron density is primarily removed from the \(\pi\)-system in both molecules. The presence of an additional methoxy group in 1,3-dimethoxybenzene allows for increased delocalization and stabilization of the resulting radical cation, which is consistent with its lower (less positive) redox potential compared to anisole, i.e., 1,3-dimethoxybenzene is easier to oxidize.
Summary¶
In this notebook, we demonstrated how to calculate redox potentials using OPI. To obtain meaningful redox potentials, a reference shift is required to relate the calculated free energy differences to an experimental scale. Such reference values can either be taken from the literature or derived by using an experimentally known redox couple as a reference, which is explicitly calculated alongside the target system. This relative referencing approach benefits from error cancellation and is commonly employed using the ferrocene/ferrocenium redox couple.