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 calculations: 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"])

To follow along, either

Worked Example

Initialization

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

# this line is not required in AMS2025+
init()

PLAMS working folder: /path/plams/ReorganizationEnergy/plams_workdir.002

Define molecule

Normally you would read it from an xyz file, but here is for convenience explicit code

# molecule = Molecule("pyrrole.xyz")


def get_molecule(input_string):
    job = AMSJob.from_input(input_string)
    return job.molecule[""]


molecule = get_molecule(
    """
System
    Atoms
        C      -1.12843000       0.00000000      -0.35463200
        C      -0.71293000       0.00000000       0.96463800
        C       0.71293000       0.00000000       0.96463800
        C       1.12843000       0.00000000      -0.35463200
        N       0.00000000       0.00000000      -1.14563200
        H       0.00000000       0.00000000      -2.15713200
        H      -2.12074000       0.00000000      -0.79100200
        H      -1.36515000       0.00000000       1.83237800
        H       1.36515000       0.00000000       1.83237800
        H       2.12074000       0.00000000      -0.79100200
    End
End
"""
)

molecule.properties.name = "pyrrole"  # normally the name of the xyz file

Setup and run job

# 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()

[05.03|13:59:19] JOB pyrrole STARTED
[05.03|13:59:19] JOB pyrrole RUNNING
[05.03|13:59:19] JOB pyrrole/go_A STARTED
[05.03|13:59:19] JOB pyrrole/go_A RUNNING
[05.03|13:59:23] JOB pyrrole/go_A FINISHED
[05.03|13:59:23] JOB pyrrole/go_A SUCCESSFUL
[05.03|13:59:23] JOB pyrrole/go_B STARTED
[05.03|13:59:23] JOB pyrrole/go_B RUNNING
[05.03|13:59:30] JOB pyrrole/go_B FINISHED
[05.03|13:59:30] JOB pyrrole/go_B SUCCESSFUL
... (PLAMS log lines truncated) ...




<scm.plams.recipes.reorganization_energy.ReorganizationEnergyResults at 0x7e91b3f3b0d0>

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("")

== Results ==

Molecule: pyrrole
State A: neutral
State B: anion

Reorganization energy: 0.473683 [eV]

|   State   | Optim Geo | Energy [eV]
|     A     |     A     | -63.801633
|     A     |     B     | -63.487503
|     B     |     A     | -61.702138
|     B     |     B     | -61.861691

Complete Python code

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

# ## Initialization

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

# this line is not required in AMS2025+
init()


# ## Define molecule
#
# Normally you would read it from an xyz file, but here is for convenience explicit code

# molecule = Molecule("pyrrole.xyz")


def get_molecule(input_string):
    job = AMSJob.from_input(input_string)
    return job.molecule[""]


molecule = get_molecule(
    """
System
    Atoms
        C      -1.12843000       0.00000000      -0.35463200
        C      -0.71293000       0.00000000       0.96463800
        C       0.71293000       0.00000000       0.96463800
        C       1.12843000       0.00000000      -0.35463200
        N       0.00000000       0.00000000      -1.14563200
        H       0.00000000       0.00000000      -2.15713200
        H      -2.12074000       0.00000000      -0.79100200
        H      -1.36515000       0.00000000       1.83237800
        H       1.36515000       0.00000000       1.83237800
        H       2.12074000       0.00000000      -0.79100200
    End
End
"""
)

molecule.properties.name = "pyrrole"  # normally the name of the xyz file


# ## Setup and run job

# 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("")