Systematic improvement of molecular excited state calculations by inclusion of nuclear quantum motion: A mode-resolved picture and the effect of molecular size

The optoelectronic properties of organic molecules are dominated by their low-lying excited states known as excitons.¹ Exciton energies are critical to several technologically relevant processes in these systems, such as singlet fission and thermally activated delayed fluorescence, which find applications in photovoltaics and LEDs, respectively.^2–5 It is therefore desirable to develop methods for the accurate prediction of exciton energies.

Typical methods for the calculation of excited state energies, such as time-dependent density functional theory (TD-DFT),⁶ coupled cluster (CC),⁷ complete active space self-consistent field (CAS-SCF),⁸ and the Bethe–Salpeter equation (BSE),^9–11 only account for electronic contributions to exciton states, typically computing vertical excitation energies at a fixed geometry of the system. However, in a recent study, Bai et al.¹² showed that the vibrational motion of molecules is responsible for a redshift of the absorption maximum compared to “static” vertical excitation energies, an effect that needs to be accounted for in order to achieve predictive accuracy. Even at 0 K, atomic nuclei vibrate with a zero-point energy of $12ℏω$ per normal mode. In a recent study on solid state organic semiconductors,¹³ two of us showed that this nuclear quantum motion can significantly change the “static” exciton energies that are commonly computed at the ground state geometry of a system and that incorporating these effects leads to improved agreement with experiment.

A number of computational studies have proposed advanced methods for accurately simulating the shape of molecular absorption spectra including vibrational effects.^14–16 Additionally, the effect of nuclear quantum fluctuations on molecular excited states is now understood to cause a redshift of exciton energies compared to their “static” values.¹² Despite these advances, there remain a number of open questions regarding the effect of quantum fluctuations on molecular excited states: (i) It remains unclear to what extent the inclusion of quantum fluctuations leads to a systematic improvement of molecular excited state calculations in terms of agreement with experiment for the exciton energy, which corresponds to the frequency of the absorption maximum in the case of a state with non-zero oscillator strength. (ii) How many vibrational normal modes significantly contribute to the renormalization of exciton energies due to quantum motion is, to the best of our knowledge, not yet understood, since the studies that have emerged so far do not provide a mode-resolved picture of this effect. (iii) Our previous work in periodic structures¹³ suggests that the change in the exciton energy due to zero-point motion depends on the size of the studied system. It is hence important to investigate whether such a trend also holds for isolated organic molecules. (iv) It is thus far not clear whether the level of theory employed for the calculation of the zero-point renormalization (ZPR) of the exciton energy has a large impact on its value or whether the magnitude of this effect is largely independent of the underlying level of calculations. This last point is particularly important for practical calculations, as such an independence would imply that one could compute the correction due to nuclear quantum effects at a cheaper level of theory than the “static” exciton energy, hence reducing the overall computational cost of the calculation.

In this work, we systematically investigate the effect of nuclear quantum motion on exciton energies of organic molecules and address the aforementioned questions. In order to accurately compute the zero-point renormalization of exciton energies, we employ a Monte Carlo sampling technique, combining TD-DFT and finite difference methods for the molecular vibrations.^17,18 This Monte Carlo method had thus far been used in the context of periodic systems¹³ and allows us to capture exciton–vibration interactions to all orders and to treat excited state energy surfaces without any harmonic assumption. It also has strong similarities to the nuclear ensemble approach,¹⁹ which has been developed for molecular systems and is also used in Ref. 12.

While both the nuclear ensemble and Monte Carlo methods capture the renormalization of exciton energies due to quantum fluctuations, they do not provide information on the individual contribution of normal modes to this effect. For example, low- and high-frequency vibrations in organic systems are known to play different roles during vibrational relaxation²⁰ and electron-transfer reactions;²¹ therefore, obtaining a mode-resolved picture for the exciton energy renormalization is critical to obtaining microscopic insights into the physics of these systems. Such mode-resolved information is available through the use of a quadratic approximation to the exciton–vibration coupling;¹⁸ however, this comes at the cost of only capturing these interactions to third order. Here, we employ the quadratic approximation for capturing the correction to exciton energies induced by nuclear quantum motion, and we assess the accuracy of this method by comparing to the more accurate Monte Carlo calculations.

