Evaluating the interactions between vibrational modes and electronic transitions using frontier orbital energy derivatives

Abstract

Vibrations affect molecular optoelectronic properties, even at zero kelvin. Accounting for these effects using computational modelling is costly, as it requires many calculations at geometries distorted from equilibrium. Here, we propose a low-cost method for identifying vibrations most strongly coupled to the electronic structure, based on using orbital energy derivatives as a diagnostic.

---

Predicting the barriers to electronic excitation and charge transfer is essential for engineering materials for solar cells,¹ mimicking redox processes in enzymes,² and building molecular electronic devices.³ Advances in electronic structure theory have enabled accurate descriptions of electronic states and the couplings between them,^4–7 facilitating in silico design and screening of potential next-generation molecular devices. However, most work on predicting molecular properties such as vertical excitation energies, polarizabilities, or single-molecule conductances is limited to equilibrium geometries, ignoring effects of nuclear motion.^4,5,8,9 This approach works well when the coupling between electronic states and nuclear motion is small and symmetric. In other words, if movement in direction +Q and −Q (which can be represented as a linear combination of normal modes Q_i) changes some property by a small amount +δ and −δ, respectively, the effect of nuclear motion can be approximated by a mere broadening (proportional to dδ/dQ) of the observed property.

Recently, Bai et al.¹⁰ demonstrated that this approximation does not always hold by showing that computed vertical excitation energies usually shift to lower energies when effects of nuclear motion are included. Alvertis and collaborators¹¹ found that this gap renormalization is present even at zero kelvin (zero-point renormalization, or ZPR), and estimated its mean value to −0.35 eV for Thiel's set⁴ of small organic molecules, with extreme cases reaching −1.36 eV. They also noted that ZPR can often be reduced to the effect of only several vibrations. A similar conclusion was reached by Lambropoulos et al., who found that ZPR in cyclo[18]carbon can be reduced to the effect of a few bond-stretching vibrational modes.¹² In organic molecular wires, bond stretching and backbone torsion play a large role, as the charge-carrier mobilities are highly dependent on the extent of conjugation modulated by these vibrations.^13–15 In more general cases, it can be challenging to understand which vibrations are important and which can be neglected.

Here, we provide a straightforward, inexpensive method for investigating the effect of individual normal modes on the electronic structure, based on calculating derivatives of frontier orbital energies with respect to molecular vibrations. We first use this method to evaluate how vertical excitation energies are modulated by specific vibrations, and then we analyse the effect of specific vibrations in a system with strong electron-vibration coupling.

In the framework of time-dependent density functional theory (TD-DFT), the first optical transition can usually be attributed to an excitation of a single electron from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO). The energy difference between the HOMO (E_H) and the LUMO (E_L) can be taken as a measure of the fundamental gap (E_gap), which is more rigorously defined as the difference between the ionization potential and the electron affinity.¹⁶ At the equilibrium geometry Q_eq = 0, the first vertical excitation energy (E_S1; also called the optical gap) will differ from the fundamental gap (E_gap) by the exciton binding energy (E_bind):

E_S1(Q_eq = 0) = E_gap(0) − E_bind(0) = E_L(0) − E_H(0) − E_bind(0) (1)

In contrast to E_gap, which is a one-electron quantity (difference between orbital energies), E_bind includes the two-body Coulombic attraction between the newly formed electron–hole pair. When the geometry is displaced along a normal mode i by some amount by Q_i, the optical gap changes by:

ΔE_S1(Q_i) = [E_L(Q_i) − ΔE_L(0)] − [E_H(Q_i) − ΔE_H(0)] − [E_bind(Q_i) − ΔE_bind(0)] (2)

which can also be written as:

ΔE_S1(Q_i) = ΔE_gap(Q_i) − ΔE_bind(Q_i) (3)

