Reduction and oxidation potentials¶

Note: This example requires AMS2023 or later.

Definitions and introduction¶

Reduction potential: A + e^– → A^– ; ΔG⁰ = –nFE⁰

Oxidation potential: A → A⁺ + e^– ; ΔG⁰ = nFE⁰

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.

Technical

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

Redox potential results¶

Results: The below script outputs the following table:

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.62                     -0.08                    2.72                     -0.09
thermodynamic_cycle      2.64                     -0.07                    2.71                     -0.10

The first two columns give the oxidation and reduction potentials incorporating vibrational effects in the calculation. The last two columns ignore the vibrational effects. Here, we are only interested in the two Ered columns (reduction potentials). We see:

Compared to the experimental E⁰ = +0.10 V, the screening method gives a higher value (+0.44 V) and the more accurate methods a lower value (-0.08 V or -0.07 V).
For this molecule, vibrational effects on the reduction potential are very small (e.g. the vibrations cause E⁰ to go from -0.09 V to -0.08 V for the direct method). Because calculating the vibrations is computationally expensive, they could have been turned off by setting vibrations = False in the below script.
The oxidation potentials (Eox) are given for demonstration purposes only. Turn them off by setting oxidation = False in the below script.

Code example¶

Download ReductionOxidationPotentials.py

#!/usr/bin/env amspython
from scm.plams import from_smiles, CRSJob, Settings
from scm.plams.recipes.redox import AMSRedoxScreeningJob, AMSRedoxDirectJob, AMSRedoxThermodynamicCycleJob


def main():
    mol = from_smiles("C1=CC(=O)C=CC1=O", forcefield="uff")  # benzoquinone
    solvent_name = "Water"  # Solvent for AMSRedoxDirectJob and AMSRedoxThermodynamicCycleJob. See the ADF Solvation documentation for which solvents are available.
    solvent_coskf = CRSJob.database() + "/Water.coskf"  # .coskf file or AMSRedoxScreeningJob

    vibrations = True  # set to False to turn off vibrations for AMSRedoxDirectJob and AMSRedoxThermodynamicCycleJob
    reduction = True
    oxidation = True  # not really relevant for benzoquinone, set to False to not calculate

    # DFT settings for AMSRedoxDirectJob and AMSRedoxThermodynamicCycleJob
    # Note: These are for demonstration purposes only. Use a better functional/settings for production purposes.
    s = Settings()
    s.input.adf.basis.type = "DZP"
    s.input.adf.xc.gga = "PBE"
    s.input.ams.GeometryOptimization.Convergence.Quality = "Basic"

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

    for job in jobs:
        job.run()
        print_results([job])

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


def print_results(jobs):
    SHE = 4.42  # standard hydrogen electrode in eV on absolute scale

    print("The experimental reduction potential of benzoquinone is +0.10 V vs. SHE")
    print(
        "{:24s} {:24s} {:24s} {:24s} {:24s}".format(
            "Jobname", "Eox(vib,rel-to-SHE)[V]", "Ered(vib,rel-to-SHE)[V]", "Eox(rel-to-SHE)[V]", "Ered(rel-to-SHE)[V]"
        )
    )

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

            try:
                Ered = job.results.get_reduction_potential(vibrations=vibrations) - SHE
                Ered = f"{Ered:.2f}"
            except:
                Ered = "N/A"
            s += f" {Ered:24s}"

        print(s)


if __name__ == "__main__":
    main()

The PLAMS recipes/classes¶

This section contains some details about the implementation.

The oxidation and reduction potentials are calculated as follows:

class AMSRedoxDirectResults(AMSRedoxParentResults):
    def get_oxidation_potential(self, vibrations=True, unit="eV"):
        Greact = self._get_energy(self.job.children["job_0"], vibrations=vibrations)
        Gprod = self._get_energy(self.job.children["job_ox"], vibrations=vibrations) + self.Gelectron
        ret = (Gprod - Greact) * Units.convert(1.0, "hartree", unit)
        return ret

    def get_reduction_potential(self, vibrations=True, unit="eV"):
        Greact = self._get_energy(self.job.children["job_0"], vibrations=vibrations) + self.Gelectron
        Gprod = self._get_energy(self.job.children["job_red"], vibrations=vibrations)
        ret = (Gprod - Greact) * Units.convert(1.0, "hartree", unit)
        ret *= -1  # deltaG = -nFE
        return ret