We apply our set of methods to the so-called Thiel set of organic molecules (sometimes also referred to as the Mülheim set), for which highly accurate exciton energies have been computed using wavefunction-based methods.²² The Thiel set of molecules consists of four broad categories of structures: unsaturated aliphatic hydrocarbons, aromatic hydrocarbons and heterocycles, carbonyl compounds, and nucleobases. The studied structures are shown in Fig. 1. In addition to the molecules that are commonly included in the Thiel set, we also study three additional molecules, anthracene, tetracene, and pentacene, which together with benzene and naphthalene that are included in the Thiel set form the first five members of the acene family. We have also excluded propanamide and pyridazine from the studied structures, as we were unable to obtain converged vibrational properties for these molecules. It has become common in the literature to compare computational methods for the calculation of exciton energies to the Thiel set,^{10–12,23–25} a path that we also take. Moreover, throughout this paper, we systematically compare our results to experiment and find that both the Monte Carlo and quadratic methods for computing the renormalization of the exciton energies due to zero-point motion of the nuclei are remarkably accurate. This allows us to draw several important conclusions on the microscopic mechanism of this phenomenon and to assess different methods of describing it accurately.

The structure of this paper is as follows: In Sec. II, we provide the necessary theoretical background for the results that follow in Sec. III. Specifically, in Subsection II A, we include a qualitative discussion of the effect of molecular vibrations on exciton energies, and based on standard expressions, we lay the foundations to discuss the Monte Carlo sampling technique and the quadratic approximation for quantifying this effect in Subsections II B and II C, respectively. In Sec. III, we compare our results from a Monte Carlo sampling of the exciton energies to experiment and to previous computational studies in Subsection III A. The effect of using different levels of electronic structure theory on the predicted magnitude of the exciton energy renormalization is examined in Subsection III B, and the impact of the molecular size on the exciton energy correction is discussed in Subsection III C. We then consider the accuracy and speed of the quadratic method in Subsection III D and proceed to use it in order to attain a mode-resolved picture of the effect of nuclear quantum motion on exciton energies in Subsection III E. Finally, we summarize our results and conclude our study in Sec. IV.

We start by presenting standard results and a qualitative discussion of the effects of molecular vibrations on exciton energies, which help inform our later discussion. In organic molecules, the highest occupied molecular orbital (HOMO) is a bonding σ or π orbital that lowers the energy of the molecule once occupied and leads to nuclei that are closer to each other. Therefore, the energy of such an orbital along a generalized nuclear coordinate u of the molecular system will look similar to a Morse potential, as visualized in Fig. 2(a) (black). The lowest unoccupied molecular orbital (LUMO) is an anti-bonding σ* or π* orbital that raises the energy of the molecule once occupied and forces the nuclei apart, its energy along u showing an exponential decay [Fig. 2(a), red]. For most conjugated organic molecules, the lowest energy singlet exciton (S₁) is formed by exciting an electron from the HOMO to the LUMO,²⁶ with energy

Here, E(S₀) is the energy of the ground state, ϵ_H and ϵ_L are the energies of the HOMO and the LUMO, respectively, J_HL is the HOMO–LUMO Coulomb integral, and K_HL is the HOMO–LUMO exchange integral. By using Eq. (1), we schematically describe the energy of an exciton along u [Fig. 2(a), blue], having assumed that the dependence of the integrals J, K on u does not alter its qualitative characteristics.

What can be observed from Fig. 2(a) is that the excited state potential energy curve (blue) has a lower curvature compared to the ground state one (black) due to the contribution of the anti-bonding orbital. Therefore, if we approximate the two within the harmonic approximation, we obtain two parabolas of different curvature, as shown in Fig. 2(b). If the structure was “frozen” at its ground state configuration u_GS, then the energy required to access the exciton state would be E_static, which is the vertical excitation energy usually computed by electronic structure calculations. However, due to the vibrational motion of the nuclei, the system can explore a distribution of configurations within the region denoted by green arrows, even at 0 K. If we excited the system from any given displaced configuration within the available region, then the excitation energy becomes smaller than E_static, since the lower parabola corresponding to the ground state has a higher curvature than the excited state one. Therefore, it is expected that inclusion of the effects of nuclear fluctuations will lead to a redshift of exciton energies. In particular, at temperature T, the absorption maximum will not be found at the energy E_static corresponding to vertical excitation but at a lower energy corresponding to the quantum mechanical vibrational average,

where $|χsu〉$ is a vibrational eigenstate on the ground state potential energy surface with energy E_s, $Z=∑se−Es/kBT$ is the partition function, and u is the nuclear displacement.

An intuitive way of representing this mean exciton energy on a diagram of the potential energy surfaces of a molecule is to plot the transition energy at a “mean value” configuration u^MV, where the vertical exciton energy is equal to the average at temperature T,²⁷ 

According to the mean-value theorem for integrals, there always exists such a configuration u^MV. We can therefore visualize the mean exciton energy, corresponding to the absorption band maximum in Fig. 3. It is also worth noting that the absorption maximum described by the quantum mechanical expectation value of Eq. (2) is distinct from the so-called 0–0 energy that is often reported in the literature, which corresponds to the energy difference between the zero-point levels of the optimized ground and excited state configurations, as visualized in Fig. 3.

