Pitfall in simulations of vibronic TD-DFT spectra: diagnosis and assessment

Abstract

In this Communication, we study the effect of spurious oscillations in the profiles of energy derivatives with respect to nuclear coordinates calculated with density functional approximations (DFAs) for formaldehyde, pyridine, and furan in their ground and electronic excited states. These spurious oscillations, which can only be removed using extensive integration grids that increase enormously the CPU cost of DFA calculations, are significant in the case of third- and fourth-order energy derivatives of the ground and excited states computed by M06-2X and ωB97X functionals. The errors in question propagate to anharmonic vibronic spectra computed under the Franck–Condon approximation, i.e., positions and intensities of vibronic transitions are affected to a large extent (shifts as significant as hundreds of cm⁻¹ were observed). On the other hand, the LC-BLYP and CAM-B3LYP functionals show a much less pronounced effect due to spurious oscillations. Based on the results presented herein, we recommend either LC-BLYP or CAM-B3LYP with integration grids (250, 974) (or larger) for numerically stable simulations of vibronic spectra including anharmonic effects.

---

Electronic UV/Vis spectroscopy is a fundamental tool for studying the electronic structure of molecules. A proper understanding of the relations between the characteristics of molecular architectures and electronic structures is pivotal for the rational design of new species presenting desired photophysical properties.¹ In particular, fluorescent organic molecules are in the limelight due to widespread applications.^2,3 Thanks to the developments in instrumentation, it is now possible to use intense optical fields and excite molecules by multiphoton absorption techniques.^1,4 Syntheses and experimental photophysical studies are frequently paralleled by electronic-structure calculations for interpreting and rationalizing the complex spectroscopic signatures in electronic UV/Vis spectra of molecules. The current popularity of computational spectroscopy tools was preceded by a few decades of electronic/vibrational–structure theory developments^5–22 and extensive method validation.^23–42 Even though some of the developed methods present “chemical accuracy”, it is not yet feasible to apply them in studies of the electronic structure of molecular systems composed of dozens of atoms. The determination of vibrationally-resolved electronic absorption/emission spectra implies a further increase in the computational cost as it requires the calculation of the gradient, and depending on the vibrational structure method and also the hessian, of the electronic excited states.²⁸ Provided that medium-sized (or large) molecules are to be studied, the introduction of electrical or mechanical anharmonicity into the simulation protocol of their vibrationally-resolved electronic spectra narrows down the panel of electronic-structure theories to the Tamm–Dancoff approximation or random-phase approximation at the Hartree–Fock or Density Functional Theory (DFT) level. The latter approach has been demonstrated to exhibit superior accuracy compared to the Hartree–Fock theory, owing to its consideration of electron correlation effects. Nevertheless, in the case of modeling electronic spectra, numerous density functional approximations (DFAs) exhibit serious flaws. The unsatisfactory performance of DFAs in the description of electronic spectra has already been well understood and improved. However, currently, the selection of the best DFA for simulations of vibrationally-resolved electronic spectra is still by no means a routine task.⁴³ Considering that simulations of harmonic-oscillator-based vibrationally-resolved electronic spectra for medium-sized molecules using DFAs have a low computational cost,⁴⁴ one may anticipate further improvements in vibrational-structure theories through the incorporation of mechanical and electrical anharmonicity. The calculation of such anharmonic corrections requires the determination of high-order energy/property derivatives.

The efficiency and robustness of DFT pave the way to applying vibrational-structure models beyond harmonic approximation. However, in our recent study we have demonstrated that the majority of developed DFAs suffer from “spurious oscillations” in the ground-state profiles of energy/property derivatives with respect to nuclear coordinates.⁴⁵ This unphysical effect was traced back to the numerical integration required to calculate the exchange–correlation energy and can be fixed using extensive integration grids that increase enormously the CPU cost of DFA calculations. Although the first pieces of evidence of the spurious oscillations were found in the profiles of energy/property derivatives with respect to intermolecular stretching modes of molecular complexes,^45,46 they are also present in the vibrations of isolated molecules. Based on the results obtained for the ground-state properties,⁴⁵ it can be assumed that vibronic spectra might also be affected by the spurious oscillations.