class AMSRedoxThermodynamicCycleResults(AMSRedoxParentResults):
    def get_oxidation_potential(self, vibrations=True, unit="eV"):
        Greact = (
            self._get_energy(self.job.children["go_0_vacuum"], vibrations=vibrations)
            + self._get_energy(self.job.children["go_0_vacuum_sp_solvated"], vibrations=False)
            + self._get_energy(self.job.children["go_0_solvated_sp_vacuum"], vibrations=False)
            - 2 * self._get_energy(self.job.children["go_0_vacuum"], vibrations=False)
        )

        Gprod = (
            self._get_energy(self.job.children["go_ox_vacuum"], vibrations=vibrations)
            + self._get_energy(self.job.children["go_ox_vacuum_sp_solvated"], vibrations=False)
            + self._get_energy(self.job.children["go_ox_solvated_sp_vacuum"], vibrations=False)
            - 2 * self._get_energy(self.job.children["go_ox_vacuum"], vibrations=False)
            + self.Gelectron
        )

        ret = (Gprod - Greact) * Units.convert(1.0, "hartree", unit)
        return ret

    def get_reduction_potential(self, vibrations=True, unit="eV"):
        Greact = (
            self._get_energy(self.job.children["go_0_vacuum"], vibrations=vibrations)
            + self._get_energy(self.job.children["go_0_vacuum_sp_solvated"], vibrations=False)
            + self._get_energy(self.job.children["go_0_solvated_sp_vacuum"], vibrations=False)
            - 2 * self._get_energy(self.job.children["go_0_vacuum"], vibrations=False)
            + self.Gelectron
        )

        Gprod = (
            self._get_energy(self.job.children["go_red_vacuum"], vibrations=vibrations)
            + self._get_energy(self.job.children["go_red_vacuum_sp_solvated"], vibrations=False)
            + self._get_energy(self.job.children["go_red_solvated_sp_vacuum"], vibrations=False)
            - 2 * self._get_energy(self.job.children["go_red_vacuum"], vibrations=False)
        )

        ret = (Gprod - Greact) * Units.convert(1.0, "hartree", unit)
        ret *= -1  # deltaG = -nFE
        return ret

where

go_0_vacuum is a geometry optimization for the neutral (0) molecule in vacuum,
go_0_vaccum_sp_solvated is a single point in implicit solvent on the previous structure,
etc.

class AMSRedoxScreeningResults(Results):
    Gelectron = -0.0375 / 27.211  # in hartree

    def get_oxidation_potential(self, vibrations=True, unit="eV"):
        """Note: vibrations cannot be disabled for AMSRedoxScreening"""
        Greact = self.job.children["activitycoef_0"].results.get_gibbs_energy()
        Gprod = self.job.children["activitycoef_ox"].results.get_gibbs_energy() + self.Gelectron
        ret = (Gprod - Greact) * Units.convert(1.0, "hartree", unit)
        return ret

    def get_reduction_potential(self, vibrations=True, unit="eV"):
        """Note: vibrations cannot be disabled for AMSRedoxScreening"""
        Greact = self.job.children["activitycoef_0"].results.get_gibbs_energy() + self.Gelectron
        Gprod = self.job.children["activitycoef_red"].results.get_gibbs_energy()
        ret = (Gprod - Greact) * Units.convert(1.0, "hartree", unit)
        ret *= -1  # deltaG = -nFE
        return ret

The complete class definitions:

import os
import shutil
from collections import OrderedDict

from scm.plams.core.basejob import MultiJob
from scm.plams.core.functions import add_to_instance
from scm.plams.core.results import Results
from scm.plams.core.settings import Settings
from scm.plams.interfaces.adfsuite.ams import AMSJob
from scm.plams.interfaces.adfsuite.crs import CRSJob, CRSResults
from scm.plams.mol.molecule import Molecule
from scm.plams.recipes.adfcosmorscompound import ADFCOSMORSCompoundJob
from scm.plams.tools.units import Units

"""
    Classes for calculating oxidation potentials and reduction potentials in implicit solvent.

    To get the results, call

    job.results.get_oxidation_potential(vibrations=False)

    or

    job.results.get_reduction_potential(vibrations=False)

    Set ``vibrations=True`` to include vibrational effects if the job was run
    with ``vibrations=True``.

    The potential is returned in V. It is on an absolute scale. To get
    potential relative to SHE, subtract 4.47 V.

    AMSRedoxDirectJob
    --------------------

    Calculates reaction energies for 
    A + e^- --> A^-     if ``reduction=True``,
    and/or
    A --> A^+ + e^-     if ``oxidation=True``.

    using direct geometry optimizations in implicit solvent. The solvent must be one supported by ADF.

    If ``vibrations`` = True, then calculate and use Gibbs free energies at room temperature.

    ``settings`` should be set to the engine settings, excluding spin polarization and implicit solvation.

    Requires an ADF license.

    AMSRedoxThermodynamicCycleJob
    -----------------------------

    Sometimes more efficient (and less accurate) alternative to AMSRedoxDirectJob. Can be
    useful if the molecule is large and if ``vibrations=True``. For small molecules just use AMSRedoxDirectJob.

    Requires an ADF license.

    AMSRedoxScreeningJob
    -----------------------

    The fastest and least accurate option. Geometries are optimized at the DFTB
    level, ADF is then used to generate a .coskf file and the solvation free
    energies are evaluated with COSMO-RS.

    Note: you must supply ``solvent_coskf``, a path to the solvent .coskf file.
    If you do not have one, you can generate one with the
    ``ADFCOSMORSCompoundJob`` recipe.

    Requires DFTB, ADF, and COSMO-RS licenses.

    Note: The functional used is different from Belic et al. This class uses the
    settings that are meant to be used for generating .coskf files.

"""