In practice, the effect of nuclear motion on the S₁ energy ΔE_S1 is determined by performing a series of excited-state calculations on geometries displaced by Q, which can either be obtained by Monte Carlo sampling (50–100 samples, in which Q for each geometry is a random linear combination of Q_i) or by finite differences (which usually involves displacements by +Q_i and −Q_i for each normal mode i). However, if the exciton binding energy is much less sensitive to geometric changes than frontier orbital energies (i.e., ΔE_bind ≪ ΔE_gap), then the changes in the optical gap due to nuclear motion ΔE_S1 can be approximated by the change in the fundamental gap ΔE_gap. We call this neglect of the geometry dependence of the two-body terms (ΔE_S1 ≈ ΔE_gap) the static binding approximation (SBA). SBA draws from early work on qualitative MO theory,¹⁷ and it is commonly used the popular Hubbard–Peierls model (see ESI†).¹⁸ SBA is attractive for two reasons. One, it simplifies the evaluation of effects of nuclear movement. When the harmonic approximation is valid (small Q), ΔE_gap can be calculated from derivatives of orbital energies with respect to normal modes, dE_H/dQ_i and dE_L/dQ_i, which can be obtained in a single ground-state calculation.¹⁹ Two, the SBA enables us to identify the vibrations most strongly coupled to the electronic structure from the magnitude of their orbital energy derivatives with respect to nuclear movement.

To investigate how well the static binding approximation holds, we tested it on the Thiel's set of small organic molecules (Fig. 1), for which the ZPR was investigated in ref. 11. We computed ΔE_gap using orbital energy derivatives and compared it to ΔE_S1 obtained using TD-DFT (B3LYP²⁰/def2-SVP), calculated for geometries modulated along each normal mode Q_i by the standard deviation of the thermal distribution σ_i(T):¹¹

where ω_i is the vibrational frequency, T is the temperature, and p_i(T) is the Bose–Einstein population at T, equal to unity at T = 0. Eqn (4) gives us a way to account for zero-point vibrations by looking at only two characteristic geometries modulated from the equilibrium by ±σ_i(0), as shown by Hele et al.¹¹

Fig. 1 Electronic-vibrational coupling calculated using ΔE_gap (vertical axis) and ΔE_S1 (horizontal axis; both in meV) for selected molecules from Thiel's set, decomposed by normal modes (triangles) and evaluated for zero-point vibrations. Black arrows in (m), (n), (q) and (r) correspond to stretching vibrations with a large ΔE_S1 relative to ΔE_gap. Results for remaining molecules are given in Fig. S1 (ESI†).

We find that SBA holds reasonably well across Thiel's set (Fig. 1 and Fig. S1–S3, ESI†): in 18 out of 30 molecules (green background), all normal modes with a notable ΔE_S1 can be diagnosed from ΔE_gap, with most of these points being close to the ΔE_gap = ΔE_S1 diagonal (usually R² > 0.8). In five further cases, most vibrations with a significant ΔE_S1 can be easily identified from ΔE_gap, but the importance of a single, usually stretching (∼1600 cm⁻¹; black arrow in Fig. 1m, n, q and r) vibration is consistently underestimated. SBA does not hold (Fig. 1r and Fig. S2f–h (ESI†); orange background) in very small molecules such as ethene or cyclopropane (Fig. S2g and h, ESI†) due to their tightly bound excitons (i.e., large two-body contributions). Finally, SBA also performs well in molecules where the first excitation does not correspond to a HOMO–LUMO transition (Fig. S3, ESI†).

Our results (Fig. 1 and Fig. S1–S3, ESI†) illustrate that in most small organic molecules only a few (usually < 15%) normal modes are strongly coupled to frontier orbitals. ΔE_gap can be used to identify the vibration most strongly coupled to frontier orbitals with high confidence (excluding very small systems), and in the majority of cases it correctly identifies vibrations strongly coupled with electronic structure, without resorting to excited-state calculations.

As E_bind decreases with molecular size,²¹ we expect the SBA to be even more accurate in molecules larger than those investigated here. Therefore, when a description of molecular properties beyond the equilibrium geometry is desired, we suggest using ΔE_gap to identify normal modes which are likely to affect the electronic structure. Limitations of the SBA are the following: (1) ΔE_gap only determines the size of the orbital coupling, and the calculation of its asymmetry requires excited state calculations. (2) The values of ΔE_gap are only valid close to the equilibrium geometry. (3) SBA will typically overestimate ΔE_gap in cases where frontier orbitals are (near) degenerate (e.g.Fig. 1o and p). (4) SBA can be used for transitions involving other orbitals, but this requires a single excited-state calculation (see ESI† Section 4 for details).

