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):
    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 (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 plams
# Compute the neutral-anion reorganization energy of pyrrole
# using ADF as computational engine

molecule = Molecule("")

# 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,

# 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("== Results ==")
print(f"Molecule: {}")
print("State A: neutral")
print("State B: anion")
print(f"Reorganization energy: {reorganization_energy:.6f} [{energy_unit}]")
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}')


To execute this PLAMS script:


[22:39:48] PLAMS working folder: /home/robert/workspace/jobs/SCMSUITE-5940/regtest/test.plams/rundir.plams.ReorganizationEnergy/plams_workdir
[22:39:48] JOB pyrrole STARTED
[22:39:48] JOB pyrrole RUNNING
[22:39:48] JOB pyrrole/go_A STARTED
[22:39:48] JOB pyrrole/go_A RUNNING
[22:39:54] JOB pyrrole/go_A FINISHED
[22:39:54] JOB pyrrole/go_A SUCCESSFUL
[22:39:54] JOB pyrrole/go_B STARTED
[22:39:54] JOB pyrrole/go_B RUNNING
[22:40:04] JOB pyrrole/go_B FINISHED
[22:40:04] JOB pyrrole/go_B SUCCESSFUL
[22:40:04] JOB pyrrole/sp_A_geo_B STARTED
[22:40:04] JOB pyrrole/sp_A_geo_B RUNNING
[22:40:06] JOB pyrrole/sp_A_geo_B FINISHED
[22:40:06] JOB pyrrole/sp_A_geo_B SUCCESSFUL
[22:40:06] JOB pyrrole/sp_B_geo_A STARTED
[22:40:06] JOB pyrrole/sp_B_geo_A RUNNING
[22:40:08] JOB pyrrole/sp_B_geo_A FINISHED
[22:40:08] JOB pyrrole/sp_B_geo_A SUCCESSFUL
[22:40:08] JOB pyrrole FINISHED
[22:40:09] JOB pyrrole SUCCESSFUL

== 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

[22:40:09] PLAMS run finished. Goodbye
Test duration in seconds: 21