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:

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