Accounting for effects of nuclear motion is particularly important in cases of strong electron-vibration coupling. Usually, in organic molecules the electronic transitions (typically v_el > 10 000 cm⁻¹) are well-separated from vibrational transitions (v_n ≤ 3200 cm⁻¹), with much broader and more intense signals. In radical cations of π-conjugated porphyrin oligomers this distinction is less clear, with both electronic and vibrational transitions producing very intense signals in the 1500–6000 cm⁻¹ region.^22,23 These unusually strong vibrational transitions can be detected using infrared (IR) spectroscopy, and they are usually called IR active vibrations (IRAVs). Recently, IRAVs were experimentally found at ∼1330, ∼1550 and ∼2080 cm⁻¹ (shaded area in Fig. 2d–f) in the butadiyne-linked porphyrin dimer 1˙⁺ (Fig. 2a), and tentatively assigned to normal modes shown in Fig. 2b. Movement along these IRAVs was attributed to the strong coupling of vibrational modes with the singly occupied molecular orbital (SOMO).²² Here, we investigate this coupling in more detail by analysing how normal modes in 1˙⁺ modulate orbital energies (Fig. 2c–f). Using the same density functional as in previous work (LC-ωPBE^22,24), we calculated the couplings between the HOMO and SOMO energies and normal modes in 1 and 1˙⁺ and compared them with computed dipole moment changes with respect to normal modes (i.e. infrared intensities) dD/dQ_i (Fig. 2c–d). Our results show that these frontier orbital energy derivatives are indeed larger in 1˙⁺ than in neutral 1 (Fig. 2c and d), in agreement with previous work. However, we also note that the contribution of orbitals up to SOMO–3 (Fig. 2e) is significant in IRAVs at 1390 cm⁻¹ and 2100 cm⁻¹, outweighing the contribution of the SOMO.

Fig. 2 (a) Structure of 1. (b) Calculated IRAVs of 1˙⁺. (c)–(d) Coupling between the computed IR intensities (positive y axis) and derivatives of frontier orbital energies (negative y axis) in (c) 1 and (d) 1˙⁺. Note the different y scales. (e) Contribution of lower-lying orbitals of 1˙⁺. (f) Modulation of the fundamental gap ΔE_gap/dQ_i in 1˙⁺ with molecular vibrations.

While the dipole moment D only depends on the occupied orbitals (within the DFT framework), a proper account of electron-vibration coupling should include the unoccupied orbitals as well. In 1˙⁺, we find that the singly unoccupied molecular orbital (SUMO) is also strongly coupled to movement along several vibrations, resulting in a very large ΔE_gap (Fig. 2f) at the IRAV positions. Therefore, a large ΔE_gap may be used as a diagnostic to identify IRAVs and vibrational modes strongly coupled to electronic states, despite being calculated at the equilibrium geometry.

An advantage of using ΔE_gap instead of the computed IR spectrum is that modes with no IR intensity (dD/dQ_i = 0) due to symmetry can also be IRAVs; indeed, some of the first IRAVs were observed in Raman-active modes.²⁵ In case of 1˙⁺, the symmetric stretch of the butadiyne linker at 2156 cm⁻¹ (see Fig. 1f and Fig. S4, ESI†) is IR forbidden, but it is strongly coupled to the ground electronic state as it decreases the conjugation between the two porphyrins.

In conclusion, we have presented a simple approach for investigating the influence of individual modes on molecular properties, based on calculating orbital energy derivatives with respect to molecular vibrations. By applying this approach to the Thiel's benchmark set of small organic molecules, we demonstrated that it can in most cases correctly identify the vibrations that are most strongly coupled to the vertical excitation energies by using only a ground-state calculation, which is significantly more efficient than using many excited-state calculations.

Molecules with unusually strong coupling between electronic and vibrational states show intense infrared active bands. Using a porphyrin dimer radical cation as an example, we showed how infrared spectrum can be decomposed into orbital contributions, finding a strong coupling between the IRAVs and the singly occupied molecular orbital (SOMO), but also a significant effect of orbitals up to SOMO–3. Finally, we showed that the sensitivity of frontier orbital energies to nuclear movement can be used to identify IRAVs. We also extended this approach to show how temperature affects the coupling between electronic states in a mixed-valence system.²⁶