The aim of this Communication is to perform a pioneering exploration of the effect of spurious oscillations on the vibrational fine structure of absorption bands in electronic absorption spectra. For the first time we study, on equal footing, the effect of the spurious oscillations on ground- and excited-state energy derivatives as well as transition property derivatives. With computational efficiency in mind, for the present study, we selected three molecules (formaldehyde, pyridine, and furan) and four DFAs (LC-BLYP,⁴⁷ CAM-B3LYP,⁴⁸ωB97X,⁴⁹ and M06-2X⁵⁰). In our previous study, LC-BLYP and CAM-B3LYP were two of the DFAs showing more robustness against spurious oscillations among the 45 functionals studied. In contrast, ωB97X and M06-2X were two of the most strongly affected by spurious oscillations.⁴⁵ Thus, these four DFAs cover both extremes of the spectrum regarding the potential magnitude of the spurious oscillations in energy/property derivatives and their propagation to vibronic spectra.

In this Communication, we analyze the grid-dependent spurious oscillations in energy derivative profiles using the procedure and in-house codes proposed in our recent study.⁴⁵ Namely, for a given DFA, we compute ground- and excited-state energies and their derivatives along the vibrational normal mode displacements (denoted as Q), hereafter labeled as P^DFA. This is done for (99, 590), (250, 974), and (750, 974) unpruned integration grids, and the results obtained with the largest grid are used as the reference for a given DFA, (a,b) referring to a grid with a and b radial and angular points, respectively. An in-house algorithm⁴⁵ detects and filters out the possible grid-related oscillations in the PES obtained with the largest grid (the reference grid). This is done in the reciprocal (frequency) space, i.e., the PES is transformed using the Discrete Fourier Transform into the so-called frequency spectrum. If present, these spurious oscillations are identified by comparison with the frequency spectrum of spurious-oscillation-free PES obtained with an ab initio method (HF and CIS for the ground and excited states, respectively, in the case of the present Communication). Subsequently, the identified spurious oscillatory bands are quantitatively filtered out using an automatically designed low-pass finite impulse response filter.^45,51 Such a filtered and spurious-oscillation-free PES is used as the reference grid, P^{DFA(large grid)}_filt, to quantify the spurious oscillations in the PES of the target DFA and integration grid, P^DFA. A thorough description of the algorithm can be found in the ESI† of our previous work.⁴⁵

We will quantify the errors due to spurious oscillations by defining the relative root-mean-square errors (RRMSE):

The RRMSEs were computed for the third and fourth energy derivatives with respect to vibrational normal modes for three molecules in their ground and electronic excited states. During RRMSE evaluations, only atomic displacements close to equilibrium geometry (for the selected electronic state) were included, and they corresponded to the displacements in the normal mode coordinate Q. In particular, displacements used in the RRMSE evaluations corresponded to ± 0.16|Q| (with a step size of 0.01|Q|) for formaldehyde, where |Q| is a normalized normal mode coordinate (calculated in atomic units), whereas for furan and pyridine, these displacements corresponded to ± 0.32|Q| (with a step size of 0.02|Q|).

In what follows, we will present the results of calculations obtained using the GAUSSIAN package⁵² and 6-311++G** atomic basis set⁵³ to compute ground and excited state energies/properties.

The results of the analysis for fourth energy derivatives are shown in Table 1, and the corresponding data for third energy derivatives can be found in Table S1 (Tables S2 and S3 contain the complementary results obtained for the second and third derivatives of transition dipole moments, ESI†). Note that the RRMSE values given for each molecule correspond to the average RRMSE for all the vibrational normal modes. Four key conclusions can be drawn from the data presented in Table 1. First, for all studied electronic states of the three molecules, the RRMSE values are significantly reduced on passing from (99, 590) to (250, 974) grid. Further extension of the number of points in the radial part of the grid does not lead to a noticeable decrease in the RRMSE. Second, RRMSE values of fourth energy derivatives in electronic excited states might be comparable in magnitude to the ground state counterparts. Indeed, for the formaldehyde studied using M06-2X functional, the excited states present larger RRMSE than the ground states. Third, CAM-B3LYP and LC-BLYP exhibit pleasingly low values of RRMSE for ground and excited states, which parallels the conclusion from our earlier work⁴⁵ devoted to ground state energy/property derivatives for other chemical systems (i.e., average RRMSE for all the studied molecules and electronic states smaller than 11% for (99, 590) grid and smaller than 2.5% for (250, 974) grid). Fourth, the magnitude of RRMSE values is highly system-dependent which hints at differences in vibrational structure (i.e., different vibrational normal mode types). The same overall general pattern of RRMSE values holds for third energy derivatives with respect to vibrational normal modes (see Table S1, ESI†).