To simplify the problem of computing the exciton band maximum energy from Eq. (2), we use the harmonic approximation for the ground state potential, and by substituting $|χsu〉$ with the wavefunction of a quantum harmonic oscillator, we obtain¹⁸ 

where E_exc(u) is the vertical excitation energy at the configuration u and

is the harmonic density at temperature T, which, in turn, is a product of Gaussian functions of width

In the above equation, atomic units and mass-weighted coordinates have been used, and ν is the index labeling the ground state vibrational modes of the studied molecule, with ω_Gν being the corresponding frequency. For non-linear molecular systems, there are six vibrational modes with zero frequency, which correspond to translations and rotations of the entire structure. These are not included in the sampling of the integral of Eq. (4).

The Monte Carlo method approximates the integral in Eq. (4) by generating configurations u that are distributed according to the harmonic density at the temperature of choice T. Then, the expectation value of Eq. (4) is computed as the simple average of the computed values of the vertical exciton energy E_exc(u) at the displaced configurations. This methodology for capturing the effect of vibrations on exciton energies has previously been applied in the context of periodic organic crystals,¹³ and here, we extend it to isolated organic molecules. This computational approach has the advantage that it makes no assumption about the shape of the excited state potential energy surface, and it also includes exciton–vibration interactions to all orders. Moreover, it relies on no adjustable parameters, apart from those in the DFT functional, which are fixed throughout the entire series of calculations. This Monte Carlo approach is conceptually very similar to the nuclear ensemble method,¹⁹ which was recently applied to the Thiel set of molecules,¹² and we compare our results to this previous study later on, in Fig. 5. For all the molecules studied, we generate 100 displaced configurations at T = 0 K and T = 300 K, and their exciton energies are computed through TD-DFT within the Tamm–Dancoff approximation²⁸ using the popular B3LYP hybrid functional²⁹ and the cc-pVDZ basis set, as implemented within the NWChem code.³⁰ We hence obtain renormalized exciton energies compared to those of a “static” TD-DFT calculation for the Thiel set of molecules. We generally find that ∼50 points are sufficient to sample the expectation value of the exciton energies (supplementary material, Sec. S4). The numerical results of these calculations (and the associated statistical uncertainties) are given in Sec. S1of the supplementary material, not only at the B3LYP level but also for two additional DFT functionals, as discussed in Subsection III B.

Since the Monte Carlo method provides an approximation to the expectation value of Eq. (2), we always compare the computed energy values to the absorption maximum, which corresponds to the same exciton state (for example, the first bright excited state) and not to the so-called 0–0 energy (see also Fig. 3). However, it is not uncommon that experimental works report the energy of the 0–0 transition. For the vast majority of systems studied here and unlike the schematic of Fig. 3 that aims to emphasize the differences between these two quantities, the energies of the 0–0 transition and the band maximum are either identical or very close to each other. Accurate explicit calculations of 0–0 energies have been reported in the literature,³¹ and these require finding the optimized geometry of the excited state and computing its Hessian matrix, which can be very challenging computationally, especially for larger molecules such as the acene series and nucleobases, which are studied here. Given that the energy of the 0–0 transition is in most cases very close or even identical to that of the band maximum, the Monte Carlo method (and also the quadratic approximation that is outlined in Subsection II C) could potentially be used to also approximate 0–0 energies at a much lower computational cost.

One can further simplify the expression of Eq. (4) by performing a quadratic expansion of the vertical excitation energy E_exc(u) in the coordinates u_ν of the vibrational normal modes, yielding

Here, E_exc(0) is the exciton energy at the ground state geometry and is equal to the static exciton energy E_static as defined previously. Substituting this expression in Eq. (4) gives

In Eq. (8), n_B(ω_ν, T) is the Bose–Einstein distribution for the vibrational quanta at temperature T. When computing the expectation value of the exciton energy in Eq. (8), all odd terms vanish due to the harmonic density |Φ(u; T)|² being an even function; thus, the resulting approximation is accurate to fourth order in u, and the exciton energy renormalization is described to lowest order by the quadratic term of the expansion. While this so-called quadratic approximation is, in principle, less accurate than the Monte Carlo approach, which we outlined previously, it has the advantage of separating the contributions of the different vibrational normal modes to the exciton energy renormalization, potentially allowing for additional microscopic insights into these effects.

Equation (8) is used for calculations with the quadratic method, which essentially involves the calculation of the second derivative of the exciton energy along each vibrational mode of the system by using the finite difference formula