__all__ = [
    "AMSRedoxScreeningJob",
    "AMSRedoxScreeningResults",
    "AMSRedoxDirectJob",
    "AMSRedoxDirectResults",
    "AMSRedoxThermodynamicCycleJob",
    "AMSRedoxThermodynamicCycleResults",
]


def _spinpol_settings(molecule, charge=0):
    sett = Settings()
    num_electrons = sum(atom.atnum for atom in molecule) - charge
    spinpol = num_electrons % 2
    if spinpol == 0:
        unrestricted = "No"
    else:
        unrestricted = "Yes"
        # sett.input.adf.OCCUPATIONS = "ElectronicTemperature=100"
    sett.input.adf.SpinPolarization = spinpol
    sett.input.adf.Unrestricted = unrestricted
    return sett


class AMSRedoxParentJob(MultiJob):
    def __init__(
        self,
        molecule: Molecule,  # Molecule
        name: str = None,
        settings: Settings = None,
        oxidation: bool = True,
        reduction: bool = False,
    ):

        MultiJob.__init__(self, children=OrderedDict(), name=name)

        self.oxidation = oxidation
        self.reduction = reduction
        self.input_molecule = molecule
        self.settings = settings or Settings()
        self.orig_charge = molecule.properties.get("charge", 0)
        self.ox_charge = self.orig_charge + 1
        self.red_charge = self.orig_charge - 1

    @staticmethod
    def solvation_settings(solvent):
        sett = Settings()
        if solvent is None:
            solvent = "vacuum"
        if isinstance(solvent, str):
            sett.input.adf.Solvation.Solv = f"name={solvent}"
        elif isinstance(solvent, tuple):
            sett.input.adf.Solvation.Solv = f"eps={solvent[0]} rad={solvent[1]}"
        else:
            raise TypeError(f"solvent is of type {type(solvent)}, expected str or 2-tuple")

        return sett

    def get_dft_0_settings(self, vibrations=False, perturbcoordinates=False):
        """all settings except task and solvation"""
        s = Settings()
        s += self.settings
        s.update(_spinpol_settings(self.input_molecule, self.orig_charge))
        s.input.ams.System.Charge = self.orig_charge
        s.input.ams.System.PerturbCoordinates = 0.01 if perturbcoordinates else 0.00
        s.input.ams.Properties.NormalModes = str(vibrations)
        return s

    def get_dft_ox_settings(self, vibrations=False, perturbcoordinates=False):
        sox = self.get_dft_0_settings(vibrations=vibrations, perturbcoordinates=perturbcoordinates)
        sox.update(_spinpol_settings(self.input_molecule, self.ox_charge))
        sox.input.ams.System.Charge = self.ox_charge
        return sox

    def get_dft_red_settings(self, vibrations=False, perturbcoordinates=False):
        sred = self.get_dft_0_settings(vibrations=vibrations, perturbcoordinates=perturbcoordinates)
        sred.update(_spinpol_settings(self.input_molecule, self.red_charge))
        sred.input.ams.System.Charge = self.red_charge
        return sred


class AMSRedoxParentResults(Results):
    Gelectron = -0.0375 / 27.211  # in hartree

    @staticmethod
    def _get_energy(job, vibrations: bool):
        if not vibrations:
            return job.results.get_energy()
        e = job.results.readrkf("Thermodynamics", "Gibbs free Energy", file="engine")
        if isinstance(e, list):
            return e[0]
        return e

    @staticmethod
    def _get_solvation_energy(job):
        dG_solvation = job.results.readrkf("Energy", "Solvation Energy (el)", "adf") + job.results.readrkf(
            "Energy", "Solvation Energy (cd)", "adf"
        )
        return dG_solvation


class CRSActivityCoefficientResults(CRSResults):
    def get_gibbs_energy(self):
        """Gibbs energy in hartree"""
        return float(self.get_results()["G solute"][1]) * 0.00159376  # kcal/mol to hartree


