Reduction and Oxidation Potentials
**********************************

ADF Redox Technical Batteries

Calculate one-electron oxidation and reduction potentials with three
PLAMS redox workflows: a fast screening method, a direct implicit-
solvent method, and a thermodynamic-cycle method.


Definitions and introduction
============================

**Reduction potential**: A + e- -> A- ; ΔG0 = -nFE0

**Oxidation potential**: A -> A+ + e- ; ΔG0 = nFE0

There are three PLAMS recipes for calculating one-electron reduction
or oxidation potentials in implicit solvent:

* "AMSRedoxDirectJob": The best method. Geometry optimizations (and
  optionally frequencies) are calculated for both neutral and
  reduced/oxidized species in implicit solvent. The solvent must be
  supported by ADF with the COSMO solvation method. Requires an ADF
  license.

* "AMSRedoxThermodynamicCycleJob": Only useful if you include the
  vibrations (frequencies) and the molecule is large (in which case it
  is faster but less accurate than AMSRedoxDirectJob). The frequencies
  are only calculated for the gasphase molecule. A thermodynamic cycle
  gives the reduction or oxidation potential. The solvent must be
  supported by ADF with the COSMO solvation method. Requires an ADF
  license.

* "AMSRedoxScreeningJob": The fastest (and least accurate) method.
  Geometry optimizations are performed at the GFN1-xTB level of
  theory. The solvation free energy is evaluated by COSMO-RS.
  Vibrational effects are always implicitly accounted for. A ".coskf"
  file for the solvent is required - this can either be obtained from
  the ADFCRS-2018 database or generated with an
  "ADFCOSMORSCompoundJob". Requires DFTB, ADF, and COSMO-RS licenses.

The "AMSRedoxScreeningJob" workflow was developed by Belic et al.,
Phys. Chem. Chem. Phys. 24, 197–210 (2022) for calculating oxidation
potentials. The paper also describes the other workflows in detail.

Note:

  The "AMSRedoxScreeningJob" uses the recommended ADF settings for
  generating .coskf files. Belic et al. used a different density
  functional. Using the class will give slightly different results
  compared to Belic et al.

The **free energy of the electron** is set to be –0.0375 eV (see Belic
et al.).

In this example the above three workflows are used to evaluate the
**reduction potential of benzoquinone in water**. The experimental
value is E⁰ = +0.10 V relative to SHE (**standard hydrogen
electrode**).

Ho et al.  (preprint pdf) evaluated E⁰ = -0.40 V or -0.28 V for
benzoquinone relative to SHE  using different computational methods.

The reduction and oxidation potentials calculated by the PLAMS classes
are given on an **absolute scale**. On this scale, the SHE is at +4.42
or +4.28 V (see Ho et al.). To get potentials relative to SHE, we
therefore subtract 4.42 V from the calculated values.

Note:

  Energy differences in eV correspond directly to potentials in V.

Tip:

  If you compare many calculated reduction potentials, use one of them
  as the reference state. See for example Huyhn et al., J. Am. Chem.
  Soc. 2016, 138, 49, 15903–15910 in the case of substituted quinones.

Important:

  These classes can only be used for

  * molecules that do not undergo big conformational changes (or
    dissociation) upon oxidation/reduction.

  * one-electron reduction/oxidation. For multiple electrons, run the
    same script with different initial charges for the molecule

  * molecules with 0 or 1 unpaired electrons

**Requires:** AMS2026 or later

-[ Related examples ]-

* Tuning the Range Separation Parameter in ADF

* Charge Transfer Integrals with ADF

* Reorganization Energy Recipe

-[ Related documentation ]-

* COSMO solvation in ADF


Reduction and oxidation potentials
==================================

This example runs the three PLAMS redox recipes for benzoquinone in
water and compares their computed oxidation/reduction potentials.


Import PLAMS recipe classes and utilities.
------------------------------------------

   from typing import Iterable

   from scm.plams import CRSJob, Settings, from_smiles, view
   from scm.plams.recipes.redox import (
       AMSRedoxDirectJob,
       AMSRedoxScreeningJob,
       AMSRedoxThermodynamicCycleJob,
   )