Hence, two exciton energies at ±δu_ν for each mode ν, as well as the “static” exciton energy E_static need to be computed. For a molecule with N atoms, this means that a total of 2 · (3N − 6) + 1 exciton energy calculations are required, which for small molecules can be significantly smaller than the number of configurations required to converge a Monte Carlo calculation, a point that we return to in Subsection III D. In principle, δu is an infinitesimal quantity that should be as small as possible in the above finite difference formula; however, in practice, one needs to choose a finite value for this parameter in order to avoid numerical divergence issues. For perfectly quadratic energy surfaces, the result would be independent of the displacement δu, and we can therefore conceptually extend the meaning of the quadratic formula to simply represent a finite displacement of the system along an individual vibrational mode. We return to this issue in Subsection III D.

Finally, it is worth transforming Eq. (8) into a different form, which is not used for practical calculations within this work but provides valuable intuition for the effects that we discuss. Let us assume the frequency of vibrational normal mode ν in the ground state to be equal to $ωGν2$⁠. If we further assume that this same normal mode is present in the excited state with a different frequency $ωEν2$⁠, then Eq. (8) becomes

This formula is identical to the one appearing in Ref. 32, where E(T) is called the first moment of the exciton. In the terms we have used here, the first moment is simply the expectation value of the exciton energy in the presence of molecular vibrations. Equation (10) demonstrates that the renormalization of the static exciton energy due to each normal mode depends on the frequency difference $ωGν2−ωEν2$ and suggests that E(T) < E_static as long as the potential energy surface of the ground state has a higher curvature than that of the excited state and therefore $ωGν2>ωEν2$⁠. This analytical result is essentially the mathematical manifestation of the intuitive picture of Fig. 2 and holds in the general multi-dimensional case, at least within the validity of the quadratic approximation.

We compare the energies of the singlet excitons obtained within our Monte Carlo approach by computing exciton energies at the B3LYP level of TD-DFT and including the effects of zero-point renormalization (ZPR) to the experimental values for the corresponding maxima. The references to the experimental works are summarized in Sec. S2 of the supplementary material. These same references were collected in the original publication on the Thiel set;²² however, the values can be slightly different for a few of the studied molecules, as we had to ensure that we compare to the band maximum in every case. We also compare the values obtained from “static” TD-DFT and from highly accurate complete active space second-order perturbation theory (CASPT2)³³ calculations in the original Thiel publication²² to experiment. The computed exciton energies always refer to the first single excitation as described within TD-DFT, and we exclude any double-excitations (bi-excitons) from our analysis, which are, however, accessible using wavefunction-based methods.²² In Fig. 4, we plot the computed vs experimental values for the three approaches for all the studied molecules, for which experimental data are given in Ref. 22. The closer a point lies to the y = x line, the better the agreement between theory and experiment. We fit a linear model to the three sets of results, as a means of visualizing the overall agreement with experiment, from where it becomes evident that the TD-DFT B3LYP+ZPR (red) results provide a significant improvement to the static TD-DFT B3LYP values (green). The numerical results of the Monte Carlo simulations for the exciton energy renormalization are summarized in Sec. S1 of the supplementary material, along with the associated statistical uncertainties. The parameters for all the linear fits are given in Sec. S2 of the supplementary material. Naturally, the Thiel values (black), which were obtained using accurate wavefunction-based methods and capture correlation effects, provide an improvement to static TD-DFT. However, from Fig. 4, it is evident that the correction to TD-DFT exciton energies induced by the quantum fluctuations of molecular vibrations is at least as significant as their correction due to the correlation effects included in the Thiel calculations, at least for the specific set of molecules we study here. From Fig. 4, it also becomes obvious that the ZPR of exciton energies can become very large, with values that can be as high as 1.362 eV in the case of pyrrole. The average correction is found to be (345 ± 67) meV, which is substantial and certainly needs to be accounted for in order for exciton energy calculations to achieve predictive accuracy. We have also investigated the effect of increasing the temperature from 0 to 300 K; however, we found that for the vast majority of molecules, the resulting difference in the exciton energy is small, as summarized in Sec. S3 of the supplementary material. This is intuitively obvious from the fact that small organic molecules are dominated by high-frequency vibrational modes, such as carbon–carbon stretching motions, which generally lie significantly above the threshold for thermal activation at room temperature.