class CRSActivityCoefficientJob(CRSJob):
    _result_type = CRSActivityCoefficientResults

    def __init__(self, solvent_coskf, solute_coskf, name=None, temperature=298.15, copy_coskf=False):
        CRSJob.__init__(self, name=name)
        self.solvent_coskf = solvent_coskf
        self.solute_coskf = solute_coskf
        self.temperature = temperature
        self.copy_coskf = copy_coskf

    def prerun(self):  # noqa F811
        self._prerun()

    def _prerun(self):
        if self.copy_coskf:
            new_solvent = f"solvent_{os.path.basename(self.solvent_coskf)}"
            new_solute = f"solute_{os.path.basename(self.solute_coskf)}"
            shutil.copy(self.solvent_coskf, os.path.join(self.path, new_solvent))
            shutil.copy(self.solute_coskf, os.path.join(self.path, new_solute))
            self.solvent_coskf = new_solvent
            self.solute_coskf = new_solute
        self.settings.input.property._h = "ACTIVITYCOEF"
        compounds = [Settings(), Settings()]
        compounds[0]._h = self.solvent_coskf
        compounds[1]._h = self.solute_coskf
        compounds[0].frac1 = 1
        compounds[1].frac1 = 0

        self.settings.input.temperature = str(self.temperature)
        self.settings.input.compound = compounds


class AMSRedoxDirectResults(AMSRedoxParentResults):
    def get_oxidation_potential(self, vibrations=True, unit="eV"):
        Greact = self._get_energy(self.job.children["job_0"], vibrations=vibrations)
        Gprod = self._get_energy(self.job.children["job_ox"], vibrations=vibrations) + self.Gelectron
        ret = (Gprod - Greact) * Units.convert(1.0, "hartree", unit)
        return ret

    def get_reduction_potential(self, vibrations=True, unit="eV"):
        Greact = self._get_energy(self.job.children["job_0"], vibrations=vibrations) + self.Gelectron
        Gprod = self._get_energy(self.job.children["job_red"], vibrations=vibrations)
        ret = (Gprod - Greact) * Units.convert(1.0, "hartree", unit)
        ret *= -1  # deltaG = -nFE
        return ret


class AMSRedoxDirectJob(AMSRedoxParentJob):
    _result_type = AMSRedoxDirectResults

    def __init__(
        self,
        molecule: Molecule,  # Molecule
        solvent: str,  # ADF pre-defined solvent
        name: str = None,
        settings: Settings = None,
        oxidation: bool = True,
        reduction: bool = False,
        vibrations: bool = True,
    ):

        AMSRedoxParentJob.__init__(
            self, name=name, settings=settings, molecule=molecule, oxidation=oxidation, reduction=reduction
        )

        self.solvent = solvent
        self.vibrations = vibrations

        s = Settings()
        s += self.settings
        s.input.ams.Task = "GeometryOptimization"
        s += self.solvation_settings(self.solvent)
        s.update(_spinpol_settings(self.input_molecule, self.orig_charge))
        s.input.ams.System.Charge = self.orig_charge
        s.input.ams.System.PerturbCoordinates = 0.01
        s.input.ams.Properties.NormalModes = str(self.vibrations)
        job_0 = AMSJob(settings=s, name="job_0", molecule=self.input_molecule)
        self.children["job_0"] = job_0

        if oxidation:
            sox = s.copy()
            sox.update(_spinpol_settings(self.input_molecule, self.ox_charge))
            sox.input.ams.System.Charge = self.ox_charge
            sox.input.ams.System.PerturbCoordinates = 0.01
            job_ox = AMSJob(settings=sox, name="job_ox")

            @add_to_instance(job_ox)
            def prerun(self):  # noqa F811
                self.molecule = job_0.results.get_main_molecule()

            self.children["job_ox"] = job_ox

        if reduction:
            sred = s.copy()
            sred.update(_spinpol_settings(self.input_molecule, self.red_charge))
            sred.input.ams.System.Charge = self.red_charge
            sred.input.ams.System.PerturbCoordinates = 0.01
            job_red = AMSJob(settings=sred, name="job_red")

            @add_to_instance(job_red)
            def prerun(self):  # noqa F811
                self.molecule = job_0.results.get_main_molecule()

            self.children["job_red"] = job_red


class AMSRedoxScreeningResults(Results):
    Gelectron = -0.0375 / 27.211  # in hartree

    def get_oxidation_potential(self, vibrations=True, unit="eV"):
        """Note: vibrations cannot be disabled for AMSRedoxScreening"""
        Greact = self.job.children["activitycoef_0"].results.get_gibbs_energy()
        Gprod = self.job.children["activitycoef_ox"].results.get_gibbs_energy() + self.Gelectron
        ret = (Gprod - Greact) * Units.convert(1.0, "hartree", unit)
        return ret

    def get_reduction_potential(self, vibrations=True, unit="eV"):
        """Note: vibrations cannot be disabled for AMSRedoxScreening"""
        Greact = self.job.children["activitycoef_0"].results.get_gibbs_energy() + self.Gelectron
        Gprod = self.job.children["activitycoef_red"].results.get_gibbs_energy()
        ret = (Gprod - Greact) * Units.convert(1.0, "hartree", unit)
        ret *= -1  # deltaG = -nFE
        return ret


