AH-FC: Adiabatic Hessian Franck-Condon
**************************************

Electronic spectra, such as absorption, emission, phosphorescence, and
ionization, may contain a vibrational structure.

In the Adiabatic Hessian Franck-Condon (AH-FC) method one needs to do
a frequency calculation both at the ground state as well as the
excited state of interest. The model makes the following assumptions:

* 1. It employs the adiabatic (Born-Oppenheimer) approximation and
     treats the nuclei as moving in an effective potential defined by
     the electronic configuration.

* 2. It works at the Franck-Condon (FC) point and assumes that the
     transition occurs at the ground state equilibrium structure for
     absorption and at the excited state equilibrium structure for
     emission.

* 3. It applies the harmonic approximation to both the ground and
     excited state potential energy surfaces.

Franck-Condon factors can be calculated for the transition between the
two electronic states, which can be done with the FCF Module,
described below. These Franck-Condon factors can then be used to
predict the relative intensities of absorption or emission lines in
the electronic spectra. Note that the Herzberg-Teller effect is not
taken into account.

In the vertical gradient Franck-Condon (VG-FC) method some extra
assumptions are made compared to the AH-FC method. This approach is
particularly effective for large molecules as it shows linear scaling
with the number of normal modes.

* See the documentation of the VG-FC method.


FCF module: Franck-Condon Factors
=================================

*fcf* is an auxiliary program which can be used to calculate Franck-
Condon factors from two vibrational mode calculations [1]. The program
assumes that the starting vibrational wavefunction is in the ground
state, as is the case for most spectroscopic experiments performed at
room temperature.

*fcf* requires an ASCII input file where the user specifies the .rkf
files from two vibrational mode calculations, carried out for two
different electronic, spin or charge states of the same molecule.
These calculations can be either numerical or analytical.

*fcf* produces a (binary) KF file TAPE61, which can be inspected using
the KF utilities. Furthermore, *fcf* writes the frequencies,
vibrational displacements and electron-phonon couplings for both
states to the standard output, including any error messages.

In AMS2022 the algorithm for the calculation of the Franck-Condon
factors has been changed in order to improve the precision, increase
the speed, and lower the memory requirements.


Theory
------

Franck-Condon factors are the squares of the overlap integrals of
vibrational wave functions. Given a transition between two electronic,
spin or charge states, the Franck-Condon factors represent the
probabilities for accompanying vibrational transitions. As such, they
can be used to predict the relative intensities of absorption or
emission lines in spectroscopy or excitation lines in transport
measurements.

