Reorganization Energy Recipe

ADF Reaction 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).

Compute Marcus reorganization energies using the dedicated PLAMS recipe.

Downloads: Notebook | Script ?

Related documentation
../_images/reorganization_energy_3_0_8c28ad6c.png

Initialization

from scm.plams import Molecule, Settings, ReorganizationEnergyJob, view, AMSJob

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

view(molecule, guess_bonds=True, width=200, height=200, direction="along_y")
image generated from notebook

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

[23.03|17:38:48] JOB pyrrole STARTED
[23.03|17:38:48] JOB pyrrole RUNNING
[23.03|17:38:48] JOB pyrrole/go_A STARTED
[23.03|17:38:48] JOB pyrrole/go_A RUNNING
[23.03|17:38:53] JOB pyrrole/go_A FINISHED
[23.03|17:38:53] JOB pyrrole/go_A SUCCESSFUL
[23.03|17:38:53] JOB pyrrole/go_B STARTED
[23.03|17:38:53] JOB pyrrole/go_B RUNNING
[23.03|17:38:58] JOB pyrrole/go_B FINISHED
[23.03|17:38:58] JOB pyrrole/go_B SUCCESSFUL
[23.03|17:38:58] JOB pyrrole/sp_A_geo_B STARTED
[23.03|17:38:58] JOB pyrrole/sp_A_geo_B RUNNING
[23.03|17:39:01] JOB pyrrole/sp_A_geo_B FINISHED
[23.03|17:39:01] JOB pyrrole/sp_A_geo_B SUCCESSFUL
[23.03|17:39:01] JOB pyrrole/sp_B_geo_A STARTED
[23.03|17:39:01] JOB pyrrole/sp_B_geo_A RUNNING
[23.03|17:39:04] JOB pyrrole/sp_B_geo_A FINISHED
[23.03|17:39:04] JOB pyrrole/sp_B_geo_A SUCCESSFUL
[23.03|17:39:04] JOB pyrrole FINISHED
[23.03|17:39:04] JOB pyrrole SUCCESSFUL

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

See also

Python Script

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

# ## Initialization

from scm.plams import Molecule, Settings, ReorganizationEnergyJob, view, AMSJob


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

view(molecule, guess_bonds=True, width=200, height=200, direction="along_y", picture_path="picture1.png")


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