class AMSRedoxScreeningJob(AMSRedoxParentJob):
    _result_type = AMSRedoxScreeningResults

    def __init__(
        self,
        molecule: Molecule,  # Molecule
        solvent_coskf: str,  # path to solvent .coskf file
        name: str = None,
        settings: Settings = None,
        oxidation: bool = True,
        reduction: bool = False,
        copy_coskf=False,
    ):
        """
        copy_coskf: bool
            Will copy the coskf files into the job's own directory before running.
        """
        AMSRedoxParentJob.__init__(
            self, name=name, settings=settings, molecule=molecule, oxidation=oxidation, reduction=reduction
        )

        self.solvent_coskf = solvent_coskf

        # dftb_go_0
        s = Settings()
        s.input.ams.Task = "GeometryOptimization"
        s.input.DFTB.Model = "GFN1-xTB"
        s.input.ams.System.Charge = self.orig_charge
        s.input.ams.System.PerturbCoordinates = 0.01
        dftb_go_0 = AMSJob(settings=s, molecule=self.input_molecule, name="dftb_go_0")
        self.children["dftb_go_0"] = dftb_go_0

        if oxidation:
            sox = s.copy()
            sox.input.ams.System.PerturbCoordinates = 0.01
            sox.input.ams.System.Charge = self.ox_charge
            dftb_go_ox = AMSJob(settings=sox, name="dftb_go_ox")

            @add_to_instance(dftb_go_ox)
            def prerun(self):  # noqa F811
                self.molecule = dftb_go_0.results.get_main_molecule()

            self.children["dftb_go_ox"] = dftb_go_ox

        if reduction:
            sred = s.copy()
            sred.input.ams.System.PerturbCoordinates = 0.01
            sred.input.ams.System.Charge = self.red_charge
            dftb_go_red = AMSJob(settings=sred, name="dftb_go_red")

            @add_to_instance(dftb_go_red)
            def prerun(self):  # noqa F811
                self.molecule = dftb_go_0.results.get_main_molecule()

            self.children["dftb_go_red"] = dftb_go_red

        # gencoskf_0
        s = self.settings.copy()
        s.update(_spinpol_settings(self.input_molecule, self.orig_charge))
        # s.input.ams.System.Charge = self.orig_charge
        gencoskf_0 = ADFCOSMORSCompoundJob(molecule=None, name="gencoskf_0", singlepoint=True, settings=s)

        @add_to_instance(gencoskf_0)
        def prerun(self):  # noqa F811
            self.input_molecule = dftb_go_0.results.get_main_molecule()  # will inherit the charge
            assert isinstance(self.input_molecule, Molecule)
            self.atomic_ion = len(self.input_molecule) == 1

        self.children["gencoskf_0"] = gencoskf_0

        if oxidation:
            sox = s.copy()
            sox.update(_spinpol_settings(self.input_molecule, self.ox_charge))
            # sox.input.ams.System.Charge = self.ox_charge
            gencoskf_ox = ADFCOSMORSCompoundJob(molecule=None, name="gencoskf_ox", singlepoint=True, settings=sox)

            @add_to_instance(gencoskf_ox)
            def prerun(self):  # noqa F811
                self.input_molecule = dftb_go_ox.results.get_main_molecule()
                self.atomic_ion = len(self.input_molecule) == 1

            self.children["gencoskf_ox"] = gencoskf_ox

        if reduction:
            sred = s.copy()
            sred.update(_spinpol_settings(self.input_molecule, self.red_charge))
            # sred.input.ams.System.Charge = self.red_charge
            gencoskf_red = ADFCOSMORSCompoundJob(molecule=None, name="gencoskf_red", singlepoint=True, settings=sred)

            @add_to_instance(gencoskf_red)
            def prerun(self):  # noqa F811
                self.input_molecule = dftb_go_red.results.get_main_molecule()
                self.atomic_ion = len(self.input_molecule) == 1

            self.children["gencoskf_red"] = gencoskf_red

        # activitycoef_0
        activitycoef_0 = CRSActivityCoefficientJob(
            name="activitycoef_0", solvent_coskf=self.solvent_coskf, solute_coskf=None, copy_coskf=copy_coskf
        )

        @add_to_instance(activitycoef_0)
        def prerun(self):  # noqa F811
            self.solute_coskf = gencoskf_0.results.coskfpath()
            self._prerun()

        self.children["activitycoef_0"] = activitycoef_0

        if oxidation:
            activitycoef_ox = CRSActivityCoefficientJob(
                name="activitycoef_ox", solvent_coskf=self.solvent_coskf, solute_coskf=None, copy_coskf=copy_coskf
            )

            @add_to_instance(activitycoef_ox)
            def prerun(self):  # noqa F811
                self.solute_coskf = gencoskf_ox.results.coskfpath()
                self._prerun()

            self.children["activitycoef_ox"] = activitycoef_ox

        if reduction:
            activitycoef_red = CRSActivityCoefficientJob(
                name="activitycoef_red", solvent_coskf=self.solvent_coskf, solute_coskf=None, copy_coskf=copy_coskf
            )

            @add_to_instance(activitycoef_red)
            def prerun(self):  # noqa F811
                self.solute_coskf = gencoskf_red.results.coskfpath()
                self._prerun()

            self.children["activitycoef_red"] = activitycoef_red