When a molecule makes a transition to another state, the equilibrium
position of the nuclei changes, and this will give rise to vibrations.
To determine which vibrational modes will be active, we first have to
express the displacement of the nuclei in the normal modes:

   \mathbf{k} =  \mathbf{L'}^T \mathbf{m}^{1/2} (\mathbf{B}_0
   \mathbf{x}_0 - \mathbf{x'}_0)

Here, \mathbf{k} is the displacement vector, \mathbf{L} is the normal
mode matrix, \mathbf{m} is a matrix with the mass of the nuclei on the
diagonal, \mathbf{B} is the zero-order axis-switching matrix and
\mathbf{x}_0  is the equilibrium position of the nuclei. For free
molecules, depending on symmetry constraints, the geometries of both
states may have been translated and/or rotated with respect to each
other. To remove the six translational and rotational degrees of
freedom, we can center the equilibrium positions around the center of
mass and rotate one of the states to provide maximum overlap. The
latter is included with the zero-order axis-switching matrix
\mathbf{B}, implemented according to [2].

When we have obtained the displacement vector, it is trivial to
calculate the dimensionless electron-phonon couplings. They are given
by:

   \mathbf{\lambda} = (\mathbf{\Gamma}/2)^{1/2} \mathbf{k}

Here, \mathbf{\Gamma} = 2 \pi \mathbf{\omega} / h is a vector
containing the reduced frequencies. [3]. The Huang-Rhys factor
\mathbf{g} is related to \mathbf{\lambda} as:

   \mathbf{g} = \mathbf{\lambda}^2

The reorganization energy per mode is

   \mathbf{E} = (h / 2 \pi) * \mathbf{\omega} \mathbf{\lambda}^2

When the displacement vector \mathbf{k}, the reduced frequencies
\mathbf{\Gamma} and \mathbf{\Gamma}', and the Duschinsky rotation
matrix \mathbf{J} = \mathbf{L'}^T \mathbf{B}_0 \mathbf{L} have been
obtained, the Franck-Condon factors can be calculated using the two-
dimensional array method of Ruhoff and Ratner [3].

There is one Franck-Condon factor for every permutation of the
vibrational quanta over both states. Since they represent transition
probabilities, all Franck-Condon factors of one state which respect to
one vibrational state of the other state must sum to one. Since the
total number of possible vibrational quanta, and hence the total
number of permutations, is infinite, in practice we will calculate the
Franck-Condon factors until those sums are sufficiently close to one.
Since the number of permutations rapidly increases with increasing
number of vibrational quanta, it is generally possible to already stop
after the sum is greater than about two thirds. The remaining one
third will be distributed over so many Franck-Condon factors that
their individual contributions are negligible.

In the limiting case of one vibrational mode, with the same frequency
in both states, the expression for the Franck-Condon factors of
transitions from the ground vibrational state to an excited
vibrational state are given by the familiar expression:

   |l_{0,n}|^2 = e^{-\lambda^2} \lambda^{2n} / n!


Algorithm
---------

The algorithm for the calculation of the FC integrals is based on the
work by Santoro et al. [4] which is here summarized. There are in
principle infinite FC factors that should be computed in order to
achieve 100% convergence in the spectrum. This number is reduced by
dividing the vibrations into classes.

A class is defined based on the number of simultaneously excited
vibrational modes, thus a class 3 state has exactly 3 modes which are
simultaneously excited, while all the remaining modes have zero
quanta. The *fcf* program first evaluates the class 1 and 2 FC up to a
fairly large vibrational quantum number and uses these results to
infer how much each mode should be excited for the calculation of the
integrals involving states of higher classes. We refer to the original
paper for the details of the algorithm.

The convergence of the results after each class can be estimated by
summing all FC factors, which should converge to 1 in the limit of a
complete set. The program is terminated when the maximum allowed class
is reached, when there is little gain from one class to the next, or
when a sufficient convergence criterion is met. It should be noted
that if one is interested in the overall shape of the spectrum then in
most cases it isn't necessary to reach a convergence that is close to
100%.


Input
-----

The input for *fcf* is keyword-oriented and is read from the standard
input. *fcf* recognizes several keywords, but only the states have to
be specified to perform the calculation. All input therefore contains
at least the following lines:

   $AMSBIN/fcf << eor
      STATE1 state1
      STATE2 state2
   eor

The spectrum is always calculated from the ground vibrational state
belonging to the potential energy surface specified by *state1* to the
vibronic states belonging to *state2*. Documentation for all keys of
*fcf*:

   $AMSBIN/fcf << eor
      State1 state1
      State2 state2
      Lambda float
      Modes first last
      Rotate [True | False]
      Translate [True | False]
      NumericalQuality ["Basic" | "Normal" | "Good" | "Excellent"]
      Prescreening
         Class1   integer
         Class2   integer
         MaxClass integer
         MaxFCI   integer
      End
      Spectrum FreqMin FreqMax NFreq
         LineShape ["Stick" | "Gaussian" | "Lorentzian" ]
         SpcMin float
         SpcMax float
         SpcLen integer
         HWHM   float
      End
      Printing
         Verbose
         FCFThresh float
      End
   eor

"Convergence"
   Type:
      Float

   Description:
      Minimum fraction of the FC factors sum to be calculated

"Lambda"
   Type:
      Float

   Default value:
      0.0

   GUI name:
      Minimum coupling

   Description:
      Optional. The minimum value of the electron-phonon coupling for
      a mode to be taken into account in the calculation. Together
      with the MODES option, this provides a way to significantly
      reduce the total number of Franck-Condon factors. As with the
      MODES option, always check if the results do not change too
      much.

"Modes"
   Type:
      Integer List

   Default value:
      [0, 0]

   Description:
      Optional. The first and last mode to be taken into account in
      the calculation. By default, all modes are taken into account.
      This option can be used to effectively specify and energy range
      for the Franck-Condon factors. When using this options, always
      check if the results (electron-phonon couplings, ground state to
      ground overlap integral, average sum of Franck-Condon factors,
      etc.) do not change too much.

"NumericalQuality"
   Type:
      Multiple Choice

   Default value:
      Normal

   Options:
      [0-0, Basic, Normal, Good, Excellent]

   Description:
      Set the quality of several technical aspects of an FC
      calculation, including details for the prescreening, the amount
      of computed integrals, the maximum integral class evaluated, as
      well as the convergence criteria.

"Printing"
   Type:
      Block

   Description:
      Optional. Printing Options.

   "FCFThresh"
      Type:
         Float

      Default value:
         0.01

      GUI name:
         FC factor threshold

      Description:
         The minimum value for the printing of a FC factor.

   "Verbose"
      Type:
         Bool

      Default value:
         No

      GUI name:
         Verbose printing

      Description:
         Increase output verbosity.

"Rotate"
   Type:
      Bool

   Default value:
      Yes

   Description:
      Recommended to be used. Rotate the geometries to maximize the
      overlap of the nuclear coordinates.

"Spectrum"
   Type:
      Block

   Description:
      Optional. Controls the details of the spectrum calculation

   "HWHM"
      Type:
         Float

      Default value:
         100.0

      Description:
         Half-Width at Half-Maximum for the spectral lineshape
         function

   "LineShape"
      Type:
         Multiple Choice

      Default value:
         Stick

      Options:
         [Stick, Gaussian, Lorentzian]

      Description:
         Type of lineshape function

   "SpcLen"
      Type:
         Integer

      Default value:
         2001

      Description:
         Number of points in the spectrum

   "SpcMax"
      Type:
         Float

      Default value:
         9500.0

      Description:
         Maximum absolute energy difference between the states
         relative to the 0-0 transition

   "SpcMin"
      Type:
         Float

      Default value:
         -500.0

      Description:
         Minimum absolute energy difference between the states
         relative to the 0-0 transition

   "Type"
      Type:
         Multiple Choice

      Default value:
         Absorption

      Options:
         [Absorption, Emission]

      Description:
         Type of spectrum

"State1"
   Type:
      String

   Description:
      Filename of the result files of a numerical or analytical
      frequency calculation for the first (initial) state.

"State2"
   Type:
      String

   Description:
      Filename of the result files of a numerical or analytical
      frequency calculation for the second (final) state.

"Translate"
   Type:
      Bool

   Default value:
      Yes

   Description:
      Recommended to be used. Move the center of mass of both
      geometries to the origin.

Only a few keys from the .rkf file are used for the calculation of the
Franck-Condon factors. Disk space usage can be significantly reduced
by extracting just these keys from the .rkf file before further
analysis. The following shell script will extract the keys from the KF
file specified by the first argument and store them in a new KF file
specified by the second argument using the *cpkf* utility:

   #!/bin/sh
   cpkf $1 $2 General Molecule Vibrations


Result: TAPE61
--------------

After a successful calculation, *fcf* produces a TAPE61 KF file. All
results are stored in the Fcf section:

 
| **contents of TAPE61**                                                                                                                                                                                                                                     | **comments**                                                                                                                                                                                                                                                                                                                                                   |
 
| lambda                                                                                                                                                                                                                                                     | the minimum value of the electron-phonon coupling                                                                                                                                                                                                                                                                                                              |
 
| translate, rotate                                                                                                                                                                                                                                          | whether the TRANSLATE and ROTATE options were specified in the input                                                                                                                                                                                                                                                                                           |
 
| natoms                                                                                                                                                                                                                                                     | number of atoms in the molecule                                                                                                                                                                                                                                                                                                                                |
 
| mass                                                                                                                                                                                                                                                       | atomic mass vector (**m**)                                                                                                                                                                                                                                                                                                                                     |
 
| xyz1, xyz2                                                                                                                                                                                                                                                 | equilibrium geometries of both states (**x**_0  and **x**_0 ')                                                                                                                                                                                                                                                                                                 |
 
| b0                                                                                                                                                                                                                                                         | zero-order axis-switching matrix matrix (**B**_0 )                                                                                                                                                                                                                                                                                                             |
 
| nmodes                                                                                                                                                                                                                                                     | number of vibrational modes with a non-zero frequency                                                                                                                                                                                                                                                                                                          |
 
| gamma1, gamma2                                                                                                                                                                                                                                             | reduced frequencies of both states (Γ and Γ')                                                                                                                                                                                                                                                                                                                  |
 
| lmat1, lmat2                                                                                                                                                                                                                                               | mass-weighted normal modes of both states (**L** and **L**')                                                                                                                                                                                                                                                                                                   |
 
| jmat                                                                                                                                                                                                                                                       | Duschinsky rotation matrix (**J**)                                                                                                                                                                                                                                                                                                                             |
 
| kvec1, kvec2                                                                                                                                                                                                                                               | displacement vectors for both states (**k** and **k'**, kvec1 is used for the calculation of the Franck-Condon factors)                                                                                                                                                                                                                                        |
 
| lambda1, lambda2                                                                                                                                                                                                                                           | electron-phonon couplings for both states (λ and λ')                                                                                                                                                                                                                                                                                                           |
 
| i0                                                                                                                                                                                                                                                         | ground state to ground state overlap integral (I_0,0 )                                                                                                                                                                                                                                                                                                         |
 

In addition to producing a binary TAPE61 file, *fcf* also writes the
frequencies, displacement vectors and electron-phonon couplings for
both states to the standard output.


FCF example absorption and fluorescence
=======================================

In this example it is assumed that the molecule has a singlet ground
state S_0 , and the interesting excited state is the lowest singlet
excited state S_1 . First one needs to do a ground state geometry
optimization, followed by a frequency calculation. For the ground
state frequency calculation we will use AMS_JOBNAME=S0. Next one needs
to do an excited state geometry optimization. Here it is assumed that
the lowest singlet excited state S_1  is of interest:

   AMS_JOBNAME=S1_GEO $AMSBIN/ams <<eor
      ...
      Task GeometryOptimization
      Engine ADF
         ...
         Excitation
            Onlysing
            Lowest 1
         End
         ExcitedGO
            State A 1
            Singlet
         End
      EndEngine
   eor

To get the frequencies for this excited state, numerical frequencies
need to be calculated, at the optimized geometry of the first excited
state:

   AMS_JOBNAME=S1 $AMSBIN/ams <<eor
      ...
      LoadSystem
         File S1_GEO.results/adf.rkf
      End
      Task SinglePoint
      Properties
         NormalModes True
      End
      Engine ADF
         ...
         Excitation
            Onlysing
            Lowest 1
         End
         ExcitedGO
            State A 1
            Singlet
         End
      EndEngine
   eor

Next for the absorption spectrum, we look at excitations from the
lowest vibrational state of the electronic ground state to the
vibrational levels of the first singlet excited state S_1  (S_1
\leftarrow  S_0 ), using the   FCF program, which calculates the
Franck-Condon factors between the vibrational modes of the two
electronic states, with input

   $AMSBIN/fcf << eor
      STATE1 S0.results/adf.rkf
      STATE2 S1.results/adf.rkf
      TRANSLATE
      ROTATE
   eor

See  the description of FCF program for more details.

For the fluorescence spectrum, we look at excitations from the lowest
vibrational state of the first singlet excited state  S_1  to the
vibrational levels of the singlet ground state state  S_0  (S_1  →
S_0 ). Input for the  FCF program is in this case:

   $AMSBIN/fcf << eor
      STATE1 S1.results/adf.rkf
      STATE2 S0.results/adf.rkf
      TRANSLATE
      ROTATE
   eor

Note that the FCF program calculates the spectrum relative to the 0-0
transition. Thus one should add to spectrum calculated with FCF the
difference in energy of the lowest vibrational state of the ground
state S_0  and the lowest vibrational state of the electronically
singlet excited state S_1 .


FCF Example phosphorescence
===========================

In this example it is assumed that the molecule has a singlet ground
state S_0 , and the interesting excited state is the lowest triplet
excited state T_1 . Emission from a triplet state to a singlet state
is spin forbidden, however, due to spin-orbit coupling such
transitions may occur. In the following we assume that the geometry of
the triplet excited state is not influenced much by spin-orbit
coupling.

First one needs to do a ground state geometry optimization, followed
by a frequency calculation, using AMS_JOBNAME=S0. Next one needs to do
an excited state geometry optimization of the lowest triplet excited
state, followed by a frequency calculation.

   AMS_JOBNAME=T1 $AMSBIN/ams <<eor
      ...
      Task GeometryOptimization
      Properties
         NormalModes True
      End
      Engine ADF
         Unrestricted
         SpinPolarization 2.0
      EndEngine
   eor

For the phosphorescence spectrum, we look at excitations from the
lowest vibrational state of the first triplet excited state  T_1  to
the vibrational levels of the singlet ground state state  S_0  (T_1  →
S_0 ). Input for the  FCF program is in this case:

   $AMSBIN/fcf << eor
      STATE1 T1.results/adf.rkf
      STATE2 S0.results/adf.rkf
      TRANSLATE
      ROTATE
   eor

Note that the FCF program calculates the spectrum relative to the 0-0
transition. Thus one should add to spectrum calculated with FCF the
difference in energy of the lowest vibrational state of the ground
state S_0  and the lowest vibrational state of the electronically
triplet excited state T_1 .

Zero field splitting (ZFS) and the radiative rate constants (i.e.
radiative phosphorescence lifetimes) could be calculated  with spin-
orbit coupled ZORA time-dependent density functional theory (ZORA-
TDDFT). With the ADF engine spin-orbit coupling can be treated self-
consistently (i.e. non perturbatively) during both the SCF and TDDFT
parts of the computation.

An alternative to the use of the unrestricted formalism to calculate
the lowest triplet excited state is to use the TDDFT formalism:

   AMS_JOBNAME=T1_GEO $AMSBIN/ams <<eor
      ...
      Task GeometryOptimization
      Engine ADF
         Excitation
            Onlytrip
            Lowest 1
         End
         ExcitedGO
            State A 1
            Triplet
         End
      EndEngine
   eor

To get the frequencies for this excited state, numerical frequencies
need to be calculated, at the optimized geometry of the first excited
state.

   AMS_JOBNAME=T1 $AMSBIN/ams <<eor
      ...
      LoadSystem
         File T1_GEO.results/adf.rkf
      End
      Task SinglePoint
      Properties
         NormalModes True
      End
      Engine ADF
         ...
         Excitation
            Onlytrip
            Lowest 1
         End
         ExcitedGO
            State A 1
            Triplet
         End
      EndEngine
   eor


Density of States computed with AH-FC method
============================================

The FCF module can also be used to evaluate the density of states
*rho* that enters the expression of the rate constant for a number of
transfer phenomena based on Fermi's Golden Rule:

   \mathbf{k} =  \frac{2\pi}{\hbar} J^2 \rho

To calculate the transfer rate from a donor molecule to an acceptor
molecule it is sufficient to employ the FCF code as described above to
evaluate the emission spectrum of the donor and the absorption
spectrum of the acceptor. In fact, it can be shown [5] that the
expression for the density of states can be written as:

   \rho = \frac{1}{\hbar}\int_{-\infty}^{+\infty}F^\mathrm{D}(\omega)
   F^\mathrm{A}(\omega) d \omega

Where *F* denotes the Franck-Condon spectrum calculated by the FCF
code, specifically the absorption for the acceptor and the emission
for the donor.

The convolution may be easily calculated by using the PLAMS library.

[1] J.S. Seldenthuis, H.S.J. van der Zant, M.A. Ratner and J.M.
    Thijssen, *Vibrational Excitations in Weakly Coupled Single-
    Molecule Junctions: A Computational Analysis*, ACS Nano 2, 1445
    (2008)

[2] G.M. Sando and K.G. Spears, *Ab Initio Computation of the
    Duschinsky Mixing of Vibrations and Nonlinear Effects*, Journal of
    Physical Chemistry A 105, 5326 (2001)

[3] P.T. Ruhoff and M.A. Ratner, *Algorithms for computing Franck-
    Condon overlap integrals*, International Journal of Quantum
    Chemistry 77, 383 (2000)

[4] F. Santoro, R. Improta, A. Lami, J. Bloino, V. Barone, *Effective
    method to compute Franck-Condon integrals for optical spectra of
    large molecules in solution*, J. Chem. Phys. 126, 084509 (2007)

[5] X. de Vries *High energy acceptor states strongly enhance exciton
    transfer between metal organic phosphorescent dyes*, Nature
    Communications 11, 1292 (2020)
