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 childern jobs.
from collections import OrderedDict
from ..core.functions import add_to_instance
from ..core.basejob import MultiJob
from ..core.results import Results
from ..core.settings import Settings
from ..mol.molecule import Molecule
from ..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):
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:
# 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:
Download ReorganizationEnergy.py
Download pyrrole.xyz
$AMSBIN/plams ReorganizationEnergy.py
Output
[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