Table 1 Average of the relative root mean square error (RRMSE) of d⁴E/dQ⁴ for all normal modes of different electronic states of formaldehyde, pyridine, and furan, computed with tested DFAs and three unpruned integration grids. For all molecular systems, 1¹A₁ denotes the electronic ground–state, and other states are selected bright excited states. In the case of furan, 1²A₂ corresponds to the lowest electronic state of its ionized form, i.e., the electronic ground–state of the furan cation

		Formaldehyde			Pyridine		Furan
		1^1A1	1^1A2	1^1B2	1^1A1	1^1B1	1^1A1	1^2A2
LC-BLYP	(99, 590)	1	1	2	20	1	48	1
	(250, 974)	0	0	2	3	0	6	0
	(750, 974)	0	0	2	3	0	6	0
CAM-B3LYP	(99, 590)	1	1	3	21	3	18	1
	(250, 974)	0	0	1	6	0	8	0
	(750, 974)	0	0	1	6	0	8	0
ωB97X	(99, 590)	11	22	16	691	460	521	72
	(250, 974)	1	6	2	274	193	332	39
	(750, 974)	1	6	1	264	153	332	39
M06-2X	(99, 590)	47	59	67	556	117	342	39
	(250, 974)	20	48	40	158	62	323	16
	(750, 974)	20	48	39	157	61	322	15

---

Now we will inspect how the RRMSE errors propagate to the vibronic spectra. To that end, we have selected M06-2X, ωB97X, and LC-BLYP functionals, the (75, 302), (99, 590), (250, 974), and (750, 974) unpruned grids, and simulated the 1¹A₁ → 1¹B₂ (S₀→S₂) electronic transition for formaldehyde under Franck–Condon (FC) approximation. Note that the UltraFine grid, a pruned (99, 590) grid that contains a few times smaller number of points than its unpruned parent, is the default grid used in the GAUSSIAN 16 program. For simplicity, we have focused our analysis on the 2¹, 4¹, and 5¹ vibronic transitions. In these simulations, we accounted for mechanical anharmonicity by adopting the perturbation theory-based method of Luis et al.⁵⁴ Note that these calculations require up to fourth energy derivatives of the electronic ground and excited states. Therefore, they are excellent tests for analyzing the effect of spurious oscillations on vibrational spectroscopic signatures. The results of these simulations are shown in Fig. 1 (M06-2X and LC-BLYP) and in Fig. S2 (ESI†) (ωB97X).

Fig. 1 Normalized PT2 FC spectra of the 1¹A₁ → 1¹B₂ (S₀→S₂) transition of formaldehyde obtained with M06–2X (top) and LC-BLYP (bottom) using various integration grids and 0.02|Q| step size for numerical differentiation. On the horizontal axis, differences with respect to the anharmonic 0–0 band are shown (ΔE_0n = E_0n − E₀₀). On the legend, next to the grid labels, the corresponding errors in the anharmonic transition energies of the 0–0 band, when compared to the reference (750, 974) grid, are given (E^err₀₀ = E₀₀ − E^(750,974)₀₀).

In the case of the M06-2X functional, the differences in the positions of the three selected FC transitions between the reference (750, 974) grid and the smaller grids may be larger than 500 cm⁻¹ (i.e., 5¹ frequency obtained with (75, 302) and (99, 590)). Indeed, for this DFA, even using the (250, 974) grid the shifts in the vibronic lines are still very large, thus making band assignments in experimental spectra of molecules difficult. Moreover, the transition intensities are also affected, albeit to a smaller extent. The overall effect of spurious oscillations has a significant impact on the spectroscopic signatures in the vibronic spectra of formaldehyde, as anticipated based on the RRMSE data presented in Table 1. The results obtained using the ωB97X functional (Fig. S2, ESI†) qualitatively parallel those already discussed for the M06-2X functional. However, the performance of the ωB97X functional is much more satisfactory when using smaller integration grids, which agrees with the analysis of the RRMSE results for formaldehyde shown in Table 1. Let us finally note that small RRMSE values found for LC-BLYP functional translate into very small differences in vibronic FC spectra between the smaller grids and the reference grid (Fig. 1).