We employ four different statistical measures in order to rigorously compare the accuracy of the different methods presented in Fig. 4, which we summarize in Table I. In particular, for our n = 24 studied molecules with computed exciton energies b_i and experimental exciton energies a_i, we compute the average values of the bias, $bias=1n∑i(bi−ai)$⁠, the root mean-squared error, $RMSE=1n∑i(bi−ai)2$⁠, as well as the relative values to these quantities, $(rel.bias)=1n∑bi−aiai$⁠, $(rel.RMSE)=1N∑i(bi−ai)2ai2$⁠. The reason that we also consider the relative quantities is that the largest errors for the static TD-DFT and Thiel approaches are found in the region of large exciton energies; hence, these molecules dominate the bias and root mean-squared error. For example, the static TD-DFT exciton energy of ethene is renormalized from 8.8 to 7.7 eV once we account for nuclear quantum effects. Using the relative bias and root mean-squared error, which are calculated by dividing by a_i and $ai2$⁠, respectively, we achieve a fairer comparison without under-representing the effect of molecules with lower exciton energies. We see from Table I that the average value of all the employed statistical measures is minimized when using TD-DFT at the B3LYP level combined with the zero-point renormalization (B3LYP+ZPR) of exciton energies due to molecular vibrations. We therefore quantitatively confirm the observation of Fig. 4: for the studied set of molecules, corrections to static exciton energies due to nuclear quantum motion are generally larger than corrections due to beyond-TD-DFT electronic effects that are included in the Thiel values. We also quantify the spread of the four statistical quantities around their average value by computing their standard deviation, which is given in Table S32 of the supplementary material. We find that the spread of the B3LYP+ZPR values is in every case comparable to that of the Thiel values and significantly smaller than that of the static TD-DFT (B3LYP) values.

Similar results to ours were presented in Ref. 12 using the nuclear ensemble method and coupled cluster calculations, although the computed values were not compared to experiment in that case. In Fig. 5, we compare our results obtained within the Monte Carlo method at the TD-DFT level to those of Ref. 12. We restrict the comparison to 14 structures for which the same exciton state has been studied, as Ref. 12 also reports the energies of some excitons that are not necessarily the lowest-lying ones, which we chose to study in this work. We see that both computational methods lead to excellent agreement with experiment, and a rigorous comparison is made in Table II. Among the reported statistical measures, the RMSE and relative RMSE are the ones that most accurately capture the accuracy of the two methods, as they do not suffer from effects such as cancellations of errors, which are evidently present in the bias and relative bias when using TD-DFT due to some values underestimating and some others overestimating experiment. This is not a problem with coupled cluster calculations in the studied molecules, since these generally predict an exciton energy value that is higher compared to TD-DFT, which also becomes apparent from the fact that the values in Ref. 12 tend to slightly overestimate experiment, as seen in Fig. 5. Since the agreement of our results within TD-DFT to experiment is comparable to the agreement of the results of Bai et al., which were obtained using coupled cluster calculations, it is once again highlighted that at least for the specific set of studied molecules, the ZPR of exciton energies provides a larger, or at least comparable, correction to the static TD-DFT exciton energies compared to the correction obtained by using more accurate electronic structure methods. Of course, it is still generally the case that one needs both an accurate description of the electronic structure and of vibrational effects in order to achieve better quality results. To demonstrate this fact, we have chosen four molecules for which TD-DFT calculations at the B3LYP level provide a poor starting point for the static exciton energy, as well as inaccurate values once ZPR is accounted for, and we now perform a Monte Carlo sampling by employing significantly more accurate coupled cluster singles and doubles (CCSD) calculations. These molecules are benzoquinone, cyclopropene, formamide, and pyrrole, and the results are summarized in Table III. In all cases, the corrected exciton energy obtained by using CCSD calculations within our Monte Carlo sampling is closer to the experimental value compared to the TD-DFT case. However, even in this case, the experimental values for benzoquinone and cyclopropene are overestimated by 0.2–0.3 eV. We attribute this behavior to the still limited treatment of correlations within CCSD compared to more accurate methods, which results in these calculations commonly overestimating excited state energies as also observed in the results of Ref. 12 shown in Fig. 5.

Ideally, one would combine accurate electronic structure calculations with our approach for accounting for nuclear quantum motion. However, for most of the studied molecules, we found that a minimum of ∼50 points is required to converge the Monte Carlo sampling of the integral in Eq. (4). While the cost of 50 TD-DFT calculations of small molecules is reasonably low, it quickly becomes very high as one moves to more accurate wavefunction-based methods. It is therefore reasonable to wonder whether the vibrational corrections to exciton energies could be computed at one (cheaper) level of theory and applied to the static exciton energy obtained at another (more accurate) level of theory. Some first conclusions on this can be drawn from the comparison of the excited state energy corrections at the TD-DFT B3LYP and CCSD levels of Table III; with the exception of pyrrole for which the B3LYP correction is 1.5 times larger than the CCSD one, the maximum deviation of the B3LYP zero-point renormalization from the CCSD value is 13%. While, of course, this is encouraging, these data points alone are not enough to draw a general conclusion. However, given this good agreement between the few B3LYP and CCSD excited state energy corrections that we compared, the fact that hybrid functionals such as B3LYP are known to lead to an accurate description of the electron–vibration interactions in organic molecules,³⁸ and also the excellent agreement between theory and experiment in Fig. 4, it seems a reasonable conclusion that hybrid functionals provide an accurate description of exciton–vibration coupling and the zero-point renormalization this induces.