This work was supported by the ERC (886506 ARO-MAT) and the UKRI (EP/X030075/1). Computational resources were provided by the Cirrus UK National Tier-2 HPC Service at EPCC (https://www.cirrus.ac.uk) funded by the University of Edinburgh and EPSRC (EP/P020267/1).

Data availability

Code for calculating orbital energy derivatives and building characteristic geometries using eqn (4) is provided at https://github.com/lisa-schroeder/mode-resolved-molecular-properties.

Conflicts of interest

There are no conflicts to declare.

References

1. R. E. Blankenship, D. M. Tiede, J. Barber, G. W. Brudvig, G. Fleming, M. Ghirardi, M. R. Gunner, W. Junge, D. M. Kramer, A. Melis, T. A. Moore, C. C. Moser, D. G. Nocera, A. J. Nozik, D. R. Ort, W. W. Parson, R. C. Prince and R. T. Sayre, Science, 2011, 332, 805–809.
2. C. C. Page, C. C. Moser, X. Chen and P. L. Dutton, Nature, 1999, 402, 47–52.
3. S. Ho Choi, B. Kim and C. D. Frisbie, Science, 2008, 320, 1482–1486.
4. M. Schreiber, M. R. Silva-Junior, S. P. A. Sauer and W. Thiel, J. Chem. Phys., 2008, 128, 134110.
5. P.-F. Loos, D. A. Matthews, F. Lipparini and D. Jacquemin, J. Chem. Phys., 2021, 154, 221103.
6. A. Kubas, F. Gajdos, A. Heck, H. Oberhofer, M. Elstner and J. Blumberger, Phys. Chem. Chem. Phys., 2015, 17, 14342–14354.
7. D. Manna, J. Blumberger, J. M. L. Martin and L. Kronik, Mol. Phys., 2018, 116, 2497–2505.
8. D. Toffoli, M. Quarin, G. Fronzoni and M. Stener, J. Phys. Chem. A, 2022, 126, 7137–7146.
9. J. O. Thomas, B. Limburg, J. K. Sowa, K. Willick, J. Baugh, G. A. D. Briggs, E. M. Gauger, H. L. Anderson and J. A. Mol, Nat. Commun., 2019, 10, 4628.
10. S. Bai, R. Mansour, L. Stojanović, J. M. Toldo and M. Barbatti, J. Mol. Model., 2020, 26, 107.
11. T. J. H. Hele, B. Monserrat and A. M. Alvertis, J. Chem. Phys., 2021, 154, 244109.
12. K. Lambropoulos, A. M. Alvertis, A. Morphis and C. Simserides, Phys. Chem. Chem. Phys., 2022, 24, 7779–7787.
13. D. Huang, H. Yao, Y. Cui, Y. Zou, F. Zhang, C. Wang, H. Shen, W. Jin, J. Zhu, Y. Diao, W. Xu, C.-A. Di and D. Zhu, J. Am. Chem. Soc., 2017, 139, 13013–13023.
14. J. Sukegawa, C. Schubert, X. Zhu, H. Tsuji, D. M. Guldi and E. Nakamura, Nat. Chem., 2014, 6, 899–905.
15. L. Matasović, B. Panić, M. Bubaš, H. Vančik, I. Biljan and I. Rončević, J. Mater. Chem. C, 2022, 10, 5433–5446.
16. E. J. Baerends, O. V. Gritsenko and R. van Meer, Phys. Chem. Chem. Phys., 2013, 15, 16408–16425.
17. K. Ruedenberg, J. Chem. Phys., 1977, 66, 375–376.
18. W. Barford, R. J. Bursill and M. Y. Lavrentiev, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 195108.
19. M. Bürkle, J. K. Viljas, T. J. Hellmuth, E. Scheer, F. Weigend, G. Schön and F. Pauly, Phys. Status Solidi B, 2013, 250, 2468–2480.
20. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
21. P. K. Nayak, Synth. Met., 2013, 174, 42–45.
22. W. J. Kendrick, M. Jirásek, M. D. Peeks, G. M. Greetham, I. V. Sazanovich, P. M. Donaldson, M. Towrie, A. W. Parker and H. L. Anderson, Chem. Sci., 2020, 11, 2112–2120.
23. D. P. Arnold, R. D. Hartnell, G. A. Heath, L. Newby and R. D. Webster, Chem. Commun., 2002, 754–755.
24. C. A. Jiménez-Hoyos, B. G. Janesko and G. E. Scuseria, Phys. Chem. Chem. Phys., 2008, 10, 6621–6629.
25. R. Österbacka, C. P. An, X. M. Jiang and Z. V. Vardeny, Science, 2000, 287, 839–842.
26. S. Kopp, A. Redman, I. Rončević, L. Schröder, L. Bogani, H. L. Anderson and C. R. Timmel, ChemRxiv, 2024, preprint DOI:10.26434/chemrxiv-2024-5kkhf.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cc02066a

---