In summary, we have studied the effect of spurious oscillations in the profiles of energy derivatives with respect to nuclear coordinates for formaldehyde, pyridine, and furan in their ground and electronic excited states. A significant effect of spurious oscillations, especially for (75, 302) and (99, 590) unpruned grids, was observed in the case of third- and fourth-order energy derivatives computed by M06-2X and ωB97X functionals (and third-order derivatives of transition moment, a quantity relevant for electrical anharmonicity effects). The errors in question propagate to anharmonic vibronic spectra computed under the Franck–Condon approximation, i.e., positions and intensities of vibronic transitions are affected to a large extent (shifts as significant as 500 cm⁻¹ were observed). On the other hand, the LC-BLYP and CAM-B3LYP functionals show much smaller errors due to spurious oscillations, which agrees with the extensive calculations for many atomic and molecular systems in their ground state.⁴⁵ Based on the results presented herein, we recommend either LC-BLYP or CAM-B3LYP with integration grids (250, 974) (or larger) for numerically stable simulations of vibronic spectra including anharmonic effects. One should not overlook, however, that range-separated functionals might not always be the best choice for simulations of vibrationally resolved spectra.⁴³

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

R. Z. and S. P. S. thank the National Science Centre (Poland) for financial support (grant no. 2018/30/E/ST4/00457) and the Wroclaw Center for Networking and Supercomputing (Poland) for computational resources. J. M. L. and E. M. acknowledge funding from Agencia Española de Investigación, “FEDER Una manera de hacer Europa” (PID2022-140666NB-C21 and PID2022-140666NB-C22), the Donostia International Physics Center (DIPC-INV-003132 and DIPC-INV-003133), and the Generalitat de Catalunya (JML) (2021SGR00623).

Notes and references

