Reorganization Energy

One of the ingredients for computing hopping rates in Marcus theory is the reorganization energy \(\lambda\), defined as

\[ \begin{align}\begin{aligned}\lambda = & (E^\text{state B}_\text{at optimal geometry of state A} - E^\text{state A}_\text{at optimal geometry of state A}) +\\& (E^\text{state A}_\text{at optimal geometry of state B} - E^\text{state B}_\text{at optimal geometry of state B})\end{aligned}\end{align} \]

where states A and B are two electronic configurations, e.g. neutral and anion (see the example usage below).

In this recipe we build a job class ReorganizationEnergyJob by extending MultiJob. Our job will perform four AMSJob calcualtions: two geometry optimizations for states A anb B, followed by two single point calculations (state A at the optimal geometry of state B and state B at the optimal geometry of state A).

In ReorganizationEnergyResults, the reorganization energy is computed by fetching and combining the results from the children jobs.

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.interfaces.adfsuite.ams import AMSJob

__all__ = ["ReorganizationEnergyJob", "ReorganizationEnergyResults"]


# using this function to pass a molecule inside a MultiJob ensures proper parallel execution
def pass_molecule(source, target):
    @add_to_instance(target)
    def prerun(self):  # noqa F811
        self.molecule = source.results.get_main_molecule()


class ReorganizationEnergyResults(Results):
    """Results class for reorganization energy."""

    def reorganization_energy(self, unit="au"):
        energies = self.get_all_energies(unit)

        reorganization_energy = (
            energies["state B geo A"]
            - energies["state A geo A"]
            + energies["state A geo B"]
            - energies["state B geo B"]
        )
        return reorganization_energy

    def get_all_energies(self, unit="au"):
        energies = {}
        energies["state A geo A"] = self.job.children["go_A"].results.get_energy(unit=unit)
        energies["state B geo B"] = self.job.children["go_B"].results.get_energy(unit=unit)
        energies["state A geo B"] = self.job.children["sp_A_for_geo_B"].results.get_energy(unit=unit)
        energies["state B geo A"] = self.job.children["sp_B_for_geo_A"].results.get_energy(unit=unit)

        return energies


class ReorganizationEnergyJob(MultiJob):
    """A class for calculating the reorganization energy using AMS.
    Given two states, A and B, the reorganization energy is defined as follows:

    reorganization energy =
    E(state B at optimal geometry for state A) -
    E(state A at optimal geometry for state A) +
    E(state A at optimal geometry for state B) -
    E(state B at optimal geometry for state B)

    This job will run two geometry optimizations and two single point calculations.
    """

    _result_type = ReorganizationEnergyResults

    def __init__(self, molecule, common_settings, settings_state_A, settings_state_B, **kwargs):
        """
        molecule: the molecule
        common_settings: a setting object for an AMSJob containing the shared settings for all the calculations
        settings_state_A: Setting object for an AMSJob containing exclusivelt the options defining the state A (e.g. charge and spin)
        settings_state_B: Setting object for an AMSJob containing exclusivelt the options defining the state B (e.g. charge and spin)
        kwargs: other options to be passed to the MultiJob constructor
        """
        MultiJob.__init__(self, children=OrderedDict(), **kwargs)

        go_settings = common_settings.copy()
        go_settings.input.ams.task = "GeometryOptimization"

        sp_settings = common_settings.copy()
        sp_settings.input.ams.task = "SinglePoint"

        # copy the settings so that we wont modify the ones provided as input by the user
        settings_state_A = settings_state_A.copy()
        settings_state_B = settings_state_B.copy()
        # In case the charge key is not specified, excplicitely set the value to 0.
        # This is to prevent the charge in molecule.properties.charge (set by get_main_molecule())
        # to be used in case of neutral systems
        for s in [settings_state_A, settings_state_B]:
            if not "charge" in s.input.ams.system:
                s.input.ams.system.charge = 0

        self.children["go_A"] = AMSJob(molecule=molecule, settings=go_settings + settings_state_A, name="go_A")
        self.children["go_B"] = AMSJob(molecule=molecule, settings=go_settings + settings_state_B, name="go_B")
        self.children["sp_A_for_geo_B"] = AMSJob(settings=sp_settings + settings_state_A, name="sp_A_geo_B")
        self.children["sp_B_for_geo_A"] = AMSJob(settings=sp_settings + settings_state_B, name="sp_B_geo_A")
        pass_molecule(self.children["go_A"], self.children["sp_B_for_geo_A"])
        pass_molecule(self.children["go_B"], self.children["sp_A_for_geo_B"])

Example usage:

#!/usr/bin/env amspython
from scm.plams import Molecule, Settings, ReorganizationEnergyJob

# Compute the neutral-anion reorganization energy of pyrrole
# using ADF as computational engine

molecule = Molecule("pyrrole.xyz")

# Generic settings of the calculation
# (for quantitatively better results, use better settings)
common_settings = Settings()
common_settings.input.adf.Basis.Type = "DZ"

# Specific settings for the neutral calculation.
# Nothing special needs to be done for the neutral calculation,
# so we just use an empty settings.
neutral_settings = Settings()

# Specific settings for the anion calculation:
anion_settings = Settings()
anion_settings.input.ams.System.Charge = -1
anion_settings.input.adf.Unrestricted = "Yes"
anion_settings.input.adf.SpinPolarization = 1

# Create and run the ReorganizationEnergyJob:
job = ReorganizationEnergyJob(
    molecule, common_settings, neutral_settings, anion_settings, name=molecule.properties.name
)
job.run()

# Fetch and print the results:

energy_unit = "eV"
energies = job.results.get_all_energies(energy_unit)
reorganization_energy = job.results.reorganization_energy(energy_unit)

print("")
print("== Results ==")
print("")
print(f"Molecule: {molecule.properties.name}")
print("State A: neutral")
print("State B: anion")
print("")
print(f"Reorganization energy: {reorganization_energy:.6f} [{energy_unit}]")
print("")
print(f"|   State   | Optim Geo | Energy [{energy_unit}]")
print(f'|     A     |     A     | {energies["state A geo A"]:.6f}')
print(f'|     A     |     B     | {energies["state A geo B"]:.6f}')
print(f'|     B     |     A     | {energies["state B geo A"]:.6f}')
print(f'|     B     |     B     | {energies["state B geo B"]:.6f}')
print("")

Note

To execute this PLAMS script: