Reorganization Energy Recipe
****************************

ADF Reaction energy

One of the ingredients for computing hopping rates in Marcus theory is
the  reorganization energy \lambda, defined as

   \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})

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.

-[ Related documentation ]-

* Reorganization energy recipe details

[image]


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: image generated from notebook][image]


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


Related documentation
---------------------

* Reorganization energy recipe details


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