class AMSRedoxThermodynamicCycleResults(AMSRedoxParentResults):
    def get_oxidation_potential(self, vibrations=True, unit="eV"):
        Greact = (
            self._get_energy(self.job.children["go_0_vacuum"], vibrations=vibrations)
            + self._get_energy(self.job.children["go_0_vacuum_sp_solvated"], vibrations=False)
            + self._get_energy(self.job.children["go_0_solvated_sp_vacuum"], vibrations=False)
            - 2 * self._get_energy(self.job.children["go_0_vacuum"], vibrations=False)
        )

        Gprod = (
            self._get_energy(self.job.children["go_ox_vacuum"], vibrations=vibrations)
            + self._get_energy(self.job.children["go_ox_vacuum_sp_solvated"], vibrations=False)
            + self._get_energy(self.job.children["go_ox_solvated_sp_vacuum"], vibrations=False)
            - 2 * self._get_energy(self.job.children["go_ox_vacuum"], vibrations=False)
            + self.Gelectron
        )

        ret = (Gprod - Greact) * Units.convert(1.0, "hartree", unit)
        return ret

    def get_reduction_potential(self, vibrations=True, unit="eV"):
        Greact = (
            self._get_energy(self.job.children["go_0_vacuum"], vibrations=vibrations)
            + self._get_energy(self.job.children["go_0_vacuum_sp_solvated"], vibrations=False)
            + self._get_energy(self.job.children["go_0_solvated_sp_vacuum"], vibrations=False)
            - 2 * self._get_energy(self.job.children["go_0_vacuum"], vibrations=False)
            + self.Gelectron
        )

        Gprod = (
            self._get_energy(self.job.children["go_red_vacuum"], vibrations=vibrations)
            + self._get_energy(self.job.children["go_red_vacuum_sp_solvated"], vibrations=False)
            + self._get_energy(self.job.children["go_red_solvated_sp_vacuum"], vibrations=False)
            - 2 * self._get_energy(self.job.children["go_red_vacuum"], vibrations=False)
        )

        ret = (Gprod - Greact) * Units.convert(1.0, "hartree", unit)
        ret *= -1  # deltaG = -nFE
        return ret