In order to check whether the vibrational correction to exciton energies changes when using different electronic structure methods in a more systematic manner, we repeated our Monte Carlo sampling using the local density approximation (LDA) and pure Hartree–Fock (HF) exchange within TD-DFT while using the vibrational modes obtained at the B3LYP level. This choice of functionals allows us to comment on the role of exact exchange for exciton–vibration interactions: in particular, LDA includes no exact exchange, while B3LYP includes 20% and HF represents the case of including the full exact exchange. Since the included fraction of exchange is generally known to have little effect on the structure of a molecule and the computed ground state vibrational modes,³⁸ using the B3LYP geometries and vibrational modes is a reasonable approximation. Therefore, all differences in the values of the exciton energy renormalization between the different functionals are purely due to the variations of the excited state wavefunction as obtained from TD-DFT.

The numerical results of the calculations employing the different functionals are presented in Sec. S1 of the supplementary material. We compare the LDA and HF energy corrections to the B3LYP ones and find that, on average, LDA predicts a redshift of the excited state energy which is 1.3 times greater than B3LYP, and HF 1.8 times greater. This is in agreement with reports on electron–vibration interactions being strongly affected by the fraction of exact exchange included in the calculation.³⁸ The results suggest that LDA provides a reasonable approximation to B3LYP for computing zero-point renormalization in the studied organic molecules, with only a few structures showing very strong deviations. As shown in Fig. S1 of the supplementary material, the largest differences between the corrections predicted by LDA and B3LYP are found for smaller molecules with less than 60 electrons. We attribute this observation to the fact that the excited states of these smaller structures have a greater electronic density in the vicinity of the localized high-frequency vibrations that can dominate in organics, hence leading to a very strong coupling that is very sensitive to the detailed structure of the excited state wavefunction. These points are revisited in Subsections III C and III E where the impact of molecular size and individual vibrational modes on zero-point renormalization are discussed, respectively.

Overall, we find that in most organic structures, hybrid functionals, such as B3LYP, provide a good approximation to the values of zero-point renormalization computed with more accurate methods, such as CCSD. For several of the studied molecules, even the complete omission of exact exchange through the use of LDA seems to also give reasonable results within 30% of the B3LYP values. Most strong deviations between electronic structure methods are found for smaller molecules, which is in a sense encouraging, since these are the ones that can commonly be studied using more expensive methods without too great of a computational cost.

We found in Ref. 13 for periodic molecular crystals that the correction to exciton energies due to nuclear quantum motion becomes more important for smaller molecules. For the single-molecule systems studied here, we find that the same trend holds, as shown in Fig. 6 for the acene and aliphatic hydrocarbons families of molecules. The results in Fig. 6 are obtained at the TD-DFT B3LYP level of theory. We believe that the reason for the observed trend is that smaller molecules have a greater electronic density in the vicinity of localized atomic motions that are activated due to quantum fluctuations, hence leading to stronger exciton–vibration interactions. Another way of reaching the same conclusion is using Hückel theory. Let us consider the example of linear polyenes, such as ethene, butadiene, hexatriene, and octatetraene, that are studied here and are included among other aliphatic hydrocarbons in Fig. 6(b), showing a decrease in the exciton energy correction with increasing system size. For a linear polyene consisting of N carbon atoms, Hückel theory predicts that the energy of the ith orbital is E_i = α + 2β cos(2πn/(N + 1)), with α, β being the Hückel parameters and i = 1, 2, …, N. From this, we find for the HOMO–LUMO gap,

In the limit of N → ∞, this gives E_gap ∝ 1/N. A similar result can be obtained for cyclic polyenes. Assuming that the same holds for any general conjugated molecule, this suggests that in larger systems, the HOMO is less bonding and the LUMO is less anti-bonding than in comparably smaller systems. Therefore, from the intuitive picture of Fig. 2, the curvature of the ground and excited state surfaces of large molecules is similar, and the difference in their normal mode frequencies is small. Hence, according to the expression (10) that results from the quadratic approximation, the correction to the exciton energy is also smaller than that for a smaller molecule. This is encouraging in the sense that corrections due to nuclear quantum motion are mostly important for small systems, for which they are cheaper to compute.

Having established the accuracy of the Monte Carlo method for computing the exciton energy renormalization due to nuclear quantum motion, we now proceed to use it as a benchmark to assess the accuracy of the quadratic method. Figure 7 visualizes the quadratic vs the Monte Carlo correction for each of the studied molecules, using TD-DFT and the B3LYP functional, along with the cc-pVDZ basis set. The better the agreement between the two methods for a given molecule, the closer the associated point lies to the y = x line (dashed line). Overall, the agreement between the two methods is good, and we propose that the quadratic method can indeed be used to make quantitative predictions of exciton energies, including vibrational renormalization effects. Particularly for smaller molecules, this has the additional advantage that it comes at a lower computational cost compared to Monte Carlo sampling (see also Sec. II and the discussion accompanying Fig. 9). The converged values for the exciton energies as predicted by the quadratic method are summarized in Sec. S5 of the supplementary material.

