{ "cells": [ { "cell_type": "markdown", "id": "71a2f48b", "metadata": {}, "source": [ "# Calculation of Redox Potentials\n", "\n", "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." ] }, { "cell_type": "markdown", "id": "d3da3e2e-870f-46e2-b1cc-eacfd4bc9f38", "metadata": {}, "source": [ "## Step 1: Installing 3rd-party Dependencies" ] }, { "cell_type": "code", "execution_count": null, "id": "9b2f99e1-beab-4ea4-8749-722c553b4329", "metadata": {}, "outputs": [], "source": [ "import sys\n", "!{sys.executable} -m pip install --upgrade py3Dmol" ] }, { "cell_type": "markdown", "id": "e86ab3fc", "metadata": {}, "source": [ "## Step 1: Import Dependencies\n", "\n", "We begin by importing all required Python modules for this tutorial. These include:\n", "\n", "- **OPI**: for setting up, running, and parsing ORCA quantum chemistry calculations.\n", "- **py3Dmol**: for interactive 3D visualization of molecular structures directly in the notebook.\n", "\n", "> **Note:** We additionally import modules for visualization/plotting like `py3Dmol`. For this, it might be necessary to install `py3Dmol` into your OPI `venv` (e.g., by activating the `.venv` and using `uv pip install py3Dmol`)." ] }, { "cell_type": "code", "execution_count": null, "id": "b4b88031", "metadata": {}, "outputs": [], "source": [ "# > Import pathlib for directory handling\n", "from pathlib import Path\n", "import shutil\n", "\n", "# > OPI imports for performing ORCA calculations and reading output\n", "from opi.core import Calculator\n", "from opi.input.structures.structure import Structure\n", "from opi.input.simple_keywords import Dft, Task, SolvationModel, Solvent, BasisSet\n", "from opi.input.simple_keywords import SimpleKeyword\n", "from opi.input.blocks import Block, BlockScf\n", "from opi.output.core import Output\n", "from opi.utils.units import AU_TO_EV\n", "\n", "# > Import py3Dmol for 3D molecular visualization\n", "import py3Dmol" ] }, { "cell_type": "markdown", "id": "62980249", "metadata": {}, "source": [ "## Step 2: Define Working Directory\n", "\n", "All actual calculations will be performed in a subfolder `redox`." ] }, { "cell_type": "code", "execution_count": null, "id": "8712d717", "metadata": {}, "outputs": [], "source": [ "# > Calculation is performed in `fock_matrix_diagonalization`\n", "working_dir = Path(\"redox\")\n", "# > The `working_dir`is automatically (re-)created\n", "shutil.rmtree(working_dir, ignore_errors=True)\n", "working_dir.mkdir()" ] }, { "cell_type": "markdown", "id": "986bd3db", "metadata": {}, "source": [ "## Step 3: Setting up the Structures\n", "\n", "We define the input structures as xyz-files:" ] }, { "cell_type": "code", "execution_count": 3, "id": "3613f8e7", "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", "text/html": [ "
\n", "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", "text/html": [ "
\n", "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# > define cartesian coordinates as python strings\n", "xyz_data_dimethoxybenzene=\"\"\"\n", "20\n", "\n", "C -1.90530 2.31127 0.12469\n", "C -3.12410 1.64361 0.01546\n", "O 0.48884 -0.42420 0.16516\n", "H -2.60202 5.51597 0.19386\n", "H -3.42166 4.33208 -0.83191\n", "C -3.15240 0.24805 -0.04572\n", "H -4.06815 2.17619 -0.02342\n", "C -1.96561 -0.48783 0.00114\n", "H -4.10509 -0.26920 -0.13016\n", "C -0.74187 0.17066 0.11006\n", "H -2.03239 -1.56917 -0.04930\n", "C -0.72242 1.56733 0.17175\n", "O -1.73499 3.66676 0.19535\n", "H -3.57809 4.26171 0.96693\n", "C -2.90794 4.46750 0.12570\n", "H 0.23029 2.08457 0.25726\n", "C 0.52053 -1.84452 0.10448\n", "H 1.56753 -2.15857 0.15754\n", "H -0.00347 -2.28803 0.95775\n", "H 0.11176 -2.20889 -0.84395\n", "\"\"\"\n", "\n", "xyz_data_anisole =\"\"\"\n", "16\n", "\n", "C -1.78955 2.19474 0.12638\n", "C -2.99242 1.49535 0.01839\n", "H 0.35762 -0.45367 0.15060\n", "H -2.56504 5.37952 0.21701\n", "H -3.33324 4.19503 -0.84863\n", "C -2.98962 0.09748 -0.03760\n", "H -3.94916 2.00459 -0.02068\n", "C -1.78685 -0.60529 0.00944\n", "H -3.92944 -0.44263 -0.11848\n", "C -0.58360 0.08793 0.11419\n", "H -1.78789 -1.69118 -0.03394\n", "C -0.58764 1.48247 0.17211\n", "O -1.65102 3.55370 0.19785\n", "H -3.52457 4.08734 0.94531\n", "C -2.84294 4.32528 0.12189\n", "H 0.35314 2.02072 0.25508\n", "\"\"\"\n", "\n", "# > Visualize the input structure - dimethoxybenzene\n", "view = py3Dmol.view(width=400, height=400)\n", "view.addModel(xyz_data_dimethoxybenzene, 'xyz')\n", "view.setStyle({}, {'stick': {}, 'sphere': {'scale': 0.3}})\n", "view.zoomTo()\n", "view.show()\n", "\n", "# > Read the structure into OPI structure\n", "structure_dimethoxybenzene = Structure.from_xyz_block(xyz_data_dimethoxybenzene)\n", "\n", "# > Visualize the input structure - anisole\n", "view = py3Dmol.view(width=400, height=400)\n", "view.addModel(xyz_data_anisole, 'xyz')\n", "view.setStyle({}, {'stick': {}, 'sphere': {'scale': 0.3}})\n", "view.zoomTo()\n", "view.show()\n", "\n", "# > Read the structure into OPI structure\n", "structure_anisole = Structure.from_xyz_block(xyz_data_anisole)" ] }, { "cell_type": "markdown", "id": "5154dda2", "metadata": {}, "source": [ "## Step 4: Defining and Running the ORCA Calculations\n", "\n", "First, we define a general function `run_calc` for handing simple keywords to an opi calculator, running the calculation, and retrieving the output." ] }, { "cell_type": "code", "execution_count": null, "id": "6d5e55e6", "metadata": {}, "outputs": [], "source": [ "def run_calc(basename : str, working_dir: Path, sk_list: list[SimpleKeyword], structure: Structure, block_list: list[Block] = [], ncores: int = 4) -> Output:\n", " \"\"\"Perform an ORCA calculation and get the output\"\"\"\n", " print(f\"Running the calculation {basename} ... \", end=\"\")\n", " # > Set up a Calculator object, the basename for the ORCA calculation is also set\n", " calc = Calculator(basename=basename, working_dir=working_dir)\n", " # > Assign structure to calculator\n", " calc.structure = structure\n", "\n", " # > Use simple keywords in calculator\n", " calc.input.add_simple_keywords(*sk_list)\n", "\n", " if block_list:\n", " for block in block_list:\n", " calc.input.add_blocks(block)\n", "\n", " # > Define number of CPUs for the calcualtion\n", " calc.input.ncores = ncores # > CPUs for this ORCA run (default: 1)\n", "\n", " # > perform the calculation\n", " calc.write_and_run()\n", " print(\"Done\")\n", "\n", " # > get the output\n", " output = calc.get_output()\n", "\n", " # > Check for normal termination, SCF, and if requested geometry optimization\n", " if not output.terminated_normally():\n", " raise RuntimeError(f\"ORCA calculation failed! Path: {working_dir}/{basename}\")\n", " if not output.scf_converged():\n", " raise RuntimeError(f\"ORCA SCF did not converge! Path: {working_dir}/{basename}!\")\n", " if Task.OPT in sk_list and not output.geometry_optimization_converged():\n", " raise RuntimeError(f\"ORCA geometry optimization did not converge! Path: {working_dir}/{basename}!\")\n", " \n", " # > Parse the output \n", " output.parse()\n", "\n", " # > return the output\n", " return output" ] }, { "cell_type": "markdown", "id": "95616a66", "metadata": {}, "source": [ "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. \n", "Here, we also define a level of theory: \n", "- Geometry optimization and frequency calculations are performed with the composite r²SCAN-3c method together with the implicit solvation model SMD. \n", "- On the optimized structures high-level single-points are performed with the hybrid DFT method ωB97X-V/def2-TZVP, also with the SMD model." ] }, { "cell_type": "code", "execution_count": 5, "id": "9d612e37", "metadata": {}, "outputs": [], "source": [ "def run_redox(basename : str, working_dir: Path, structure: Structure, solvent: Solvent, ncores: int = 4) -> tuple[Output, Output, Output, Output]:\n", " \"\"\"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\"\"\"\n", " # Define a level of theory for geometry optimization (r²SCAN-3c)\n", " opt_sk_list = [\n", " Dft.R2SCAN_3C, # > r²SCAN-3c method (Comes with a predefined basis set)\n", " Task.OPT, # > Perform the geometry optimization\n", " Task.FREQ, # > After geometry optimization perform a frequency calculation\n", " SolvationModel.SMD(solvent) # SMD solvation model\n", " # > Note that the order of OPT and FREQ does not matter! FREQ is always performed after OPT! \n", " ]\n", "\n", " # Defin a level of theory for single-point calculation (wB97X-V/def2-TZVP)\n", " sp_sk_list = [\n", " Dft.WB97X_V, # > wB97X-V\n", " BasisSet.DEF2_TZVP,\n", " SolvationModel.SMD(solvent), # SMD solvation model\n", " Task.SP\n", " ]\n", "\n", " # SCF stability performance block for single-point calculations\n", " block_list = [ BlockScf(stabperform=True) ]\n", "\n", " # Optimize the structures\n", "\n", " # > Neutral species\n", " opt_neutral = run_calc(basename=f\"{basename}_neutral_opt\", working_dir=working_dir, sk_list=opt_sk_list, structure=structure, ncores=ncores)\n", "\n", " # > Charged species\n", " # > Adjust the charge\n", " structure.charge += 1\n", " # > Adjust the electron-spin multiplicity\n", " structure.set_ls_multiplicity() # equal to structure.multiplicity = 2\n", " opt_cation = run_calc(basename=f\"{basename}_cation_opt\", working_dir=working_dir, sk_list=opt_sk_list, structure=structure, ncores=ncores)\n", "\n", " # Perform high-level single-point calculations\n", "\n", " # > Neutral\n", " 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)\n", "\n", " # > Charged\n", " 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)\n", "\n", " # > return outputs \n", " return opt_cation, sp_cation, opt_neutral, sp_neutral" ] }, { "cell_type": "markdown", "id": "183c60f3", "metadata": {}, "source": [ "The function `run_redox` is then applied to dimethoxybenzene and anisole and the actual calculations are performed:" ] }, { "cell_type": "code", "execution_count": 6, "id": "8282dff7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Running the calculation dimethoxybenzene_neutral_opt ... Done\n", "Running the calculation dimethoxybenzene_cation_opt ... Done\n", "Running the calculation dimethoxybenzene_neutral_sp ... Done\n", "Running the calculation dimethoxybenzene_cation_sp ... Done\n", "Running the calculation anisole_neutral_opt ... Done\n", "Running the calculation anisole_cation_opt ... Done\n", "Running the calculation anisole_neutral_sp ... Done\n", "Running the calculation anisole_cation_sp ... Done\n" ] } ], "source": [ "\n", "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)\n", "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)" ] }, { "cell_type": "markdown", "id": "39f3334e", "metadata": {}, "source": [ "## Step 5: Process the Results\n", "\n", "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.\n", "\n", "The calculated absolute free energies are then inserted into the Nernst equation to obtain standard redox potentials:\n", "\n", "$$\n", "E^\\circ = -\\frac{\\Delta G^\\circ_{\\mathrm{R/O}}}{nF} - E^\\circ_{\\mathrm{ABS}}(\\mathrm{REF})\n", "$$\n", "\n", "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.\n", "\n", "To report meaningful redox potentials, a reference electrode ($E^\\circ_{\\mathrm{ABS}}(\\mathrm{REF})$) is required. Typically, there are two approaches:\n", "1. Take a pre-computed shift from the literature, e.g., 4.422 V for SCE in MeCN (Synlett 2016, 27, 714-723)\n", "2. Or use an already known experimental redox potential, calculate its absolute redox potential, and use this as a reference. \n", "\n", "Below, we first apply an absolute reference shift and then determine the redox potential of anisole relative to dimethoxybenzene.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "abcff9b1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The calculated reduction potential (vs SCE, literature shift) of 1,3-dimethoxybenzene is 1.38 V\n", "The calculated reduction potential (vs SCE, literature shift) of anisole is 1.45 V\n", "The calculated reduction potential (vs SCE, referenced to exp. dimethoxybenzene) for anisole is 1.62 V\n", "Experimental E°(anisole) vs SCE: 1.77 V (J. Org. Chem. 2014, 79, 9297–9304)\n", "Experimental E°(1,3-dimethoxybenzene) vs SCE: 1.55 V (J. Org. Chem. 2014, 79, 9297–9304)\n" ] } ], "source": [ "# Absolute free energy difference (oxidation) in eV:\n", "# ΔG_ox = G(cation) - G(neutral)\n", "# For a one-electron couple, E_red(abs) = ΔG_ox (in V) because E_red = -ΔG_red and ΔG_red = -ΔG_ox.\n", "\n", "SCE_ABS_MECN = 4.422 # V, Synlett 2016, 27, 714-723\n", "\n", "# --- 1,3-dimethoxybenzene ---\n", "G_ox = dimethoxybenzene_sp_cation.get_final_energy() + dimethoxybenzene_opt_cation.get_free_energy_delta() \n", "G_red = dimethoxybenzene_sp_neutral.get_final_energy() + dimethoxybenzene_opt_neutral.get_free_energy_delta() \n", "\n", "delta_g_ox_ev_dim = (G_ox - G_red) * AU_TO_EV # oxidation free energy in eV\n", "E_red_abs_dim = delta_g_ox_ev_dim # for 1e-, because E = -ΔG_red and ΔG_red = -ΔG_ox\n", "E_red_vs_SCE_dim = E_red_abs_dim - SCE_ABS_MECN # substract literature shift\n", "\n", "print(f\"The calculated reduction potential (vs SCE, literature shift) of 1,3-dimethoxybenzene is {E_red_vs_SCE_dim:.2f} V\")\n", "\n", "# --- anisole ---\n", "G_ox = anisole_sp_cation.get_final_energy() + anisole_opt_cation.get_free_energy_delta()\n", "G_red = anisole_sp_neutral.get_final_energy() + anisole_opt_neutral.get_free_energy_delta()\n", "\n", "delta_g_ox_ev_anis = (G_ox - G_red) * AU_TO_EV # oxidation free energy in eV\n", "E_red_abs_anis = delta_g_ox_ev_anis # for 1e-, because E = -ΔG_red and ΔG_red = -ΔG_ox\n", "E_red_vs_SCE_anis = E_red_abs_anis - SCE_ABS_MECN # substract literature shift\n", "\n", "print(f\"The calculated reduction potential (vs SCE, literature shift) of anisole is {E_red_vs_SCE_anis:.2f} V\")\n", "\n", "# --- anisole vs dimethoxybenzene using experimental anchor ---\n", "E_exp_dim_vs_SCE = 1.549 # V, J. Org. Chem. 2014, 79, 9297-9304\n", "E_calc_anis_vs_SCE = (E_red_abs_anis - E_red_abs_dim) + E_exp_dim_vs_SCE\n", "\n", "print(f\"The calculated reduction potential (vs SCE, referenced to exp. dimethoxybenzene) for anisole is {E_calc_anis_vs_SCE:.2f} V\")\n", "\n", "# --- Experimental references ---\n", "print(\"Experimental E°(anisole) vs SCE: 1.77 V in MeCN (J. Org. Chem. 2014, 79, 9297-9304)\")\n", "print(\"Experimental E°(1,3-dimethoxybenzene) vs SCE: 1.55 V in MeCN (J. Org. Chem. 2014, 79, 9297-9304)\")\n" ] }, { "cell_type": "markdown", "id": "cd52930f", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "id": "68d4271e", "metadata": {}, "source": [ "## Step 7: Visualization \n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 14, "id": "bce13dff", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "anisole_neutral_sp\n", "The id of the HOMO is 28\n" ] }, { "data": { "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", "text/html": [ "
\n", "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "dimethoxybenzene_neutral_sp\n", "The id of the HOMO is 36\n" ] }, { "data": { "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", "text/html": [ "
\n", "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# > to keep the size of this notebook small the resolution is not too large\n", "def view_homo(output: Output):\n", " homo_id = output.get_homo().index\n", " print(output.basename)\n", " print(f\"The id of the HOMO is {homo_id}\")\n", " cube_output = output.plot_mo(homo_id,resolution=40)\n", " cube_data = cube_output.cube\n", "\n", " view = py3Dmol.view(width=500, height=500)\n", " view.addModel(output.get_structure().to_xyz_block(), \"xyz\")\n", " view.setStyle({'stick': {'radius': 0.1}, 'sphere': {'scale': 0.2}})\n", " view.addVolumetricData(cube_data, \"cube\", {\"isoval\": 0.05, \"color\": \"blue\", \"opacity\": 0.8})\n", " view.addVolumetricData(cube_data, \"cube\", {\"isoval\": -0.05, \"color\": \"red\", \"opacity\": 0.8})\n", " view.zoomTo()\n", " view.show()\n", "\n", "view_homo(anisole_sp_neutral)\n", "view_homo(dimethoxybenzene_sp_neutral)" ] }, { "cell_type": "markdown", "id": "c3524b26", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "id": "6fd66184", "metadata": {}, "source": [ "## Summary\n", "\n", "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. " ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.3" } }, "nbformat": 4, "nbformat_minor": 5 }