1. M. Pawlicki, H. A. Collins, R. G. Denning and H. L. Anderson, Angew. Chem., Int. Ed., 2009, 48, 3244–3266.
2. D. Frath, J. Massue, G. Ulrich and R. Ziessel, Angew. Chem., Int. Ed., 2014, 53, 2290–2310.
3. Q. Zheng, G. Xu and P. Prasad, Chem. – Eur. J., 2008, 14, 5812–5819.
4. M. Drobizhev, N. S. Makarov, S. E. Tillo, T. E. Hughes and A. Rebane, Nat. Methods, 2011, 8, 393–399.
5. R. M. Dreizler and E. K. U. Gross, Density functional theory, Springer-Verlag, Berlin Heidelberg, 1990.
6. D. R. E. Yarkony, Modern electronic structure theory, World Scientific Publishing Company, 1995.
7. M. E. Casida, in Recent advances in density functional methods, ed. D. P. Chong, World Scientific, Singapore, 1995, part I, pp. 155–192.
8. A. Rubio and M. Marques, Phys. Chem. Chem. Phys., 2009, 11, 4436.
9. Computational strategies for spectroscopy, ed. V. E. Barone, John Wiley & Sons, Ltd, 2011.
10. R. J. Bartlett and M. Musiał, Rev. Mod. Phys., 2007, 79, 291–352.
11. J. Olsen and P. Jørgensen, J. Chem. Phys., 1985, 82, 3235–3264.
12. F. Plasser, M. Wormit and A. Dreuw, J. Chem. Phys., 2014, 141, 024106.
13. S. M. Parker, D. Rappoport and F. Furche, J. Chem. Theory Comput., 2018, 14, 807–819.
14. A. Dreuw and M. Head-Gordon, Chem. Rev., 2005, 105, 4009.
15. C. Hättig, Adv. Quantum Chem., 2005, 50, 37–60.
16. G. Fischer, Vibronic coupling: The interaction between the electronic and nuclear motions, Academic Press, 1984.
17. J. Cerezo, J. Zuniga, A. Requena, F. J. Avila Ferrer and F. Santoro, J. Chem. Theory Comput., 2013, 9, 4947–4958.
18. D. W. Silverstein and L. Jensen, J. Chem. Phys., 2012, 136, 064111.
19. P. Macak, Y. Luo and H. Ågren, Chem. Phys. Lett., 2000, 330, 447–456.
20. M. Cossi and V. Barone, J. Chem. Phys., 2001, 115, 4708–4717.
21. R. Improta, G. Scalmani, M. J. Frisch and V. Barone, J. Chem. Phys., 2007, 127, 074504.
22. F. Kossoski, Y. Damour and P.-F. Loos, J. Phys. Chem. Lett., 2022, 13, 4342–4349.
23. L. Goerigk and S. Grimme, J. Phys. Chem. A, 2009, 113, 767–776.
24. M. Casanova-Páez and L. Goerigk, J. Chem. Phys., 2020, 153, 064106.
25. L. Goerigk, J. Moellmann and S. Grimme, Phys. Chem. Chem. Phys., 2009, 11, 4611–4620.
26. F. J. Avila Ferrer, J. Cerezo, E. Stendardo, R. Improta and F. Santoro, J. Chem. Theory Comput., 2013, 9, 2072–2082.
27. A. Petrone, J. Cerezo, F. J. A. Ferrer, G. Donati, R. Improta, N. Rega and F. Santoro, J. Phys. Chem. A, 2015, 119, 5426–5438.
28. F. J. A. Ferrer and F. Santoro, Phys. Chem. Chem. Phys., 2012, 14, 13549–13563.
29. F. J. A. Ferrer, R. Improta, F. Santoro and V. Barone, Phys. Chem. Chem. Phys., 2011, 13, 17007–17012.
30. R. Improta, V. Barone and F. Santoro, Angew. Chem., Int. Ed., 2007, 119, 409–412.
31. T. J. Zuehlsdorff, A. Montoya-Castillo, J. A. Napoli, T. E. Markland and C. M. Isborn, J. Chem. Phys., 2019, 151, 074111.
32. M. Dierksen and S. Grimme, J. Phys. Chem. A, 2004, 108, 10225–10237.
33. M. Dierksen and S. Grimme, J. Chem. Phys., 2004, 120, 3544–3554.
34. J. P. Götze, B. Karasulu and W. Thiel, J. Chem. Phys., 2013, 139, 234108.
35. C. Suellen, R. G. Freitas, P.-F. Loos and D. Jacquemin, J. Chem. Theory Comput., 2019, 15, 4581–4590.
36. R. Sarkar, P.-F. Loos, M. Boggio-Pasqua and D. Jacquemin, J. Chem. Theory Comput., 2022, 18, 2418–2436.
37. B. L. Guennic and D. Jacquemin, Acc. Chem. Res., 2015, 48, 530–537.
38. S. Chibani, A. D. Laurent, B. Le Guennic and D. Jacquemin, J. Chem. Theory Comput., 2014, 10, 4574–4582.
39. A. V. Marenich, C. J. Cramer, D. G. Truhlar, C. A. Guido, B. Mennucci, G. Scalmani and M. J. Frisch, Chem. Sci., 2011, 2, 2143–2161.
40. C. A. Guido, S. Knecht, J. Kongsted and B. Mennucci, J. Chem. Theory Comput., 2013, 9, 2209–2220.
41. J. Bednarska, R. Zaleśny, W. Bartkowiak, B. Ośmiałowski, M. Medved’ and D. Jacquemin, J. Chem. Theory Comput., 2017, 13, 4347–4356.
42. M. J. G. Peach, P. Benfield, T. Helgaker and D. J. Tozer, J. Chem. Phys., 2008, 128, 044118.
43. A. Charaf-Eddin, A. Planchat, B. Mennucci, C. Adamo and D. Jacquemin, J. Chem. Theory Comput., 2013, 9, 2749–2760.
44. T. Petrenko and F. Neese, J. Chem. Phys., 2012, 137, 234107.
45. S. P. Sitkiewicz, R. Zaleśny, E. Ramos-Cordoba, J. Luis and E. Matito, J. Phys. Chem. Lett., 2022, 13, 5963–5968.
46. R. Zalesny, M. Medved’, S. P. Sitkiewicz, E. Matito and J. M. Luis, J. Chem. Theory Comput., 2019, 15, 3570–3579.
47. H. Iikura, T. Tsuneda, T. Yanai and K. Hirao, J. Chem. Phys., 2001, 115, 3540–3544.
48. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57.
49. J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2008, 128, 084106.
50. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
51. J. McClellan and T. Parks, IEEE Trans. Circuit Theory, 1973, 20, 697–701.
52. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 Revision C.01, Gaussian Inc., Wallingford CT, 2016.
53. R. Ditchfield, W. J. Hehre and J. A. Pople, J. Chem. Phys., 1971, 54, 724–728.
54. J. M. Luis, D. M. Bishop and B. Kirtman, J. Chem. Phys., 2004, 120, 813–822.

---

Footnote

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

---