Build benzoquinone and define solvent settings.
-----------------------------------------------

   molecule = from_smiles("C1=CC(=O)C=CC1=O", forcefield="uff")
   view(molecule, width=260, height=220)
   solvent_name = "Water"
   solvent_coskf = CRSJob.database() + "/Water.coskf"
   reduction = True
   oxidation = True
   vibrations = True


Define ADF settings used by the direct and thermodynamic-cycle workflows.
-------------------------------------------------------------------------

   adf_settings = Settings()
   adf_settings.input.adf.basis.type = "DZP"
   adf_settings.input.adf.xc.gga = "PBE"
   adf_settings.input.ams.GeometryOptimization.Convergence.Quality = "Basic"


Create helper printing function and instantiate jobs.
-----------------------------------------------------

   def print_results(jobs: Iterable) -> None:
       she = 4.42  # eV, absolute SHE scale
       print("The experimental reduction potential of benzoquinone is +0.10 V vs. SHE")
       print(
           f"{'Jobname':24s} {'Eox(vib,rel-to-SHE)[V]':24s} {'Ered(vib,rel-to-SHE)[V]':24s} "
           f"{'Eox(rel-to-SHE)[V]':24s} {'Ered(rel-to-SHE)[V]':24s}"
       )

       for job in jobs:
           row = f"{job.name:24s}"
           for vib in [True, False]:
               try:
                   eox = f"{job.results.get_oxidation_potential(vibrations=vib) - she:.2f}"
               except Exception:
                   eox = "N/A"
               row += f" {eox:24s}"

               try:
                   ered = f"{job.results.get_reduction_potential(vibrations=vib) - she:.2f}"
               except Exception:
                   ered = "N/A"
               row += f" {ered:24s}"
           print(row)


   jobs = [
       AMSRedoxScreeningJob(
           name="quick_screening",
           molecule=molecule,
           reduction=reduction,
           oxidation=oxidation,
           solvent_coskf=solvent_coskf,
       ),
       AMSRedoxDirectJob(
           name="direct_best_method",
           molecule=molecule,
           settings=adf_settings,
           reduction=reduction,
           oxidation=oxidation,
           vibrations=vibrations,
           solvent=solvent_name,
       ),
       AMSRedoxThermodynamicCycleJob(
           name="thermodynamic_cycle",
           molecule=molecule,
           settings=adf_settings,
           reduction=reduction,
           oxidation=oxidation,
           vibrations=vibrations,
           solvent=solvent_name,
       ),
   ]


Run all workflows and print per-workflow plus summary tables
------------------------------------------------------------

   for job in jobs:
       job.run()

   for job in jobs:
       job.ok()

   [02.04|12:29:07] JOB quick_screening STARTED
   [02.04|12:29:07] JOB quick_screening RUNNING
   [02.04|12:29:07] JOB quick_screening/dftb_go_0 STARTED
   [02.04|12:29:07] JOB quick_screening/dftb_go_0 RUNNING
   [02.04|12:29:08] JOB quick_screening/dftb_go_0 FINISHED
   [02.04|12:29:08] JOB quick_screening/dftb_go_0 SUCCESSFUL
   [02.04|12:29:08] JOB quick_screening/dftb_go_ox STARTED
   [02.04|12:29:08] JOB quick_screening/dftb_go_ox RUNNING
   [02.04|12:29:09] JOB quick_screening/dftb_go_ox FINISHED
   [02.04|12:29:09] JOB quick_screening/dftb_go_ox SUCCESSFUL
   ... output trimmed ....
   [02.04|12:51:11] JOB thermodynamic_cycle/go_red_solvated STARTED
   [02.04|12:51:11] JOB thermodynamic_cycle/go_red_solvated RUNNING
   [02.04|12:54:17] JOB thermodynamic_cycle/go_red_solvated FINISHED
   [02.04|12:54:17] JOB thermodynamic_cycle/go_red_solvated SUCCESSFUL
   [02.04|12:54:17] JOB thermodynamic_cycle/go_red_solvated_sp_vacuum STARTED
   [02.04|12:54:17] JOB thermodynamic_cycle/go_red_solvated_sp_vacuum RUNNING
   [02.04|12:54:24] JOB thermodynamic_cycle/go_red_solvated_sp_vacuum FINISHED
   [02.04|12:54:24] JOB thermodynamic_cycle/go_red_solvated_sp_vacuum SUCCESSFUL
   [02.04|12:54:24] JOB thermodynamic_cycle FINISHED
   [02.04|12:54:24] JOB thermodynamic_cycle SUCCESSFUL

   print("Final summary:")
   print_results(jobs)

   Final summary:
   The experimental reduction potential of benzoquinone is +0.10 V vs. SHE
   Jobname                  Eox(vib,rel-to-SHE)[V]   Ered(vib,rel-to-SHE)[V]  Eox(rel-to-SHE)[V]       Ered(rel-to-SHE)[V]
   quick_screening          2.95                     0.44                     2.95                     0.44
   direct_best_method       2.61                     -0.08                    2.72                     -0.08
   thermodynamic_cycle      2.63                     -0.07                    2.71                     -0.11


