Download The ezSpectra Suite: An easy-to-use Toolkit for Spectroscopy ... and more Exercises Chemistry in PDF only on Docsity!
The ezSpectra Suite: An easy-to-use Toolkit for
Spectroscopy Modeling
Samer Gozema^ and Anna I. Krylovb
a (^) Department of Chemistry, Georgia State University, Atlanta, GA 30303, USA
b (^) Department of Chemistry, University of Southern California, Los Angeles, CA 90089, USA
March 1, 2021
Abstract A molecule’s spectrum encodes information about its structure and electronic proper- ties. It is a unique fingerprint that can serve as a molecular ID. Quantum chemistry cal- culations provide key ingredients for interpreting spectra, but modeling the spectra rarely ends there; it requires additional steps that entail combined treatments of electronic and nuclear degrees of freedom and account for specifics of the experimental setup (light energy, polarization, averaging over molecular orientations, temperature, etc.). This Software Fo- cus article describes the ezSpectra suite, which currently comprises two stand-alone open- source codes: ezFCF and ezDyson. ezFCF calculates Franck-Condon factors, which yield vibrational progressions for polyatomic molecules within the double-harmonic approxi- mation. ezDyson calculates absolute cross-sections for photodetachment/photoionization processes and photoelectron angular distributions using Dyson orbitals computed by a quantum chemistry program.
1 Introduction
“The spectrum is a message from the molecule” – Unknown
Spectroscopy is the most powerful tool for probing matter. Despite the mind-boggling sophistication of modern light sources and spectroscopic techniques, conceptually it remains just as simple as in the time of Newton’s experiments: we shine light at a sample and back comes the spectrum with information about the internal structure of the sample. Spectroscopy allows us to communicate with molecules. Each spectrum is a message from the molecule, reporting on its structure. However, the message is coded. To translate an in- tricate pattern of spectral lines and intensities into what the nuclei and electrons are doing, we need a cipher. Theory provides the key to decipher the message (see Fig. 1). By using quantum chemical calculations, we can obtain essential parameters of the molecular structure. Using these parameters as an input, we can then simulate the spectra. By comparing the com- puted and experimentally measured spectra, we can confirm molecular identity, assign spectral features, and relate the spectra to such concepts as equilibrium bond lengths, vibrational fre- quencies, molecular orbitals, and so on.
Figure 1: Spectroscopic measurements (red panel) yield encrypted messages, which can be decoded using electronic structure calculations (blue panels) to reveal detailed information about molecular structure (black panel).
Different types of spectroscopies probe different types of transitions and report on different molecular properties. To appreciate the scope and variety of modern spectroscopy, consider the vast range of energies and intensities of light used,^1 in addition to nonlinear,^2 angle-resolved,3, 4
Hamiltonian (full molecular Hamiltonian minus the nuclear kinetic energy operator), such that they depend only parametrically on the nuclear coordinates (as signified by ’;’ instead of ’,’):
Pi′f ′′ ∝
ψi(r; R)χi′ (R) ˆμ ψf (r; R)χf ′′ (R) dr dR
Here indices i′^ and f ′′^ denote the vibrational states of the two electronic states (which, within the Born-Oppenheimer approximation, are determined by the nuclear problem with the potential defined by the electronic Schr¨odinger equation). Integration over the electronic coordinates yields the electronic transition dipole moment for the i → f transition
μif (R) =
ψi(r; R) ˆμ ψf (r; R) dr, (3)
so that the probability of transition can be written as
Pi′f ′′ ∝
χi′ (R) μif (R) χf ′′ (R) dR
Eq. (4) is the starting point for understanding the spectrum. It contains electronic factor μif (R), which encodes the information about the electronic and vibrational wave-functions of the two states. Within the Condon approximation,^11 it is assumed that μif (R) depends weakly on the nuclear coordinates and, therefore, can be evaluated at a fixed nuclear geometry (e.g., at the equilibrium geometry Re of the initial state). In this case, the probability of the transition simply becomes: Pi′f ′′ ∝ μif (Re)^2
χi′ (R) χf ′′ (R) dR
The overlap integrals between the nuclear wave-functions of the two states are called the Franck- Condon factors (FCFs)11–13—their squares are the intensities due to vibrational progressions (sometimes the intensities are called FCFs). The stick spectrum can then be assembled from the individual FCFs as follows:
I(ω) =
i′,f ′′
Pi′f ′′ · δ(ω − Ei′f ′′ ), (6)
where δ is the Kroneker δ-function and Ei′f ′′ is the energy difference between the two vibronic
states. The macroscopic observables, such as cross sections, Einstein coefficients, etc, can be computed from I(ω) (see, for example, Ref. 14).
Figure 2: Factors determining the spectral line shape for electronic transitions. The ground state (red), the first excited state (blue), and the second excited state (green) are represented by one-dimensional harmonic potentials along a normal mode Q. The change in the excited- state equilibrium geometry relative to the ground state is denoted by ∆Q′^ and ∆Q′′. The main progression is determined by the overlap of the initial state ν = 0 wave function (red translucent area curve) and each of the excited-state vibrational wave functions (blue translu- cent area curves). If T 6 = 0, hot bands may appear, such as peaks resulting from the ν = 1 transitions indicated by dotted lines on the right. The total spectrum is an overlay of the main progression and hot bands. Instrumental broadening and interactions with the environment (e.g., inhomogeneous broadening in condensed phase) results in a broadened spectrum like the one shown by black dashed lines on the right. The definition of electronic vertical (Ev′ ) and adiabatic (Ea ee′) excitation energies are shown for the first excited state.
Fig. 2 illustrates how Eqs. (5) and (6) give rise to the spectral profile due to transitions from the initial state (red) to the two final states (blue and green). The μif for each state determines the magnitude of the electronic contribution to the cross-section, and the FCFs determine the probability of transition to different vibrational levels within each electronic state. As per Eq. (5), FCFs are simply the overlap integrals between the initial-state vibrational wave function (red translucent curve in Fig. 2) and the final-state vibrational wave functions (blue translucent curves). A large overlap leads to a large FCF, giving rise to a high probability of transition to that vibrational level. For two states represented by two identical one-dimensional harmonic potentials, only the diagonal transitions (∆ν=0) are allowed by virtue of the orthogonality of the Hermite polynomials, leading to a spectrum with no vibrational progression (e.g., for
molecules show that the maximum intensity often corresponds to the 00 transition, such that λmax represents the adiabatic excitation energy^15 (more on this below). Hence, the most reliable way to compare theory and experiment is to complement quantum-chemistry calculations of the electronic energy differences by calculation of the FCFs. The above treatment describes both the transitions between the electronically bound states (i.e., the ground and excited states of a molecule) and the transitions to the continuum, when the energy of light (Ehν ) is sufficient to eject an electron. In the latter case, the spectra are measured as the yield of photo-electrons at different kinetic energy (Ek), as per the following energy balance: Ek = Ehν − ∆Ei′f ′′. (7)
In the experimental community, it is common to distinguish between photoionization (electrons ejected from a neutral molecule) and photodetachment (electrons ejected from an anion), but, for the sake of brevity, we will use the term photoionization to refer to both. For the transitions between bound electronic states, Eq. (5) yields the total absorption cross section, and μif is an electronic transition dipole moment element. For the transitions into the continuum, the probability of ejecting an electron is also determined by Eq. (5), but μif is replaced by the photoelectron dipole matrix element μi∞:^16
μi∞(R) =
ψni (r 1 , r 2 ,... , rn; R) ˆμ Aˆ [ψ fn −^1 (r 1 , r 2 ,... , rn− 1 ; R)ψkel (rn)]^ dr 1 dr 2... drn, (8)
where n is the number of electrons in the initial state, Aˆ is an anti-symmetrizing operator, ψ fn −^1 is the wave-function of the ionized (detached) target, and ψkel is the wave-function of the ejected electron with wave vector k. The ezSpectra suite currently comprises two stand-alone open-source codes (distributed un- der the GPL license): ezFCF and ezDyson. The former computes FCFs (the overlap integrals in Eq. (5)), and the latter computes properties related to the photoelectron dipole matrix element in Eq. (8). Both are written in C++ and use a built-in parser to read xml-formatted input files. The source code, precompiled executables, manuals, and examples of input and out- put files are available for download at http://iopenshell.usc.edu/downloads/.^17 We also provide Python scripts for automatic preparation of the xml input files from raw outputs of quantum-
chemistry programs. The scripts currently support Q-Chem, Cfour/ACESII, Molpro, ORCA, GAMESS, and Gaussian outputs.
2 ezFCF
2.1 Software capabilities and features
Figure 3: ezFCF input parameters, options, and output. The grey box indicates selected features implemented in the code. The output contains information listed in yellow box.
Overview: ezFCF, formerly known as ezSpectrum, computes vibrational progressions (stick spectra) for electronic excitation or photoelectron spectra for polyatomic molecules within the double-harmonic approximation. Required input: ezFCF requires equilibrium geometries, harmonic frequencies, and normal- mode vectors for each electronic state. These quantities can be computed using an ab initio electronic structure package (Fig. 3). Additionally, the user should specify the temperature (which depends on experimental conditions) and the adiabatic excitation or ionization energy (which can be either computed or taken from the experiment). What the code computes: ezFCF computes FCFs, the overlaps between the multi-dimensional vibrational wave functions of the initial (usually, ground) state and the final (excited or ionized) states, as per Eq. (5). In the current version, both the initial and the final states are described
(b)
{ ”}q
{ ’}q
{ ”}q
{ ’}q
{ ”}^ { ’}qq
{ ”}q
{ ’}q
{ ”}q
{ ’}q
{ ”}q
{ ’}q
(a) (c)
Q’
Q” Q’
Q”
Q’
Q” Q’
Q”
Q’
Q”
Q’
Q”
Figure 4: The effect of rotations of the normal coordinates on FCFs within the parallel-mode approximation. (a) The correct overlap between wave functions of the lower (Q′) and upper (Q′′) surfaces. (b) The overlap when the lower normal coordinates are rotated to coincide with the upper coordinates. (c) The overlap when the upper normal coordinates are rotated to coincide with the lower coordinates. Reproduced with permission from Ref. 21.
The validity of the parallel-mode approximation (and the alignment of the normal modes) can be assessed by checking the normal-modes overlap matrix, which should be close to diagonal. The order of the normal modes can be flipped, which is manifested by near-unity off-diagonal values in the overlap matrix; this does not mean that the parallel mode approximation is invalid, as the diagonal form of the overlap matrix can be easily restored by the manual reordering of the normal modes to match. If the normal modes are rotated (i.e., when the values of the diagonal elements are significantly less than one, even after the modes are properly ordered), it is recommended to use the normal modes of the target state instead of the normal modes of the initial state. As explained in Fig. 4, the initial vibrational state (ν 1 ′ = ν 2 ′... = ν N′ = 0), which gives rise to the main vibrational progressions, has a fully symmetric wave-function. Consequently, projecting such a state into the rotated coordinate system (normal modes of the target state Q′′) does not affect its nodal structure and, therefore, introduces relatively small errors in the overlaps. In contrast, projecting the vibrationally excited state into another coordinate system can significantly affect the overlaps because of the incorrect representation of the nodal structure of the target state. To further reduce the cost of the calculations, the maximum number of vibrational exci- tations in the target state are truncated at a user-defined threshold (this is justified because the overlaps for highly excited states eventually become small) and by excluding some modes (such as those that have small ∆Q and small change in frequencies) from the calculation. The
number of quanta in the initial state depends on the temperature: for T =0 K, ν max′ =0, but for finite T more states need to be included to describe hot bands. The populations of higher vibrational states are given by the Boltzmann factors. Fig. 2 illustrates what is computed in ezFCF. If there are multiple final states (as in the figure), ezFCF reads the input parameters for each state and computes the overall spectrum by adding up the FCFs for all states, properly shifted by the respective excitation or ionization energies. Using the equilibrium geometries, harmonic frequencies, and normal mode vectors of the initial and final state(s) provided by ab initio calculations, ezFCF first aligns the equilibrium geometry of the initial and final states to maximize the overlap of the normal modes. The displacements between the initial state and each final state (∆Qs in Fig. 2) along each normal mode are then determined, using either the normal modes of the initial or the final state, as specified by the user (see Fig. 4). If T =0 K, only excitations from the ν = 0 vibrational level of the initial state are considered. At finite temperatures, the thermal (Boltzmann) distribution is used to populate higher vibrational levels of the initial state (up to νmax) and hot bands are computed as well. ezFCF prints the energies and the FCFs (and their squares) for each vibronic transition, so that the result can be plotted as a stick spectrum, as per Eq. (6). The overall spectrum (printed in the ezFCF output) contains the superposition of the transitions from ν = 0 and all hot bands (scaled by Boltzmann factors), as shown in Fig. 2, and can be directly compared with the vibrational progressions in experimental spectra. The details of the algorithm, which uses a recurrence relation to evaluate the overlap in- tegrals accounting for Duschinsky rotations,^19 is described in the ezFCF manual and in more details in Refs. 22 and 23. Such calculations are more expensive than parallel-mode approxi- mation and scale unfavorably with the system size. The cost can be controlled by limiting the number of Franck-Condon active modes and vibrational excitation levels in each mode. We note that double-harmonic approximation fails in systems with large-amplitude soft modes, such as molecular clusters; in such cases more sophisticated treatments are needed.^24
2.2 Examples
To illustrate the theory, we consider several representative examples from the literature. Fig. 5 compares the experimental photoelectron spectrum of the phenolate anion with the stick
action spectrum of the HBDI anion.^28 In the condensed phase, the spectra are broadened due to the interactions with the solvent (in addition to instrumental broadening); yet, they often retain discernible features of the underlying vibrational structure. An example, shown in Fig. 6, is flavin. The excitation and emission spectra of flavin were measured at high resolution by Dick and co-workers for lumiflavin in helium nanodroplets (Fig. 6A).^29 The accompanying calculations reproduced reasonably well the Franck-Condon progression and ascribed the strongest progression-forming mode to an in-plane bending vibration.^29 The UV/Vis spectrum in aqueous solution is broad and featureless, but in many proteins the spectra are better resolved and show three distinct bands in the S 0 →S 1 transition (black spectrum in Fig. 6B). The computed spectrum (FCFs broadened by an empirically fit Gaussians, red spectra in Fig. 6B) reproduces these features and assigns them to the main progression-forming mode.^30 The spectrum features two intense vibronic bands at ca. 2.6 eV and 2.8 eV, as well as a third, less-intense peak at ca. 2.95 eV. Neither of the two more intense bands is close to the vertical one; they are shifted by about 0.15 -0.40 eV relative to Ev.^30 The magnitude of the zero-point energy correction, -0.065 eV, is also noticeable. The large difference between Ev^ and a properly computed λmax has been pointed out in several studies on flavins using various levels of theory.30– These examples illustrate that a common practice of interpreting the band maxima in photoelectron or UV/Vis absorption spectra as vertical ionization or excitation energies is not justified for polyatomic molecules, in which λmax often corresponds to the 00 transition, and is, therefore, much closer to the adiabatic ionization/excitation energies. In molecules discussed here, the adiabatic and vertical energy differ by as much as 0.1-0.4 eV, which is clearly important for quantitative comparisons between theory and experiment. The reason that the conclusions based on the one-dimensional model (often referred to as the Franck- Condon principle) do not hold up in polyatomic molecules is that small structural relaxation along many degrees of freedom adds up to a large energy value.15, 26^ Hence, despite a large energetic difference between Ev^ and Ea, the actual displacements along Franck-Condon active modes are sufficiently small, so that the 00 transition remains the most intense one. We note that the smallness of structural changes can be attributed to the delocalized nature of molecular orbitals in polyatomic molecules.^15
Figure 6: A. Experimental and computed (stick spectra) Franck-Condon pattern of the exci- tation and emission spectra of flavin.^29 B. Experimental (black, for flavin mononucleotide in iLOV^33 ) and computed (red) spectra for flavin in condensed phase. The FCFs were computed for a gas-phase lumiflavin model using ezFCF with B3LYP/cc-pVTZ (red stick spectra and solid red line broadened spectrum) and by Davari et. al using a slightly different approach (dashed red line).^33 The adiabatic and vertical excitation energies are indicated with vertical dashed lines. Reproduced with permission from Refs. 29 and 33.
Figure 7: A schematic illustration of FCFs revealing the differences between excited-state meth- ods that are not captured by vertical energy calculations. Different excited-state geometries (compare panels A and B) give rise to different FCF intensity pattern and band position. Dif- ferent excited-state curvatures (compare panels A and C) give rise to different spacings between the vibrational peaks.
Figure 8: ezDyson input parameters, options, and output. The grey box indicates options and features implemented in the code, which the user controls. The main information printed in the output are in yellow.
be computed. The user also should specify whether the calculation is for randomly oriented molecules or for a molecule aligned and fixed in space. In the latter case, the user should indicate the direction of the laser polarization in the molecular frame. The quality of the calculations is sensitive to the treatment of the photoelectron wave function. For absolute cross sections, correct spin and orbital degeneracy factors also need to be provided by the user. What the code computes: Absolute photoionization cross sections, PADs, and the β anisotropy parameter, which are all functions of the photoelectron dipole matrix element from Eq. (8). In what follows we briefly introduce these quantities and explain how they are computed in ezDyson. Eq. (8) defines the photoelectron dipole matrix element μi∞ as an integral over coordinates of all n electrons (represented by r). Because electrons are indistinguishable and assuming strong orthogonality (i.e., that the photoelectron wave function is orthogonal to all orbitals of the molecule), the integration can be done in two steps, giving rise to
μi∞ =
φd(r 1 )μψkel (r 1 )dr 1 ≡ 〈φd|μ|ψkel 〉, (12)
where φd^ is obtained by integrating the initial n-electron and final n−1-electron wave functions:
φd(r 1 ) =
N
ψni (r 1 ,... , rn)ψn f −^1 (r 2 ,... , rn)dr 2... drn. (13)
For the sake of brevity, here we have dropped the implicit dependence on nuclear coordinates R. φd(r), which is a function of only one electron, contains all the information about the two many-body wave functions ψin and ψn f −^1 needed to evaluate μi∞; it is called a Dyson orbital or generalized overlap.38, 39^ Within the Koopmans’ framework, φd^ is simply the Hartree-Fock canonical orbital vacated by the electron upon ionization. When computed using correlated wave functions (such as EOM-CC,40–43^ configuration interaction,^44 Green’s functions and propa- gator methods,45, 46^ CASSCF,^47 and CASPT2^48 ), Dyson orbitals account for electron correlation and orbital relaxation effects that are missing in Koopmans’ picture of ionization while retaining the conceptual simplicity of a one-electron picture. Among numerous correlated approaches,
The plot of β as a function of the energy of ionizing radiation (or photoelectron kinetic energy) encodes the information about the shape (i.e., angular momentum and the diffuseness) of the orbital from which the electron is ejected. Fig. 9 shows the β values for H−^ and F−^ as a function of the photoelectron kinetic energy computed with ezDyson. While β is constant for H−^ at all energies (consistent with a p-wave that is parallel to the polarization of light), F−^ has a more complex energy-dependent β profile. At low energies, ionization of fluoride’s p-orbital results in a photoelectron wave function dominated by an s-wave (and, therefore, isotropic). In contrast, at high energies, the photoelectron wave function is dominated by a d -wave. At intermediate energies, cross-terms arising from the combination of s-waves and d -wave result in a perpendicular distribution of electrons, with β = −1 at some point.^63
Figure 9: β values for H−^ and F−^ as a function of energy of the photoelectron kinetic energy.
The absolute photoionization cross section is a measure of the total probability of an electron being ejected in any direction. It is obtained by integrating the differential cross section over all solid angles (Ωk). The calculation of energy-dependent photoionization cross sections is useful for assigning the electronic origins of photoelectron spectra (e.g., to identify contributions from direct ionization and auto-ionizing metastable states). The ability to accurately compute absolute cross sections is also important in applications of photoelectron spectroscopy as an analytic tool for quantitative chemical analysis.^64 The differential cross section is computed from μi∞ using16, 65 dσ dΩk^ =
4 π^2 kE c 〈φ
d|μ|ψelk 〉 (^2) =^4 π^2 kE c |μi∞(θ, φ)|
where k is the magnitude of the photoelectron wave vector k, E is the energy of the ionizing radiation, and c is the speed of light. θ and φ are the polar and azimuthal angles of the photoelectron relative to the polarization of light. Eq. (15) is the starting point for computing total cross sections, PADs, and β anisotropy. ezDyson evaluates the integral in Eq. (15) by numerical integration on a grid using Dyson or- bitals computed by a quantum chemistry program (Dyson orbitals are specified as an expansion over the atom-centered Gaussian basis functions). The treatment of the photoelectron wave function ψelk is one of the key challenges for quanti- tative modeling of photoionization and photodetachment. Currently, ezDyson uses two simple models—plane wave and Coulomb wave. The limitations of these models are well documented and can be addressed by more sophisticated approaches.^66 In the plane-wave representation, the photoelectron wave-function is
ψelk (r) = (^) (2π^1 ) 3 / 2 eik·r, (16)
where (2π)−^3 /^2 is the continuum normalization for plane waves.^67 ezDyson provides two imple- mentations for computing the respective matrix elements: one uses Eq. (16) as-is and another uses a partial-wave expansion in which the plane wave is expressed as a sum of spherical waves:^68
eik·r^ = 4π
∑^ ∞
l=
∑^ l m=−l
ilRl(kr)Ylm(r)Ylm(k). (17)
By using the partial-wave expansion, one can analyze the wave function of the ejected photoelectron in terms of angular momentum quantum numbers. Each spherical wave is char- acterized by its energy (E = 2 km^2 ) and angular momentum (l, m), and is a product of a radial function Rl(kr) and spherical harmonic functions Ylm. Ylm(r) depends on the position vector r and, therefore, determines the shape of the plane wave in position space, while Ylm(kˆ) de- termines the direction in which the plane wave travels.^68 By providing the coefficients of each partial wave, ezDyson enables quantitative application of the mixed s&p model of Sanov to polyatomic molecules.63, 69 In practice, the expansion in Eq. (17) is truncated at some lmax (specified by the user). In