As already mentioned in Subsection II C, when computing the quadratic correction to the exciton energies, one needs to make a choice for the displacement δu appearing in Eq. (8). We find that generally a displacement of δu = σ, where σ is the standard deviation of the thermal quantum distribution appearing in Eq. (6), leads to results that are in very close agreement with the computed Monte Carlo values. Figures 8(a) and 8(b) show the quadratic correction of pyrazine and tetracene, respectively, comparing the quadratic values at different values of δu (red crosses) with the Monte Carlo correction (blue line) and its associated statistical uncertainty (blue shaded region). The variation of the quadratic correction due to changes in δu is generally comparable to the statistical uncertainty of the Monte Carlo method.

We also explore the computational cost of the quadratic method and compare to that of a Monte Carlo sampling. In Fig. 9, we visualize the central processing unit (CPU) hours required to compute the correction to the static exciton energies that is induced by the zero-point motion of the different molecules. For the Monte Carlo sampling, the calculation always includes taking an average over 100 configurations as mentioned previously, which we find leads to converged results and small statistical uncertainties to the final exciton energies (supplementary material, Sec. S4). Overall, the quadratic method is cheaper to use for smaller structures with fewer electrons, where we also argue that the correction to exciton energies due to vibrations is larger (see Fig. 6 and discussion). Therefore, given also the accuracy of this method and the mode-resolved information it provides (see Subsection III E), we propose that it is the better alternative once one approaches the limit of smaller molecules. For larger structures with several vibrational normal modes, we find that the Monte Carlo method provides a faster estimate of the redshift and should thus be preferred.

While the quadratic method for computing the renormalization of exciton energies due to quantum fluctuations is, in principle, less accurate than performing a Monte Carlo sampling, it is a cheaper and simpler method that is easier to implement computationally. However, perhaps its most important advantage is the fact that it separates the contributions of the individual vibrational normal modes to the exciton energy renormalization, as revealed by Eq. (8). In Figs. 10(b) and 11(b), we visualize the mode-resolved renormalization of the exciton energy of pyrazine and ethene, respectively, at T = 0 K. From these figures, it becomes evident that for both molecules, it is only a few modes that dominate the zero-point renormalization of the exciton energy. In both cases, the most important mode for the exciton energy renormalization is highlighted in a red circle. For pyrazine, this dominant mode is responsible for 35% of the exciton energy renormalization, while for ethene for 56%.

In order to gain a deeper microscopic understanding of the effect of the specific vibrational modes on the exciton states of these molecules, we plot the undistorted “static” exciton wavefunctions (represented by the transition density) in Figs. 10(a) and 11(a), respectively, together with the wavefunction at a typical distorted configuration (δu = σ) along these dominant modes, visualized in Figs. 10(d) and 11(d). For comparison, we also plot the exciton wavefunction when the molecule is displaced along a mode, which is weakly coupled to the exciton in Figs. 10(c) and 11(c), together with the displacement pattern of these motions. These weakly coupled vibrational modes are highlighted in green in Figs. 10(b) and 11(b). It is evident that while the displacement of the dominant modes leads to a significant change of the excitonic wavefunction from its form at the equilibrium geometry, this is not the case for weakly coupled modes.

The dominating modes for exciton renormalization can sometimes be qualitatively explained from a chemically intuitive perspective by considering the conventional skeletal structures of the ground and excited state molecules and [from Eq. (10)] which stretching modes are likely to be far weaker (lower frequency) in the excited state compared to the ground state. For pyrazine in Fig. 10, as an aromatic molecule similar to benzene, it has two principle resonance structures as shown in Fig. 10(e). The HOMO and LUMO are both π orbitals and excitation to S₁ can be qualitatively drawn as breaking a π bond, leading to a (singlet) biradical. The biradical³⁹ has numerous possible resonance structures, three of which are drawn in Fig. 10(e). These quinoidal structures suggest that in the S₁ state, for a given nitrogen, it will be easier to elongate one C–N bond while shortening the other to form a quinoidal-like biradical. This corresponds to the dominant stretching vibration in Fig. 10(d).

Similarly, for ethene in Fig. 11, in the ground S₀ state, the molecule is planar due to the π bond, although this causes some small steric hindrance between hydrogen atoms of different carbons. Upon exciting to the S₁ state, the π bond is broken such that it is much easier to rotate around the central C–C bond, as shown in Fig. 11(e), and we would expect the frequency of this mode to drastically reduce, leading to a large renormalization by Eq. (10). This rotation is the dominating mode as shown in Fig. 11(d).