class AMSRedoxThermodynamicCycleJob(AMSRedoxParentJob):
    _result_type = AMSRedoxThermodynamicCycleResults

    def __init__(
        self,
        molecule: Molecule,  # Molecule
        name: str = None,
        solvent: str = "Water",
        settings: Settings = None,
        oxidation: bool = True,
        reduction: bool = False,
        vibrations: bool = False,
    ):

        AMSRedoxParentJob.__init__(
            self, name=name, settings=settings, molecule=molecule, oxidation=oxidation, reduction=reduction
        )

        self.vibrations = vibrations
        self.solvent = solvent

        # go_0_vacuum
        s = self.get_dft_0_settings(vibrations=self.vibrations, perturbcoordinates=True)
        s.input.ams.Task = "GeometryOptimization"
        go_0_vacuum = AMSJob(settings=s, molecule=self.input_molecule, name="go_0_vacuum")

        @add_to_instance(go_0_vacuum)
        def prerun(self):  # noqa F811
            self.molecule = self.parent.input_molecule

        self.children["go_0_vacuum"] = go_0_vacuum

        # go_0_vacuum_sp_solvated
        ss = self.get_dft_0_settings(vibrations=False, perturbcoordinates=False)
        ss.input.ams.Task = "SinglePoint"
        ss += self.solvation_settings(self.solvent)
        go_0_vacuum_sp_solvated = AMSJob(settings=ss, molecule=None, name="go_0_vacuum_sp_solvated")

        @add_to_instance(go_0_vacuum_sp_solvated)
        def prerun(self):  # noqa F811
            self.molecule = go_0_vacuum.results.get_main_molecule()

        self.children["go_0_vacuum_sp_solvated"] = go_0_vacuum_sp_solvated

        # go_0_solvated
        s = self.get_dft_0_settings(vibrations=self.vibrations, perturbcoordinates=True)
        s.input.ams.Task = "GeometryOptimization"
        s += self.solvation_settings(self.solvent)
        go_0_solvated = AMSJob(settings=s, molecule=None, name="go_0_solvated")

        @add_to_instance(go_0_solvated)
        def prerun(self):  # noqa F811
            self.molecule = go_0_vacuum.results.get_main_molecule()

        self.children["go_0_solvated"] = go_0_solvated

        # go_0_solvated_sp_vacuum
        ss = self.get_dft_0_settings(vibrations=False, perturbcoordinates=False)
        ss.input.ams.Task = "SinglePoint"
        go_0_solvated_sp_vacuum = AMSJob(settings=ss, molecule=None, name="go_0_solvated_sp_vacuum")

        @add_to_instance(go_0_solvated_sp_vacuum)
        def prerun(self):  # noqa F811
            self.molecule = go_0_solvated.results.get_main_molecule()

        self.children["go_0_solvated_sp_vacuum"] = go_0_solvated_sp_vacuum

        if oxidation:
            # go_ox_vacuum
            s = self.get_dft_ox_settings(vibrations=self.vibrations, perturbcoordinates=True)
            s.input.ams.Task = "GeometryOptimization"
            go_ox_vacuum = AMSJob(settings=s, molecule=self.input_molecule, name="go_ox_vacuum")

            @add_to_instance(go_ox_vacuum)
            def prerun(self):  # noqa F811
                self.molecule = go_0_vacuum.results.get_main_molecule()

            self.children["go_ox_vacuum"] = go_ox_vacuum

            # go_ox_vacuum_sp_solvated
            ss = self.get_dft_ox_settings(vibrations=False, perturbcoordinates=False)
            ss.input.ams.Task = "SinglePoint"
            ss += self.solvation_settings(self.solvent)
            go_ox_vacuum_sp_solvated = AMSJob(settings=ss, molecule=None, name="go_ox_vacuum_sp_solvated")

            @add_to_instance(go_ox_vacuum_sp_solvated)
            def prerun(self):  # noqa F811
                self.molecule = go_ox_vacuum.results.get_main_molecule()

            self.children["go_ox_vacuum_sp_solvated"] = go_ox_vacuum_sp_solvated

            # go_ox_solvated
            s = self.get_dft_ox_settings(vibrations=self.vibrations, perturbcoordinates=True)
            s.input.ams.Task = "GeometryOptimization"
            s += self.solvation_settings(self.solvent)
            go_ox_solvated = AMSJob(settings=s, molecule=None, name="go_ox_solvated")

            @add_to_instance(go_ox_solvated)
            def prerun(self):  # noqa F811
                self.molecule = go_ox_vacuum.results.get_main_molecule()

            self.children["go_ox_solvated"] = go_ox_solvated

            # go_ox_solvated_sp_vacuum
            ss = self.get_dft_ox_settings(vibrations=False, perturbcoordinates=False)
            ss.input.ams.Task = "SinglePoint"
            go_ox_solvated_sp_vacuum = AMSJob(settings=ss, molecule=None, name="go_ox_solvated_sp_vacuum")

            @add_to_instance(go_ox_solvated_sp_vacuum)
            def prerun(self):  # noqa F811
                self.molecule = go_ox_solvated.results.get_main_molecule()

            self.children["go_ox_solvated_sp_vacuum"] = go_ox_solvated_sp_vacuum

        if reduction:
            # go_red_vacuum
            s = self.get_dft_red_settings(vibrations=self.vibrations, perturbcoordinates=True)
            s.input.ams.Task = "GeometryOptimization"
            go_red_vacuum = AMSJob(settings=s, molecule=self.input_molecule, name="go_red_vacuum")

            @add_to_instance(go_red_vacuum)
            def prerun(self):  # noqa F811
                self.molecule = go_0_vacuum.results.get_main_molecule()

            self.children["go_red_vacuum"] = go_red_vacuum

            # go_red_vacuum_sp_solvated
            ss = self.get_dft_red_settings(vibrations=False, perturbcoordinates=False)
            ss.input.ams.Task = "SinglePoint"
            ss += self.solvation_settings(self.solvent)
            go_red_vacuum_sp_solvated = AMSJob(settings=ss, molecule=None, name="go_red_vacuum_sp_solvated")

            @add_to_instance(go_red_vacuum_sp_solvated)
            def prerun(self):  # noqa F811
                self.molecule = go_red_vacuum.results.get_main_molecule()

            self.children["go_red_vacuum_sp_solvated"] = go_red_vacuum_sp_solvated

            # go_red_solvated
            s = self.get_dft_red_settings(vibrations=self.vibrations, perturbcoordinates=True)
            s.input.ams.Task = "GeometryOptimization"
            s += self.solvation_settings(self.solvent)
            go_red_solvated = AMSJob(settings=s, molecule=None, name="go_red_solvated")

            @add_to_instance(go_red_solvated)
            def prerun(self):  # noqa F811
                self.molecule = go_red_vacuum.results.get_main_molecule()

            self.children["go_red_solvated"] = go_red_solvated

            # go_red_solvated_sp_vacuum
            ss = self.get_dft_red_settings(vibrations=False, perturbcoordinates=False)
            ss.input.ams.Task = "SinglePoint"
            go_red_solvated_sp_vacuum = AMSJob(settings=ss, molecule=None, name="go_red_solvated_sp_vacuum")

            @add_to_instance(go_red_solvated_sp_vacuum)
            def prerun(self):  # noqa F811
                self.molecule = go_red_solvated.results.get_main_molecule()

            self.children["go_red_solvated_sp_vacuum"] = go_red_solvated_sp_vacuum