See also
========


Related examples
----------------

* Tuning the Range Separation Parameter in ADF

* Charge Transfer Integrals with ADF

* Reorganization Energy Recipe


Related documentation
---------------------

* COSMO solvation in ADF


Python Script
=============

   #!/usr/bin/env python
   # coding: utf-8

   # ## Reduction and oxidation potentials
   #
   # This example runs the three PLAMS redox recipes for benzoquinone in water and compares their computed oxidation/reduction potentials.
   #

   # ### Import PLAMS recipe classes and utilities.
   #

   from typing import Iterable

   from scm.plams import CRSJob, Settings, from_smiles, view
   from scm.plams.recipes.redox import (
       AMSRedoxDirectJob,
       AMSRedoxScreeningJob,
       AMSRedoxThermodynamicCycleJob,
   )


   # ### Build benzoquinone and define solvent settings.
   #

   molecule = from_smiles("C1=CC(=O)C=CC1=O", forcefield="uff")
   view(molecule, width=260, height=220)
   solvent_name = "Water"
   solvent_coskf = CRSJob.database() + "/Water.coskf"
   reduction = True
   oxidation = True
   vibrations = True


   # ### Define ADF settings used by the direct and thermodynamic-cycle workflows.
   #

   adf_settings = Settings()
   adf_settings.input.adf.basis.type = "DZP"
   adf_settings.input.adf.xc.gga = "PBE"
   adf_settings.input.ams.GeometryOptimization.Convergence.Quality = "Basic"


   # ### Create helper printing function and instantiate jobs.


   def print_results(jobs: Iterable) -> None:
       she = 4.42  # eV, absolute SHE scale
       print("The experimental reduction potential of benzoquinone is +0.10 V vs. SHE")
       print(
           f"{'Jobname':24s} {'Eox(vib,rel-to-SHE)[V]':24s} {'Ered(vib,rel-to-SHE)[V]':24s} "
           f"{'Eox(rel-to-SHE)[V]':24s} {'Ered(rel-to-SHE)[V]':24s}"
       )

       for job in jobs:
           row = f"{job.name:24s}"
           for vib in [True, False]:
               try:
                   eox = f"{job.results.get_oxidation_potential(vibrations=vib) - she:.2f}"
               except Exception:
                   eox = "N/A"
               row += f" {eox:24s}"

               try:
                   ered = f"{job.results.get_reduction_potential(vibrations=vib) - she:.2f}"
               except Exception:
                   ered = "N/A"
               row += f" {ered:24s}"
           print(row)


   jobs = [
       AMSRedoxScreeningJob(
           name="quick_screening",
           molecule=molecule,
           reduction=reduction,
           oxidation=oxidation,
           solvent_coskf=solvent_coskf,
       ),
       AMSRedoxDirectJob(
           name="direct_best_method",
           molecule=molecule,
           settings=adf_settings,
           reduction=reduction,
           oxidation=oxidation,
           vibrations=vibrations,
           solvent=solvent_name,
       ),
       AMSRedoxThermodynamicCycleJob(
           name="thermodynamic_cycle",
           molecule=molecule,
           settings=adf_settings,
           reduction=reduction,
           oxidation=oxidation,
           vibrations=vibrations,
           solvent=solvent_name,
       ),
   ]


   # ### Run all workflows and print per-workflow plus summary tables

   for job in jobs:
       job.run()

   for job in jobs:
       job.ok()


   print("Final summary:")
   print_results(jobs)
