Vibrationally resolved electronic spectra

To calculate vibrational effects on the electronic excitations (Uv/vis, X-ray), one needs to do a frequency calculation both at the ground state as well as the excited state of interest. Next Franck-Condon factors need to be calculated for the transition between the two electronic states, which can be done with the FCF program, described below. These Franck-Condon factor 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.

FCF program: Franck-Condon Factors

fcf is an auxiliary program which can be used to calculate Franck-Condon factors from two vibrational mode calculations [418].

fcf requires an ASCII input file where the user specifies the TAPE21 files from two adf 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. The number of vibrational quanta that have to be taken into account for both states in the evaluation of the Franck-Condon factors have to be specified.

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 too the standard output, including any error messages.

Introduction

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 [419].

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. [420]. 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[420].

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!\]

Input

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

$ADFBIN/fcf << eor
STATES state1 state2
QUANTA l1 l2
eor
STATES state1 state2
The filenames of two TAPE21 files resulting from a numerical or analytical frequency calculation. The calculations must have been performed on the same molecule, i.e. the type, mass and order of occurrence of all the atoms (or fragments) has to be the same in both files.
(optional) MODES first last
The first and last mode to be taken into account in the calculation. If this option is omitted, 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.
(optional) LAMBDA lambda
The minimum value of the electron-phonon coupling for a mode to be taken into account in the calculation. The default value is zero. 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.
QUANTA l1 l2
The maximum number of vibrational quanta to be taken into account for both states. Franck-Condon factors will be calculated for every permutation of up to and including l1/l2 quanta over the vibrational modes.
(optional) TRANSLATE
Move the center of mass of both geometries to the origin.
(optional) ROTATE
Rotate the geometries to maximize the overlap of the nuclear coordinates.
(optional) SPECTRUM freqmin freqmax nfreq

If SPECTRUM is included the vibrational spectrum is calculated. A histogram of the spectrum is calculated for the frequency range that is provided on input. The three parameters that define the frequency range are:

freqmin
minimum frequency for which the spectrum is calculated.
freqmax
maximum frequency for which the spectrum is calculated.
nfreq
number of frequencies for which the spectrum is calculated.

Only a few keys from the TAPE21 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 TAPE21 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 "Geometry%nr of atoms" "Geometry%xyz" "Geometry%nr of atomtypes" \
"Geometry%fragment and atomtype index" "Geometry%atomtype" "Geometry%mass" \
"Freq%Frequencies" "Freq%Normalmodes"

Result: TAPE61

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

contents of TAPE61 comments
firstmode, lastmode the first and last vibrational mode taken into account
lambda the minimum value of the electron-phonon coupling
maxl1, maxl2 maximum level (or maximum number of vibrational quanta) in both states
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 (x0 and x0 ‘)
b0 zero-order axis-switching matrix matrix (B0 )
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 λ’)
maxp1, maxp2 maximum number of permutations of maxl1/maxl2 quanta over the vibrational modes
i0 ground state to ground state overlap integral (I0,0 )
freq1, freq2 frequencies of every permutation of the vibrational quanta for both states
fcf maxp1 by maxp2 Franck-Condon factor matrix
fcfsum1, fcfsum2 average sum of the Franck-Condon factors for both states

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.

Example absorption and fluorescence

In this example it is assumed that the molecule has a singlet ground state S0 , and the interesting excited state is the lowest singlet excited state S1 . First one needs to a the ground state geometry optimization, followed by a frequency calculation. The TAPE21 of the ground state frequency calculation will be called s0.t21. Next one needs to do an excited state geometry optimization. Here it is assumed that the lowest singlet excited state S1 is of interest:

EXCITATION
   Onlysing
   Lowest 1
END
EXCITEDGO
   State A 1
   Singlet
END
GEOMETRY
END

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

EXCITATION
ONLYSING
lowest 1
END
EXCITEDGO
    State A 1
    Singlet
END
GEOMETRY
    frequencies
END
...
mv TAPE21 s1.t21

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 S1 (S1 \(leftarrow\) S0 ), using the FCF program, which calculates the Franck-Condon factors between the vibrational modes of the two electronic states, with input

STATES s0.t21 s1.t21
QUANTA 0 5
SPECTRUM 0 10000 1001
TRANSLATE
ROTATE

The number of vibrational quanta for the excited state should be larger in case of small molecules. 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 S1 to the vibrational levels of the singlet ground state state S0 (S1 → S0 ). Input for the FCF program is in this case:

STATES s0.t21 s1.t21
QUANTA 5 0
SPECTRUM -10000 0 1001
TRANSLATE
ROTATE

The number of vibrational quanta for the ground state should be larger in case of small molecules.

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 S0 and the lowest vibrational state of the electronically singlet excited state S1 .

Example phosphorescence

In this example it is assumed that the molecule has a singlet ground state S0 , and the interesting excited state is the lowest triplet excited state T1 . 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 a the ground state geometry optimization, followed by a frequency calculation. The TAPE21 of the ground state frequency calculation will be called s0.t21. Next one needs to do an excited state geometry optimization of the lowest triplet excited state, followed by a frequency calculation.

CHARGE 0.0 2.0
UNRESTRICTED
GEOMETRY
END
AnalyticalFreq
End
...
mv TAPE21 t1.t21

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

STATES s0.t21 t1.t21
QUANTA 5 0
SPECTRUM -10000 0 1001
TRANSLATE
ROTATE

The number of vibrational quanta for the ground state should be larger in case of small molecules.

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 S0 and the lowest vibrational state of the electronically triplet excited state T1 .

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). In ADF spin-orbit coupling can be treated self-consistently (i.e. non perturbatively) during both the SCF and TDDFT parts of the computation, see the section on excitation energies and spin-orbit coupling.

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

EXCITATION
   Onlytrip
   Lowest 1
END
EXCITEDGO
   State A 1
   Triplet
END
GEOMETRY
END

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

EXCITATION
   Onlytrip
   Lowest 1
END
EXCITEDGO
   State A 1
   Triplet
END
GEOMETRY
   frequencies
END
...
mv TAPE21 t1.t21