Details about individual calculations¶

The paper by Belic et al. describes four ways to calculate the oxidation potential which are described in the tabs below, however we have extended them to allow for calculation of the reduction potential. In general, \(g_p^i\) denotes the optimised geometry of the molecule in phase \(p\) (solvent [sol] or gaseous [gas]) and in state \(i\) (oxidised [+]/reduced [-] or neutral [0]), \(G_P^I(g_p^i)\) denotes the Gibbs free energy of the geometry \(g_p^i\) calculated in phase \(P\) and state \(I\), similarly \(E_P(g_p^i)\) denotes the bond energy of the geometry \(g_p^i\) calculated in phase \(P\). The electron Gibbs free energy \(G_{gas}(e^-)=0.0375 \text{eV}\) is also used in these calculations.

Direct method of calculating oxidation potential using the COSMO solvation model.

\[\Delta G^{DC}_{COSMO} = G_{sol}^\pm(g_{sol}^\pm) \pm G_{gas}(e^-) - G_{sol}^0(g_{sol}^0)\]

The following steps are calculated

Step	Task	Structure	Method	I	P	Frequencies
1	GO		DFT	0	sol	Yes
2	GO		DFT	\(\pm\)	sol	Yes

Thermodynamic cycle using either the COSMO or the COSMO-RS solvation model.

\[\Delta G^{TC}_{COSMO/COSMO-RS} = G_{sol}^\pm(g_{gas}^\pm) \pm G_{gas}(e^-) - G_{sol}^0(g_{gas}^0)\]

where we have

\[G_{sol}^i = G_{gas}^i(g_{gas}^i) + \Delta G_{sol}^i(g_{gas}^i) + [E_{gas}(g_{sol}^i) - E_{gas}(g_{gas}^i)].\]

The following steps are calculated

Step	Task	Structure	Method	I	P	Frequencies
1	GO		DFT	0	gas	Yes
2	SP	\(g^0_{gas}\), step 1	DFT	0	sol	No
3	GO		DFT	0	sol	No
4	SP	\(g^0_{sol}\), step 3	DFT	0	gas	No
5	GO		DFT	\(\pm\)	gas	Yes
6	SP	\(g^\pm_{gas}\), step 5	DFT	\(\pm\)	sol	No
7	GO		DFT	\(\pm\)	sol	No
8	SP	\(g^\pm_{sol}\), step 7	DFT	\(\pm\)	gas	No

This method also implements the thermodynamic cycle, but uses lower theory methods (DFTB) to greatly increase speed of calculation. This method uses COSMO-RS as the solvation method.

\[\Delta G_{COSMO-RS}^{screening} = G_{sol, CRS}^\pm(g_{gas}^\pm) \pm G_{gas}(e^-) - G_{sol, CRS}^0(g_{gas}^0)\]

where we have

\[G_{sol, CRS}^i = E_{sol}^i(g_{gas}^i) + \Delta G_{sol, CRS}^i(g_{gas}^i])\]

Step	Task	Structure	Method	I	P	Frequencies
1	GO		DFTB	0	gas	No
2	SP	\(g^0_{gas}\), step 1	DFT	0	sol	No
3	GO		DFTB	\(\pm\)	gas	No
4	SP	\(g^\pm_{gas}\), step 3	DFT	\(\pm\)	sol	No

Electronic Structure

ADF

Periodic DFT

DFTB & MOPAC

Interatomic Potentials

ReaxFF

Machine Learning Potentials

Force Fields

kMC and Microkinetics

Bumblebee: OLED stacks

Fluid Thermodynamics

COSMO-RS

Workflows and Utilities

OLED workflows

ChemTraYzer2

Conformers

Reactions Discovery

AMS Driver

Properties

PES Exploration

Molecular Dynamics

Monte Carlo

Interfaces

ParAMS

PLAMS

GUI

VASP

Downloads

Windows

Mac

Linux

Documentation

Overview

Tutorials

Installation Manual

Brochures

Other Resources

Changelog

Workshops

Knowledgebank

FAQ

Pricing and licensing