It is not always true that a subset of modes can be identified to dominate the exciton energy renormalization; however, some patterns can still be observed. Perhaps the most obvious one is a ring-breathing mode of frequency in the vicinity of 1600 cm⁻¹, which is largely responsible for the exciton energy renormalization in several cyclic compounds, which can be intuitively understood as in the case of pyrazine in Fig. 10(e), where this motion contributes 35% of the renormalization, and similarly for other cyclic structures: benzoquinone (20%), pyridine (38%), pyrimidine (39%), tetrazine (62%), and triazine (30%). For each of the studied molecules, we give the percentage of contribution of the mode that most strongly couples to the exciton energy renormalization, along with the frequency of this motion, in Sec. S5 of the supplementary material. Overall, in most molecules, high-frequency modes dominate the renormalization of exciton energies, with the average frequency of the most strongly coupled mode being $ω̄dominant=(1024±121)cm−1$⁠. This dominance of high-frequency modes is consistent with our observation of Subsection III C that zero-point renormalization is for the molecules studied here, more prominent in smaller molecular structures. The excited states of these systems have a greater electronic density in the vicinity of localized high-frequency motions, such as carbon–carbon stretches; hence, the corresponding wavefunctions show stronger variations upon displacement of the atoms. It is for the same reasons that larger differences in the computed energy corrections between the different functionals tend to appear for smaller molecular structures, as seen in Subsection III B; the strong exciton–vibration interactions in these systems can be very sensitive to the local structure of the excited state wavefunctions, which depends on the level of electronic structure theory.

In this work, we have presented an intuitive picture for the redshift of exciton energies that is caused by nuclear quantum motion and two computational approaches for estimating its magnitude. We find that this redshift can indeed be very substantial, reaching values of more than 1 eV in some cases. We rigorously compare our results for the exciton energies of a large set of molecules, obtained within a Monte Carlo framework and time-dependent density functional theory, to experiment and previous computational studies and find that accounting for the effects of nuclear quantum motion in this manner leads to predictive accuracy for exciton energies. We show that the magnitude of the redshift provided by quantum fluctuations is typically larger for smaller molecules, a result that we explain within a Hückel theory picture. Moreover, we find that while the predicted redshift caused by vibrations depends on the level of theory used for the exciton calculations, TD-DFT employing hybrid functionals should provide a good compromise between accuracy and computational cost for estimating its magnitude. Naturally, further improvement of the results can be achieved by combining our Monte Carlo approach with more accurate methods of electronic structure calculations, such as coupled cluster, as we also demonstrate for a subset of the studied molecules. Additionally, we employ a quadratic approximation for computing the magnitude of exciton energy renormalization due to zero-point motion, which allows us to disentangle the contribution of individual normal modes to this effect. While this method is conceptually easier to implement compared to a Monte Carlo sampling, it is, in principle, less accurate as well. Nevertheless, we find that it performs very well for the diverse set of molecules that are studied here, and it also provides a cheaper way of computing the renormalization of exciton energies in smaller molecules, where this effect is also most relevant. By using the quadratic method, we find that for several molecules, the renormalization of the exciton energy is dominated by a few normal modes of vibration, with a ring-breathing motion playing a prominent role in the redshift of the exciton energy in several cyclic compounds. Overall, our study provides critical microscopic insights into the effect of nuclear quantum motion on exciton energies at equilibrium and emphasizes its importance for achieving predictive accuracy.

See the supplementary material for the full results of the Monte Carlo sampling using different functionals, details for the comparison of these results to experiment, and convergence tests and results at 300 K. Detailed data on the vibrational modes dominating the quadratic correction of the exciton energy are also provided.

The authors acknowledge useful discussions with Stuart Althorpe (Cambridge) and thank the anonymous reviewers for their constructive feedback. T.J.H.H. acknowledges a Royal Society University Research Fellowship (Ref. No. URF\R1\201502). B.M. acknowledges support from the Gianna Angelopoulos Programme for Science, Technology, and Innovation. A.M.A. acknowledges support from the Engineering and Physical Sciences Research Council for funding under Grant No. EP/L015552/1 and the Winton Programme for the Physics of Sustainability. Part of the calculations were performed using resources provided by the Cambridge Tier-2 system operated by the University of Cambridge Research Computing Service (http://www.hpc.cam.ac.uk) and funded by EPSRC Tier-2 capital Grant No. EP/P020259/1.

The data that support the findings of this study are openly available in Apollo, the University of Cambridge repository, at https://doi.org/10.17863/CAM.70585.⁴⁰