 View Online  Export Citation RESEARCH ARTICLE | MAY 28 2024 The eXact integral simplified time-dependent density functional theory (XsTD-DFT)  Marc de Wergifosse ; Stefan Grimme J. Chem. Phys. 160, 204110 (2024) https://doi.org/10.1063/5.0206380 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp/article/160/20/204110/3295137/The-eXact-integral-simplified-time-dependent https://pubs.aip.org/aip/jcp/article/160/20/204110/3295137/The-eXact-integral-simplified-time-dependent?pdfCoverIconEvent=cite javascript:; https://orcid.org/0000-0002-9564-7303 javascript:; https://orcid.org/0000-0002-5844-4371 https://crossmark.crossref.org/dialog/?doi=10.1063/5.0206380&domain=pdf&date_stamp=2024-05-28 https://doi.org/10.1063/5.0206380 https://servedbyadbutler.com/redirect.spark?MID=176720&plid=2372063&setID=592934&channelID=0&CID=872267&banID=521836446&PID=0&textadID=0&tc=1&scheduleID=2290748&adSize=1640x440&data_keys=%7B%22%22%3A%22%22%7D&matches=%5B%22inurl%3A%5C%2Fjcp%22%5D&mt=1716965977119991&spr=1&referrer=http%3A%2F%2Fpubs.aip.org%2Faip%2Fjcp%2Farticle-pdf%2Fdoi%2F10.1063%2F5.0206380%2F19968721%2F204110_1_5.0206380.pdf&hc=a236fb4db726e01c370caf5a343d64dbe4945c49&location= The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp The eXact integral simplified time-dependent density functional theory (XsTD-DFT) Cite as: J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 Submitted: 1 March 2024 • Accepted: 5 May 2024 • Published Online: 28 May 2024 Marc de Wergifosse1 ,2,a) and Stefan Grimme2 AFFILIATIONS 1 Theoretical Chemistry Group, Molecular Chemistry, Materials and Catalysis Division (MOST), Institute of Condensed Matter and Nanosciences, Université Catholique de Louvain, Place Louis Pasteur 1, B-1348 Louvain-la-Neuve, Belgium 2Mulliken Center for Theoretical Chemistry, Institut für Physikalische und Theoretische Chemie, Beringstr. 4, 53115 Bonn, Germany Note: This paper is part of the 2024 JCP Emerging Investigators Special Collection. a)Author to whom correspondence should be addressed: marc.dewergifosse@uclouvain.be ABSTRACT In the framework of simplified quantum chemistry methods, we introduce the eXact integral simplified time-dependent density functional theory (XsTD-DFT). This method is based on the simplified time-dependent density functional theory (sTD-DFT), where all semi-empirical two-electron integrals are replaced by exact one- and two-center two-electron integrals, while other approximations from sTD-DFT are kept. The performance of this new parameter-free XsTD-DFT method was benchmarked on excited state and (non)linear response properties, including ultra-violet/visible absorption, first hyperpolarizability, and two-photon absorption (2PA). For a set of 77 molecules, the results from the XsTDA approach were compared to the TDA data. XsTDA/B3LYP excitation energies only deviate on average by 0.14 eV from TDA while drastically cutting computational costs by a factor of 20 or more depending on the energy threshold chosen. The absolute deviations of excitation energies with respect to the full scheme are decreasing with increasing system size, showing the suitability of XsTDA/XsTD-DFT to treat large systems. Comparing XsTDA and its predecessor sTDA, the new scheme generally improves excitation energies and oscillator strengths, in particular, for charge transfer states. TD-DFT first hyperpolarizability frequency dispersions for a set of push-pull π-conjugated molecules are faithfully reproduced by XsTD-DFT, while the previous sTD-DFT method provides redshifted resonance energy positions. Excellent performance with respect to the experiment is observed for the 2PA spectrum of the enhanced green fluorescent protein. The obtained robust accuracy similar to TD-DFT at a fraction of the computational cost opens the way for a plethora of applications for large systems and in high throughput screening studies. Published under an exclusive license by AIP Publishing. https://doi.org/10.1063/5.0206380 I. INTRODUCTION The computation of molecular properties involving the inter- action of light with matter using quantum chemistry (QC) methods remains challenging for large systems that are of timely interest. An ideal method should enable the treatment of open shell systems that may present multi-configurational character1,2 and be able to com- pute excitations in vertical and adiabatic regimes. Considering linear and nonlinear optical response properties, the computation of such higher-order quantities is computationally involved and often needs a proper treatment of electron correlation effects, frequency disper- sion, and explicit account of the environment.3,4 Such methods exist but only for small systems. If used with care, the time-dependent density functional the- ory (TD-DFT)5–8 provides a reasonable route for the computation of excited states and response properties for medium-sized systems. However, it introduces substantial errors for several types of excita- tions involving charge-transfer (CT), double excitation, or Rydberg characters when using standard exchange–correlation (XC) func- tionals.9 Even so, in many applications, it works well for low-lying valence states.9 Considering large systems with a high density of excited states, it is not really necessary to describe individual excita- tions with high accuracy, and a TD-DFT treatment, including only singly excited configurations, might be sufficient. A decade ago, the simplified time-dependent density func- tional theory (sTD-DFT) and its variant using the Tamm–Dancoff J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-1 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp https://doi.org/10.1063/5.0206380 https://pubs.aip.org/action/showCitFormats?type=show&doi=10.1063/5.0206380 https://crossmark.crossref.org/dialog/?doi=10.1063/5.0206380&domain=pdf&date_stamp=2024-May-28 https://doi.org/10.1063/5.0206380 https://orcid.org/0000-0002-9564-7303 https://orcid.org/0000-0002-5844-4371 mailto:mailto:marc.dewergifosse@uclouvain.be https://doi.org/10.1063/5.0206380 The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp approximation (sTDA) were introduced by one of us.10,11 These approaches were extended to range-separated hybrid (RSH) func- tionals by Risthaus et al.12 These methods were originally designed to evaluate the ultra-violet/visible (UV/vis) and circular dichro- ism (CD) spectra of large molecules (<1000 atoms) for which it was not possible to use conventional TD-DFT. In 2016, these methods were extended to ultra-large systems with the extended tight-binding (xTB) variant,13 enabling the computation of excited states for organic systems with up to 5000 atoms and, more recently, to compounds with 4d and 5d metals as well as heavy p-block ele- ments.14 Since then, the reach of these simplified QC (sQC) methods was expanded to response properties allowing the evaluation of the polarizability,15 optical rotation,16 excited state absorption (ESA),17 first hyperpolarizability,15 and two-photon absorption (2PA).18 FIG. 1. Chemical structures of molecules 1 to 36 from the benchmark set of 77 molecules (1–17,11 18–20,12 21–46,13 47–57,34 and 58–7712). J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-2 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp A spin-flip variant19 was also implemented and analysis tools involv- ing natural transition/response orbitals to interpret excited states20 and response properties.21 Here, sTD-DFT/sTDA methods also inspired alternative schemes to speed-up the excited state calculations. In 2016, Rüger et al.22 proposed the TD-DFT+TB method. The TD-DFT+TB scheme, which is also based on a DFT ground state, uses a similar monopole approximation in the linear response treat- ment but considers local XC functionals instead of hybrid ones. Asadi-Aghbolaghi et al.23 showed that the TD-DFT+TB method is up to 100 times faster than TD-DFT to evaluate excited states of large gold and silver plasmonic nanoparticles with less than 0.15 eV of error with respect to TD-DFT for excitation energies. Recently, Havenridge et al.24 extended this method to the efficient compu- tation of the analytical excited state gradient using the Z-vector method. Giannone and Della Sala25 proposed the TD-DFT-as method, TD-DFT using the resolution of identity (RI) with only one s-type Gaussian basis function per atom. They argued that FIG. 2. Chemical structures of molecules 37–77 from the benchmark set of 77 molecules (1–17,11 18–20,12 21–46,13 47–57,34 and 58–7712). J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-3 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp sTD-DFT/sTDA/TD-DFT+TB can be seen as approximations to the TD-DFT-as method, where three-index RI two-electron integrals are replaced by Löwdin approximated ones. The TD-DFT-as approach was recently extended to hybrid XC functionals with the TD-DFT- ris method.26 The use of a minimal auxiliary basis set with TD-DFT reduced the cost by around two orders of magnitude while reproduc- ing the full scheme exceptionally well. Hehn et al.27 added periodic boundary conditions to the sTDA scheme by splitting the Coulomb operator into a semiempirical short-range and an exact long-range contribution. A simplified GW/Bethe–Salpether equation (BSE) was also proposed by Cho et al.28 in 2022. In a perspective article,29 we recently discussed future chal- lenges for sQC methods. One idea was to improve the global accuracy of sQC schemes by going beyond its monopole approxi- mation for the two-electron integrals. We explored this multipole approach but found an alternative route, which is the sub- ject of this work. The sTD-DFT/sTDA approximated molec- ular orbital (MO) two-electron integrals use atomic transition charges obtained from a Löwdin orthogonalization procedure and the Mataga–Nishimoto–Ohno–Klopman (MNOK)30–32 damped Coulomb operator, a function of the interatomic distance. To recover this, instead of using the MNOK operator, Löwdin- orthogonalized two-electron integrals (λαλα∣λβλβ) are approxi- mated by atomic orbital (AO) two-electron integrals (αα∣ββ). By removing the semi-empirical nature of the two-electron integrals, we arrive to a more “ab initio” version of the sTD-DFT method where two-electron integrals in the zero differential overlap (ZDO) approx- imation33 are computed exactly. This new method is introduced as the eXact integral simplified time-dependent density functional theory (XsTD-DFT), and its TDA variant (XsTDA) is also presented. FIG. 3. Chemical structure of further considered systems: carbamazepine, lamotrigine, and diazepam [2,2]-paracyclophane (PCP); six push-pull π-conjugated molecules; four FP chromphores (HBI, HBDI, CFP, and PYP); and the eGFP chromophore with its first shell of residues. J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-4 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp To assess the performance of XsTDA for the evaluation of excited state energies and oscillator strengths, we benchmarked them with respect to the TDA results for a set of 77 molecules taken from earlier studies.11–13,34 Figures 1 and 2 show the chemical struc- tures of these systems. Note that this set includes the DYE12 test set,34 which is composed of 12 large organic dyes and for which back-corrected “experimental” vertical excitation energies are avail- able. This set also contains the mixed CT test set from Risthaus et al.12 for which the reference SCS-CC2/def2-TZVP(-f) excitation energies are available. In addition, UV–vis absorption spectra were computed and compared to the experiment35–37 for the three anti- epileptic drugs, carbamazepine, lamotrigine, and diazepam. The first singlet ESA spectrum of [2,2]-paracyclophane (PCP) computed at the XsTDA/B3LYP level of theory is examined with respect to the EOM-EE-CCSD reference one from Krylov and coworkers.38 The XsTD-DFT/XsTDA methods inherit from all previous implementations in the stda program,39 and thus, we also bench- mark nonlinear optical properties. To assess the performance of the XsTD-DFT method to compute the first hyperpolarizability, six push-pull π-conjugated molecules were employed to assess the per- formance of the XsTD-DFT method. Finally for 2PA, we investigate the spectra for a set of fluorescent protein chromophores40 and for eGFP with a model structure, including the chromophore and its first shell of surrounding residues.41 Figure 3 shows these additional structures. This paper is organized as follows: first, we introduce the the- oretical background for the XsTD-DFT/XsTDA methods for global hybrid XC functionals, and then, we provide computational details and discuss results before providing conclusions and outlooks. II. THEORETICAL BACKGROUND The XsTD-DFT/XsTDA methods are introduced by first recall- ing the density-matrix-based TD-DFT formalism8,42–47 and then showing how simplifications were originally applied to give rise to the sTD-DFT/sTDA framework10,11,13,15–17,29 and finally how to remove its semi-empiricism. In the following, p, q, r, s indices refer to general molecular orbitals (MOs); i, j, k, l to occupied; a, b, c, d to unoccupied molecular orbitals; α, β, γ, δ to atomic orbitals (AOs); and A and B to atoms. Casida’s TD-DFT equations5 are usually used to compute excited states, [(A B B A ) − ω(1 0 0 −1 )](X Y ) = 0 (1) or linear response properties in the presence of a perturbation, [(A B B A ) − ω(1 0 0 −1 )](Xζ(ω) Yζ(ω) ) = −(μζ μζ ), (2) where X and Y are excitation and de-excitation vectors, respec- tively, Xζ(ω) and Yζ(ω) are frequency-dependent linear-response vectors, and the perturbation is μζ,ai = ⟨ϕa∣μ̂ζ ∣ϕi⟩. For a hybrid exchange–correlation functional, A and B supermatrix elements are defined for a given amount ax of Fock exchange as Aia,jb = δijδab(ϵa − ϵi) + 2(ia∣jb) − ax(ij∣ab) + (1 − ax)(ia∣ fXC∣jb) (3) and Bia,jb = 2(ia∣bj) − ax(ib∣aj) + (1 − ax)(ia∣ fXC∣bj), (4) where (ia∣ jb), (ia∣bj), and (ib∣aj) are exchange-type integrals, (ij∣ab) are Coulomb-type two-electron integrals, and (ia∣ fXC∣ jb) and (ia∣ fXC∣bj) are the response of the exchange–correlation functional. To compute excited states using the Tamm–Dancoff approximation, the B supermatrix is neglected, giving rise to a simple Hermitian eigenvalue problem, AX = ωX. (5) The sTD-DFT/sTDA framework applied three approximations to these equations. First, the terms involving the response of the exchange–correlation functional (ia∣ fXC∣ jb) and (ia∣ fXC∣bj) are neglected in the expressions of A and B supermatrices. Second, the expansion space of configuration state functions (CSFs) is truncated considering a single cutoff energy Ethresh. In this procedure, the active MO space is determined for MO energies in between ϵmin = ϵLUMO − 2(1 + 0.8ax)Ethresh (6) and ϵmax = ϵHOMO + 2(1 + 0.8ax)Ethresh. (7) Then, the primary CSFs (P-CSFs) space is selected for i→ a exci- tations when Aia,ia ≤ Ethresh. Remaining j→ b excitations for which Ajb,jb > Ethresh are only selected as secondary CSF (S-CSF) if E(2)jb = P−CSFs ∑ ia ∣Aia,jb∣2 Ajb,jb − Aia,ia > 10−4Eh. (8) The full space of CSFs accounted for in the sTD-DFT/sTDA proce- dures for a given Ethresh is the sum of P-CSFs and S-CSFs. The last approximation regards two-electron integrals that are approximated by a monopole approximation to balance the computational cost and accuracy. MO two-electron integrals should be normally evaluated by the four-index transformation of AO two-electron integrals, (ia∣jb) = ∑ αβγδ C∗iαCaβC∗jγCbδ(αβ∣γδ). (9) This costly transformation can be a computational bottleneck in conventional “ab initio” quantum chemistry. Considering an orthogonalized Löwdin basis λα = ∑β χβS−1/2 βα , MO two-electron integrals are equivalently computed as (ia∣jb) = ∑ αβγδ Clow∗ iα Clow aβ Clow∗ jγ Clow bδ (λαλβ∣λγλδ). (10) Considering the ZDO approximation,33 MO two-electron integrals are approximated as (ia∣jb) ≈∑ αβ Clow∗ iα Clow aα Clow∗ jβ Clow bβ (λαλα∣λβλβ). (11) Assuming that Löwdin-orthogonalized orbitals (LOs) are atom- centered and, thus that the Löwdin-orthogonalized coefficients do J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-5 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp belong to a particular atom, transition charges can be collected for atom A as Qia A =∑ α∈A Clow∗ αi Clow αa . (12) The simplified expression for MO two-electron integrals reads (ia∣jb) ≈∑ AB Qia A Q jb B (AA∣BB), (13) where (AA∣BB) = ∑ α∈A,β∈B (λαλα∣λβλβ). (14) This expression implies that Löwdin-orthogonalized two-electron integrals are available for atoms A and B. Equation (13) presents a simple Coulomb law where (AA∣BB) is approximated by the MNOK30–32 damped Coulomb operator that does not depend on LOs λα anymore. For exchange integrals, it takes the form, (AA∣BB)K = ⎛ ⎝ 1 (RAB)yK + ( ηA+ηB 2 ) −yK ⎞ ⎠ 1 yK , (15) while for Coulomb integrals, it reads ax(AA∣BB)J = ⎛ ⎝ 1 (RAB)yJ + (ax ηA+ηB 2 ) −yJ ⎞ ⎠ 1 yJ , (16) where yK and yJ are globally fitted parameters for the whole range of ax, RAB is the interatomic distance, and ηA is the chemical hardness of atom A. This also means that the onsite integral (AA∣AA) is sim- ply the chemical hardness ηA of A multiplied by ax considering the type of integrals. The MNOK approximation eliminates the depen- dence on LOs of Eq. (14). These expressions are quite beneficial in terms of computational cost by first pre-computing (ia∣BB) =∑ A Qia A(AA∣BB) (17) and then taking the dot product that only scales with the number of atoms when constructing A and B supermatrices, (ia∣jb) ≈∑ B (ia∣BB)Q jb B . (18) This constitutes the basis of the sTD-DFT/sTDA schemes. Now, we introduce the XsTD-DFT/XsTDA methods. It is common in semi-empirical molecular orbital models48 to replace Löwdin-orthogonalized two-electron integrals by AO ones accord- ing to (ia∣jb) ≈∑ αβ Clow∗ iα Clow aα Clow∗ jβ Clow bβ (αα∣ββ). (19) Thus, AO transition charges can be collected for each basis functions as Qia α =∑ α Clow∗ αi Clow αa . (20) Expression (19) becomes (ia∣jb) ≈∑ αβ Qia α Q jb β (αα∣ββ). (21) In the same way as in sTD-DFT/sTDA schemes, its evaluation becomes computationally efficient by pre-computing (ia∣ββ) =∑ α Qia α (αα∣ββ) (22) and then taking the dot product that scales with the number of AOs, (ia∣jb) ≈∑ β (ia∣ββ)Q jb β . (23) In XsTD-DFT/XsTDA, we still consider the original three simpli- fications of the sTD-DFT method, but two-electron integrals to generate A and B supermatrices are now computed using one- and two-center AO two-electron integrals with Eq. (23). Thus, we removed part of the semiempirical nature of the sTD-DFT/sTDA methods. Note that in our implementation, the CSFs space is still deter- mined with the MNOK integrals for efficiency without any notice- able changes in the results. Compared to the sTD-DFT, XsTD-DFT is computationally more involved because the construction of the A and B supermatrices scales with the number of AOs rather than the number of atoms. Figure 4 shows the computation wall times as a function of the number of basis functions to compute the 20 first excited states for the benchmark set of 77 molecules (Figs. 1 and 2) at the B3LYP/6-31+G(d) level of theory on 8 CPUs (Intel Xeon CPU E5-2660 v4, 3.2 GHz). Wall times are also reported for the sB3LYP and XsB3LYP calculations considering FIG. 4. B3LYP, sB3LYP, and XsB3LYP computation wall times as a function of the number of basis functions for the set of 77 molecules. TDA calculations computed 20 excited states, while sTDA and XsTDA ones considered Ethresh. = 10 eV. J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-6 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp Ethresh. = 10 eV performed on an 8 core desktop computer (Intel core i7-6700, 3.4 GHz) with Q-Chem 5.1.49 Typically for Ethresh. = 10 eV, more than 1000 excited states are computed for these systems. Note that the number of excited states computed depends on the system and Ethresh., and thus, it is not possible to give formal scaling factors for simplified methods. As an example, for molecule 50, 1204 and 1121 states were computed at sB3LYP and XsB3LYP levels of theory in 14 s and 3.85 min, respectively, while the full TD-DFT calculation took 26.5 min to compute 20 states. Using a smaller Ethresh. value of 7 eV (stda program39 default), the XsTDA calculation computed “only” 223 excited states in 23 s. This is more than 60 time faster with respect to the full scheme. For compound 35, the TDA calculation took 25.5 min, while sTDA and XsTDA at Ethresh. = 10 eV only took 4 and 11 s, respectively, corresponding to a speed-up of 382 and 139, respectively. For this benchmark set and Ethresh. = 10 eV, the average speed-up is around 100 for sTDA and 20 for XsTDA. These factors can be increased by using a smaller Ethresh. value. III. COMPUTATIONAL DETAILS The XsTD-DFT/XsTDA methods were implemented in a development version of the freely available stda program,39 which is now interfaced with the Libcint library50 to enable the compu- tation of two-electron AO integrals. This implementation benefits from the whole range of linear and nonlinear optical properties already implemented at the sTD-DFT/sTDA levels of theory in the stda program,39 including UV–vis absorption,10,11 CD,10,11 polar- izability,15 optical rotation,16 ESA,17 first hyperpolarizability,15 and two-photon absorption.18 To test the performance of the XsTDA method to compute excited state energies and oscillator strengths, geometries of 77 molecules (Figs. 1 and 2) were taken from previous studies.11–13,34 For 10 and 47–57, back-corrected “experimental” vertical excita- tion energies are taken from Ref. 34. Molecules 69–77 are parts of the mixed CT test set12 for which reference excitation energies for transitions to CT states are available at the SCS-CC2/def2-TZVP(-f) level of theory. Additional RI-CC2/aug-cc-pVDZ excitation energies were computed for 72 for different distances between benzene and TCNE molecular plans. To assess the performance of the XsTDA method to com- pute UV–vis absorption spectra, we first compare it to the exper- imental spectrum of the largest molecule of the benchmark set of 77 molecules, i.e., molecule 50 (perylene-3,4,9,10-tetracarboxylic bisimide) that was taken from Langhals.51 Second, the structure of ferrocene was optimized at the ωB97X-D3/def2-TZVP level of the- ory. Its experimental absorption spectrum was taken from Ref. 10 for comparison. Third, we selected the three drugs, carbamazepine, lam- otrigine, and diazepam, the structures of which were optimized at the ωB97X-D3/6-311G(d) level of theory. The experimental UV–vis absorption spectra were taken from Refs. 35–37, respectively. PCP geometry and the reference EOM-EE-CCSD ESA spectrum were taken from Ref. 38. The reference TD-DFT/TDA excited state calculations were performed with Q-CHEM49 using B3LYP, PBE0, BHandHLYP, and M06-2X functionals and the 6-31+G(d) basis set for the set of 77 molecules. BHandHLYP/6-31+G(d) TD-DFT/TDA calcula- tions were also performed with Q-CHEM49 for carbamazepine, lamotrigine, and diazepam, while sTDA and XsTDA computations employed the same functional/basis set combination. For the set of 77 molecules, an energy threshold of 10 eV was employed to trun- cate the CSFs space, while the default value of 7 eV was used for other compounds. For ferrocene, the 20 first excited states were com- puted at the PBE0/def2-TZVP/TD-DFT/TDA level of theory and compared to the sTDA and XsTDA results using Ethresh. = 7 eV. Note that to simulate intensity borrowing by vibronic coupling, forbidden transition oscillator strengths were set to 5 × 10−4 as it was originally proposed by one of us.10 To assess the performance of the XsTD-DFT method to reproduce frequency-dependent hyper-Rayleigh scattering first hyperpolarizability (βHRS) values, the structures of six push-pull π-conjugated molecules, which were originally used to bench- mark the sTD-DFT method,15 were taken from de Wergi- fosse and Champagne.41 TD-DFT and sTD-DFT BHandHLYP/ 6-31+G(d) βHRS frequency dispersions were taken from our previ- ous study de Wergifosse and Grimme15 The XsTD-DFT βHRS values were computed at the same level. For 2PA, we first used the same set of FP chromophores that we employed to benchmark the sTD-DFT method.18 These struc- tures were taken from Ref. 40 and reference excitation energies as well as 2PA cross sections (σ2PA) at the RI-CC2/d-aug(-d)-cc-pVDZ level for HBI, HBDI, and PYP and with the aug-cc-pVDZ basis set for CFP as well as EOM-CCSD/d-aug(-d)-cc-pVDZ level,52 when available. Second, we took eGFP chromophore and its first shell of surrounding residues from de Wergifosse et al.53 This structure was successfully used to evaluate the 2PA spectrum of eGFP at the sTD- DFT-xTB level of theory.18 The experimental 2PA spectrum was recorded by Drobizhev et al.54 Regarding FP chromophores, TD- DFT 2PA cross sections and excited state energies were computed using Dalton2018.055,56 with the 6-31+G(d) basis set and B3LYP. The sTD-DFT and XsTD-DFT calculations were performed with the same XC functionals and basis set. For eGFP, the B3LYP/6-31G(d) calculations were performed at both the sTD-DFT and XsTD-DFT levels of theory. Input molecular orbitals and their energies for the sTD-DFT and XsTD-DFT calculations used by the stda program39 were com- puted with Q-Chem 5.1.49 All geometry optimizations were also per- formed with Q-Chem.49 RI-CC2 excitation energies were computed with Turbomole 7.6.57 The βHRS values are given in atomic units con- sidering the Taylor series convention,58 while 2PA strengths ⟨δ2PA⟩ are in a.u. and σ2PA in GM (Göppert–Mayer) units. IV. RESULTS AND DISCUSSIONS A. Excitation energies, oscillator strengths, and UV–vis spectra Figure 5 shows the excitation energy and oscillator strength correlation plots for the two lowest-lying singlet excited states of the set of 77 molecules for sTDA and XsTDA with respect to TDA using B3LYP, PBE0, BHandHLYP, and M06-2X XC functionals. For each functional, Table I presents statistical analysis of the excitation energies and oscillator strengths considering regular TDA results as the reference. For oscillator strengths, statistical analysis are given by orders of magnitude (f ∈ [0.0; 0.1[, [0.1; 1.0[, and [1.0; 3.0]). Note that the first range of magnitude includes the order from 0.01 to 0.1 but J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-7 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp FIG. 5. For the two lowest-lying singlet excited states of the set of 77 molecules, correlation graphs for sTDA and XsTDA excitation energies with respect to TDA ones as well as for oscillator strengths using B3LYP, PBE0, BHandHLYP, and M06-2X exchange–correlation functionals and the 6-31+G(d) basis set. Linear regression data are also provided. TABLE I. For B3LYP, PBE0, BHandHLYP, and M06-2X exchange–correlation functionals, mean absolute deviations (MAD), difference between maximum and minimum absolute deviations (AD), and mean deviations (MD) among the set of 77 molecules for sTDA and XsTDA excitation energies and oscillator strengths with respect to TDA ones. The sign of MD indicates a systematic overestimation or underestimation. E (eV) sB3LYP XsB3LYP sPBE0 XsPBE0 sBHandHLYP XsBHandHLYP sM06-2X XsM06-2X MAD 0.16 0.14 0.24 0.16 0.61 0.24 0.74 0.35 max AD - min AD 0.70 0.70 0.92 0.74 2.18 1.42 2.45 1.60 MD −0.15 −0.08 −0.24 −0.10 −0.60 −0.18 −0.74 −0.35 f ∈ [0.0; 0.1[ sB3LYP XsB3LYP sPBE0 XsPBE0 sBHandHLYP XsBHandHLYP sM06-2X XsM06-2X MAD 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 max AD - min AD 0.06 0.11 −0.05 −0.08 0.42 0.33 0.22 0.07 MD +0.00 +0.00 +0.00 +0.00 +0.01 +0.00 +0.00 0.00 f ∈ [0.1; 1.0[ sB3LYP XsB3LYP sPBE0 XsPBE0 sBHandHLYP XsBHandHLYP sM06-2X XsM06-2X MAD 0.04 0.05 0.03 0.04 0.05 0.05 0.10 0.06 max AD - min AD 0.15 0.19 0.15 0.19 0.11 0.19 0.34 0.14 MD −0.01 +0.02 −0.01 +0.03 +0.07 +0.03 −0.09 −0.00 f ∈ [1.0; 3.0] sB3LYP XsB3LYP sPBE0 XsPBE0 sBHandHLYP XsBHandHLYP sM06-2X XsM06-2X MAD 0.09 0.06 0.12 0.05 0.22 0.08 0.27 0.04 max AD - min AD 0.29 0.20 0.33 0.22 0.37 0.26 0.46 0.20 MD −0.07 −0.00 −0.11 +0.02 −0.22 +0.07 −0.27 +0.02 J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-8 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp FIG. 6. Excitation energy absolute deviations of both sTDA and XsTDA methods with respect to TDA calculations as a function of the number of electrons for B3LYP, PBE0, BHandHLYP, and M06-2X functionals for the set of 77 molecules. also lower oscillator strengths for convenience. The statistical anal- ysis includes mean absolute deviations (MAD), differences between maximum and minimum absolute deviations (AD) representing the error spread, and mean deviations (MD) for which the sign indicates systematic overestimation or underestimation. For B3LYP, both sTDA and XsTDA excitation energies devi- ate only slightly with respect to TDA ones with MAD values of 0.16 and 0.14 eV for sTDA and XsTDA, respectively. This shows that the sTDA method is particularly well parameterized for small amounts of Fock exchange as in B3LYP. The error spreads have TABLE II. Excitation energies (eV), oscillator strengths, and types of transition for the two lowest-lying singlet excited states of uracil obtained with BHandHLYP/6-31+G(d) TDA, sTDA, and XsTDA schemes. Uracil TDA sTDA XsTDA ΔE (eV) f Transition ΔE (eV) f Transition ΔE (eV) f Transition 1 5.38 0.00 n→ π∗ 4.52 0.00 π → Ry 5.03 0.00 n→ π∗ 2 5.91 0.29 π → π∗ 5.19 0.00 n→ π∗ 5.36 0.01 π → Ry 3 5.40 0.00 π → Ry 5.86 0.32 π → π∗ 4 5.51 0.24 π → π∗ J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-9 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp the same order of magnitude for both methods. Oscillator strength MADs are small for their respective ranges of magnitude with both the methods. The XsB3LYP scheme slightly outperforms sB3LYP for oscillator strengths larger or equal to 1.0 with respect to the B3LYP results, as clearly shown in the correlation graph (Fig. 5). Considering PBE0 with 25% Fock exchange, the MAD is increased to 0.24 eV for sTDA, while it remains almost constant for XsTDA with 0.16 eV. Increasing the amount of Fock exchange in the functional globally increases the MAD for the sTDA scheme as confirmed by the sBHandHLYP and sM06-2X results. Excitation energy deviations remain smaller with XsTDA than for sTDA. The sM06-2X MAD of 0.74 eV is more than two times larger than the XsM06-2X MAD of 0.35 eV. Similar trends are observed for oscilla- tor strengths especially in the [1.0; 3.0] range for which the XsTDA method clearly outperforms sTDA (Fig. 5). Perylene (9) was one of the sTDA “outliers” pointed out by one of us10 with an absolute energy deviation larger than 0.50 eV at the sPBE0/TZVP level with respect to the experiment. For its first low-lying excited state, the TDA/PBE0/6-31+G(d) excitation energy is 0.37 eV smaller than the experimental value of 3.44 eV, while the XsPBE0 value deviates by only −0.27 eV. This is a significant improvement with respect to the sTDA absolute deviation of 0.60 eV. Figure 6 presents the sTDA and XsTDA excitation energy abso- lute deviations (ADs) with respect to regular TDA as a function of the number of electrons in the molecule. Considering all four func- tionals, ADs are diminishing with increasing number of electrons for both sTDA and XsTDA with respect to TDA. ADs are globally lower for large systems when using the XsTDA method rather than the sTDA one. This supports the suitability of the XsTDA scheme to compute excitation energies for low-lying excited states of large systems. For the DYE12 benchmark set, Tables S1 and S2 compare the computed sTDA and XsTDA excitation energies to TDA as well as to the back-corrected experiment. The XsTDA method reproduced TDA excitation energies well with MADs of about 0.10 eV. With respect to the back-corrected experiment, the best result is provided by the XsB3LYP scheme with a MAD of 0.23 eV, while TDA yields a value of 0.24 eV. Note that the best sTDA comparison with respect to the experiment is obtained with the BHandHLYP functional with a MAD of 0.23 eV. Clearly, the sTDA parametrization and error can- cellations play a role here, as indicated by the TDA MAD of 0.37 eV with respect to the experiment. As originally pointed out by one of us,10 while π → π∗ tran- sitions are relatively well-treated by sTDA, larger discrepancies are expected for n→ π∗ and d → d transitions due to the applied monopole approximation, which fails for such cases. To illus- trate this behavior and how the XsTDA scheme deals with these types of excitations, we first selected uracil (4) that presents a dipole-forbidden n→ π∗ excitation and a π → π∗ transition as the two first low-lying singlet excited states. Table II presents the excited state energies and oscillator strengths, as obtained by the BHandHLYP/6-31+G(d) TDA, sTDA, and XsTDA schemes. First, both sTDA and XsTDA schemes yield π → Ry unwanted low-lying “ghost” states that are not present in the TDA calcula- tion. Second, the π → π∗ transition is better treated by the XsTDA method with only −0.05 eV of deviation with respect to the full scheme, while the sTDA error is −0.40 eV. For, the n→ π∗ transi- tion, relatively large deviations are observed with both sTDA and XsTDA schemes (−0.19 and −0.35 eV). As stressed by one of us,10 n→ π∗ transition are badly described by sTD-DFT/sTDA because of wrong localized exchange-type contributions due to the monopole approximation. This is inherent to any ZDO-based methods, such as the XsTD-DFT/XsTDA formalism, since one-center exchange inte- grals [(αβ∣αβ), α and β ∈ A] are excluded.59 Thus, the XsTD-DFT/ XsTDA formalism presents the same problem as the sTD- DFT/sTDA scheme to describe local n→ π∗ transitions. A possi- ble correction could be to use the INDO (intermediate neglect of differential overlap) approximation60 instead of the ZDO one. As a second illustrative case, we selected the ferrocene molecule that was already used in the original sTDA publication10 to illus- trate the relative poor treatment of d → d transitions by the sTDA method. Figure 7 shows the experimental absorption spectrum of ferrocene in comparison with the PBE0/def2-TZVP TDA, sTDA, and XsTDA spectra. The experimental spectrum presents three main absorption bands where the two lowest ones are related to dipole-forbidden metal-centered d → d transitions. These two bands are poorly described by the sTDA method, which are significantly improved with the XsTDA scheme. The rest of the XsTDA spectrum is globally blueshifted with respect to the experiment. Treating electronic transitions to Rydberg states correctly is probably out of the scope of XsTD-DFT/XsTDA schemes. Neverthe- less, among the test set, we have the first one photon-active transition n→ Ry(3s) of acetone (3) with an experimental excitation energy of 6.35 eV and a reference theoretical oscillator strength of 0.037 (MRCI1). At the B3LYP/6-31+G(d) level, the experimental excita- tion energy is largely overestimated with a value of 6.91 eV, but the FIG. 7. Experimental absorption spectrum of ferrocene in comparison with the PBE0/def2-TZVP TDA, sTDA, and XsTDA spectra. Note that the oscilla- tor strengths for forbidden transitions are set to 5 × 10−4 to simulate vibronic couplings. J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-10 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp FIG. 8. Excitation energies (left panel) computed at TDA, sTDA, and XsTDA levels using B3LYP, PBE0, BHandHLYP, and M06-2X exchange–correlation functionals and the 6-31+G(d) basis set in comparison with reference SCS-CC2/def2-TZVP(-f) values for the mixed CT test set.12 The mean absolute deviations (right panel) for all these levels of theory with respect to SCS-CC2. MRCI oscillator strength is well reproduced with a value of 0.035. Only one pair of natural transition orbitals (NTOs) contributes to this transition, confirming its n→ Ry(3s) character (see Fig. S2). The treatment is surprisingly improved with the XsB3LYP scheme giving an excitation energy of 6.33 eV and f = 0.037. The shape of NTOs at this level are similar to B3LYP ones. The sB3LYP level of theory is also providing a better energy value with respect to the full scheme with 6.46 eV ( f = 0.035). Hybrid XC functionals do not provide a good description for electron transitions to CT states. Usually, increasing the amount of Fock exchange in the functionals improves the treatment, but still large deviations are observed.10,61 Employing RSH functionals is the usual solution for a better treatment of CT states but is not yet implemented with XsTD-DFT/XsTDA schemes. To assess the performance of the XsTDA scheme to treat CT states, the mixed CT test set from Risthaus et al.12 was taken. Figure 8 shows the comparison between the TDA, sTDA, and XsTDA excitation ener- gies to CT states using B3LYP, PBE0, BHandHLYP, and M06-2X exchange–correlation functionals and 6-31+G(d) basis set and ref- erence SCS-CC2/def2-TZVP(-f) values. MADs for these excitation energies with respect to SCS-CC2 are also reported. As expected, the XC functional with a large amount of HF exchange agrees best with the reference. MADs for BHandHLYP and M06-2X are 0.56 and 0.48 eV, respectively. XsTDA reproduces the TDA results almost perfectly, while sTDA on the contrary presents larger MADs of 1.0 eV (sBHandHLYP) and 0.98 eV, respectively. Regarding inter- molecular CT states for systems 72–77, all reference excitation ener- gies are systematically underestimated by TDA, a behavior almost perfectly reproduced by the XsTDA scheme. Molecule 72 is composed of one benzene unit and one TCNE unit, for which their molecular planes are parallel to each other. Figure 9 shows the first low-lying excited state energy of this sys- tem shown in the inset of Fig. 9 for NTOs as a function of the distance (R) between both molecular planes evaluated at the TDA, sTDA, and XsTDA [BHandHLYP/6-31+G(d)] levels of theory with respect to the reference RI-CC2/aug-cc-pVDZ. A reference curve is also provided and computed as ECT = IP′ − EA′ − 1/R, where IP′ refers to the KS-DFT HOMO energy of benzene and EA′, the KS-DFT LUMO energy of TCNE. All TDA, sTDA, and XsTDA FIG. 9. Comparison of CT excitation energies computed at TDA, sTDA, and XsTDA [BHandHLYP/6-31+G(d)] levels of theory with respect to the RI-CC2/aug-cc-pVDZ reference for molecule 72 as a function of the distance (R) between both parallel molecular plan of benzene and TCNE. A reference analytical single-particle curve is also provided. The inset presents the natural transition orbitals obtained at the TDA/BHandHLYP/6-31+G(d) level of theory for the first low-lying state. Only one pair contributes to the transition with a weight of 0.99. J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-11 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp FIG. 10. UV–vis experimental absorption spectra35–37 of carbamazepine in methanol, lamotrigine, and diazepam in water as well as BHandHLYP/6-31+G(d) TDA, sTDA, and XsTDA spectra in vacuum. results asymptotically tend to an excitation energy value of around 1.3 eV lower than the RI-CC2 one. This is expected for hybrid XC functionals because of the integer discontinuity problem and the self-interaction error.10,62–64 At large distance, this CT state is com- posed of benzene+ and TCNE−, two charges with a 1/R interaction energy. In TD-DFT, this interaction is described by the −(ii∣aa) term in Eq. (1) but scaled inappropriately by ax, hence not providing the expected 1/R decay. By design, the sTDA scheme was corrected for this behavior by asymptotically approaching the 1/R contribu- tion directly into the semi-empirical integral [see Eq. (16)].10 This is illustrated in Fig. 9 showing that the sTDA scheme reproduces the ECT = IP′ − EA′ − 1/R curve well, while the XsTDA scheme matches almost perfectly with the TDA results. Figure 10 shows the experimental absorption spectra35–37 of three anti-epileptic drugs, carbamazepine, lamotrigine, and diazepam, recorded in methanol for the carbamazepine and in water for lamotrigine and diazepam. They are compared to the BHandHLYP/6-31+G(d) TDA, sTDA, and XsTDA vertical spectra. FIG. 11. Comparison of the first singlet ESA spectrum for PCP computed at XsB3LYP/TDA/6–311++G(d,p) to reference the EOM-EE-CCSD/6–311++G(d,p) spectra for both x and (y + z) polarizations. An energy shift of 0.9 eV was applied to the XsB3LYP spectrum. As it was observed for the benchmark set, sBHandHLYP/TDA exci- tation energies are redshifted with respect to the BHandHLYP/TDA ones, which is corrected by XsTDA and the TDA spectra are well reproduced. However, both the TDA and XsTDA spectra are sys- tematically blueshifted with respect to the experiment because of a missing solvent, relaxation, and vibrational effects. PCP is a model for studying benzene excimer. Figure 11 shows its first singlet ESA spectra computed at both EOM-EE- CCSD and XsB3LYP levels of theory. The XsB3LYP reproduces well the intensity of the main absorption band around 2.38 eV with respect to the EOM-EE-CCSD one. As it could be expected using the B3LYP XC functional, some features of the reference spectrum are not present, such as some lower energy peaks below 2 eV. Nevertheless, the comparison is rather good for such low-level treatment. B. Second-harmonic generation of push-pull π-conjugated systems One recurring problem with the sTD-DFT scheme to evaluate βHRS is the systematically redshifted resonance enhancements with respect to the full scheme for push-pull π-conjugated molecules. To circumvent this, we advocated refitting the MNOK two-electron integral parameters.15,29,65 Figure 12 shows the βHRS frequency dispersions obtained at BHandHLYP/6-31+G(d) TD-DFT, sTD-DFT, and XsTD-DFT levels of theory for six push-pull π-conjugated molecules. The XsTD-DFT method reproduces TD-DFT βHRS frequency disper- sions almost perfectly and, in particular, the energy position of the resonance enhancement, clearly outperforming the sTD-DFT scheme. This is easily explainable as these resonance enhancements are due to the CT states for which XsTD-DFT is particularly good at reproducing the TD-DFT excitation energies. The computational cost to evaluate the first hyperpolarizability at the XsTD-DFT level is very low with respect to the full scheme. Computing one single βHRS value for m-3 with GAMESS66,67 took 30.9 min15 on an 8 cores desk- top computer (Intel core i7-6700, 3.40 GHz), while the XsTD-DFT calculation finished in only 5.2 min for 13 βHRS values. The sTD-DFT calculation took 1.9 min. J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-12 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp FIG. 12. βHRS frequency dispersions obtained at BHandHLYP/6-31+G(d) TD-DFT, sTD-DFT, and XsTD-DFT levels of theory for the six push-pull π-conjugated molecules. For sTD-DFT and XsTD-DFT calculations, Ethresh. = 15 eV was used. C. Two-photon absorption Recently, we showed that calculating absolute 2PA cross sec- tions [depending on (2ω)2 and the half width at half maximum for a given line shape] is inherently difficult for TD-DFT and de facto for simplified methods based on TD-DFT.18 The error for the excitation energy directly impacts the 2PA cross section. Regard- ing 2PA strengths [independent of (2ωge)2] among a set of push- pull π-conjugated systems, Zaleśny and coworkers68 demonstrated that global hybrid functionals are not well reproducing trends but have better absolute values, while RSH functionals underestimate J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-13 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp FIG. 13. 2PA spectra for HBI, HBDI, CFP, and PYP computed at B3LYP/6-31+G(d) TD-DFT, sTD-DFT, and XsTD-DFT in comparison to reference RI-CC2/d-aug(-d)-cc- pVDZ or aug-cc-pVDZ excitation energies and 2PA cross sections for the lowest energy excited state from Beerepoot et al.40 and EOM-CCSD/d-aug(-d)-cc-pVDZ reference calculations from 52, when available. RI-CC2 2PA strengths but better correlate with them. RSH func- tionals will be implemented soon at the XsTD-DFT level. Note that currently the XsTD-DFT method is extensively benchmarked for 2PA on the benchmark set from Zaleśny and coworkers68 as well as on compounds from the QUEST database69–71 at UCLouvain. Figure 13 compare the B3LYP TD-DFT, sTD-DFT, and XsTD- DFT 2PA spectra to RI-CC2 and EOM-CCSD excitation energy and σ2PA for the lowest excited state of the four FP chromophores HBI, HBDI, CFP, and PYP. Ethresh = 9 eV was employed for the simplified calculations, and a logarithmic scale was used to highlight small 2PA cross sections. Figure S1 shows that 2PA spectra are almost identical using larger basis sets (aug-cc-pVDZ and d-aug-cc-pVDZ), and thus, the 6-31+G(d) basis set was selected. TD-DFT 2PA spectral shapes are globally well-reproduced by both simplified methods. The sTD-DFT 2PA spectra are slightly redshifted with respect to the TD-DFT ones, while XsTD-DFT excitation energies are blueshifted. Considering the first bright 2P-induced excitation energy for the four systems with respect to the RI-CC2 results, MAD of 0.22, 0.26, and 0.09 eV are obtained for TD-DFT, sTD-DFT, and XsTD-DFT, respectively. Comparing orders of magnitudes of 2PA strengths with respect to RI-CC2 ones, MADs for log σ2PA are 0.34 (B3LYP), 0.18 (sB3LYP), and 0.24 (XsB3LYP). This shows that the XsB3LYP excitation ener- gies are in very good agreement with the RI-CC2 ones, but 2PA strengths are slightly better reproduced by the sB3LYP scheme. Note that for the first bright 2PA excitation of HBDI and PYP, all the results are redshifted with respect to EOM-CCSD excitation energies. Figure 14 shows the experimental 2PA spectrum recorded by Drobizhev et al.54 as well as the 2PA B3LYP/6-31G(d) sTD-DFT and XsTD-DFT spectra obtained for an eGFP model system (the depro- tonated chromophore and its first shell of surrounding residues with 359 atoms). This model system was successfully used to reproduce the experimental 2PA spectrum at the sTD-DFT-xTB level of the- ory.18 Note that the XsB3LYP 2PA spectrum was globally redshifted by 0.60 eV. Both sB3LYP and XsB3LYP spectra compare very well with the experiment. J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-14 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp FIG. 14. Experimental 2PA spectrum from Drobizhev et al.54 in comparison with the sB3LYP and XsB3LYP ones computed with the 6-31G(d) basis set, Ethresh = 7 eV and N = 2 [see Eq. (16) Ref. 18]. An energy shift of −0.60 eV was applied to the XsTD-DFT spectrum. V. CONCLUSIONS This contribution introduced XsTD-DFT, a computationally efficient method designed to calculate linear and nonlinear optical properties of large systems. It is based on the sTD-DFT scheme for which semi-empirical integrals are replaced by exact one- and two- center AO two-electron integrals. This new parameter-free method still uses the ZDO approximation with Löwdin-orthogonalized LCAO coefficients to compute MO two-electron integrals. This modification still keeps the XsTD-DFT method computationally efficient with respect to TD-DFT by drastically cutting computa- tional costs by a factor of 20 or more depending on the single energy threshold used. The main difference appears in the construction of A and B TD-DFT supermatrices, where the XsTD-DFT method scales with the number of AOs instead of the number of atoms. Compared to sTD-DFT, these changes were implemented to robustly repro- duce the TD-DFT results. A TDA variant named XsTDA was also implemented. Both schemes are currently limited to global hybrid XC functionals but can be extended to the RSH level. The XsTDA scheme was tested for vertical excitation energies of a set of 77 molecules and compared to the conventional TDA results using B3LYP, PBE0, BHandHLYP, and M06-2X XC func- tionals. With respect to the previous sTDA scheme, XsTDA is more robust at reproducing TDA results especially when using XC func- tionals with large amount of Fock exchange. Oscillator strengths are also improved with XsTDA, especially for large f values. Interestingly, we observed that excitation energy absolute deviations decrease with increasing size of the system, a behavior most wel- comed for a method designed to treat large systems. In 2013, one of us10 showed that sTDA performs badly when computing n→ π∗ and metal-centered d → d transitions. The treatment of n→ π∗ tran- sitions is not much improved by the XsTDA scheme because of the inherent ZDO approximation that neglects one-center exchange integrals. An INDO treatment might solve this problem. Regarding d → d transitions, the XsTDA method was able to predict the lowest energy metal-centered d → d excitations for ferrocene. Surprisingly, the first Rydberg-type n→ Ry(3s) transition of acetone (3) was well reproduced by the XsTDA method. Regarding CT states, the XsTDA method robustly reproduces the TDA results, outperforming sTDA, especially for a large amount of Fock exchange. The XsTD-DFT scheme performs equivalently well for NLO properties. For a set push-pull π-conjugated systems, we showed that the sTD-DFT method systematically redshifts βHRS resonance enhancements in comparison with the full scheme.15,65 The XsTD- DFT method corrects this behavior thanks to its consistent treat- ment of CT states with respect to TD-DFT. Both simplified schemes reproduce the TD-DFT 2PA spectra well for a set of FP chro- mophores. Note that the XsTD-DFT method provided the best comparison with respect to RI-CC2 excitation energies for the first bright 2PA excited states. Finally, the experimental 2PA spectrum of eGFP was well reproduced by both sB3LYP and XsB3LYP methods using a structural model system. In summary, we can conclude that the XsTD-DFT/XsTDA methods are robust alternatives to TD-DFT/TDA to compute at a fraction of their computational costs excited state and nonlinear optical response properties. With these new methods, all-atom QC calculations of large systems can now be performed with a similar accuracy as the full scheme. Such methods can also be useful for high throughput screening studies considering large ensembles of molecules. To complement the range of applications of Xs methods, extra implementations are currently ongoing at UCLouvain. The XsTDA/XsTD-DFT excited state gradient is currently being imple- mented and support for range-separated hybrids will be available soon. The extension of XsTD-DFT/XsTDA schemes to compute core to valence excitations is an exciting prospect for which we have already gathered interesting results. An extension of XsTD-DFT/ XsTDA schemes to use tight-binding ground states for very large systems will be investigated in a near future. A TD-DFT-ris ver- sion of XsTD-DFT may be investigated to further reduce the computational cost. SUPPLEMENTARY MATERIAL The supplementary material provides basis set effects on the HBI 2PA spectrum, NTOs for the first one photon active elec- tronic transition of acetone, and for B3LYP, PBE0, BHandHLYP, and M06-2X exchange–correlation functionals, excitation energies for the DYE12 benchmark set and absolute deviations at both sTDA and XsTDA levels from TDA. J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-15 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp https://doi.org/10.60893/figshare.jcp.c.7218456 The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp ACKNOWLEDGMENTS M. de W. dedicates this work to his late, beloved partner Coralie Leroy. The authors thank Marilù G. Maraldi for performing the RI-CC2 calculations. This work was supported by the DFG in the framework of the project “Theoretical studies of nonlinear opti- cal properties of fluorescent proteins by novel low-cost quantum chemistry methods” (Grant No. 450959503). AUTHOR DECLARATIONS Conflict of Interest The authors have no conflicts to disclose. Author Contributions Marc de Wergifosse: Conceptualization (lead); Funding acquisi- tion (equal); Investigation (lead); Methodology (lead); Software (lead); Writing – original draft (lead); Writing – review & editing (equal). Stefan Grimme: Funding acquisition (equal); Methodol- ogy (supporting); Software (supporting); Writing – original draft (supporting); Writing – review & editing (equal). DATA AVAILABILITY The data that support the findings of this study are available from the corresponding author upon reasonable request. REFERENCES 1S. Grimme, “Calculation of the electronic spectra of large molecules,” in Reviews in Computational Chemistry (John Wiley & Sons, Ltd, 2004), Chap. III, pp. 153–218. 2S. Ghosh, P. Verma, C. J. Cramer, L. Gagliardi, and D. G. Truhlar, Chem. Rev. 118, 7249 (2018). 3D. M. Bishop and P. Norman, “Nonlinear optical materials,” in Handbook of Advanced Electronic and Photonic Materials and Devices, edited by H. S. Nalwa (Academic Press, Burlington, 2001). 4M. G. Papadopoulos, A. J. Sadlej, and J. Leszczynski, Non-linear Optical Properties of Matter (Springer, 2006). 5M. E. Casida, Time-Dependent Density Functional Response Theory for Molecules (World Scientific, 1995), pp. 155–192. 6E. K. U. Gross, J. F. Dobson, and M. Petersilka, “Density functional theory of time-dependent phenomena,” in Density Functional Theory II: Relativistic and Time Dependent Extensions, edited by R. F. Nalewajski (Springer Berlin Heidelberg, Berlin, Heidelberg, 1996), pp. 81–172. 7R. Bauernschmitt and R. Ahlrichs, Chem. Phys. Lett. 256, 454 (1996). 8F. Furche, J. Chem. Phys. 114, 5982 (2001). 9C. Adamo and D. Jacquemin, Chem. Soc. Rev. 42, 845 (2013). 10S. Grimme, J. Chem. Phys. 138, 244104 (2013). 11C. Bannwarth and S. Grimme, Comput. Theor. Chem. 1040–1041, 45 (2014). 12T. Risthaus, A. Hansen, and S. Grimme, Phys. Chem. Chem. Phys. 16, 14408 (2014). 13S. Grimme and C. Bannwarth, J. Chem. Phys. 145, 054103 (2016). 14J. Seibert, J. Pisarek, S. Schmitz, C. Bannwarth, and S. Grimme, Mol. Phys. 117, 1104 (2019). 15M. de Wergifosse and S. Grimme, J. Chem. Phys. 149, 024108 (2018). 16M. de Wergifosse, J. Seibert, and S. Grimme, J. Chem. Phys. 153, 084116 (2020). 17M. de Wergifosse and S. Grimme, J. Chem. Phys. 150, 094112 (2019). 18M. de Wergifosse, P. Beaujean, and S. Grimme, J. Phys. Chem. A 126, 7534 (2022). 19M. de Wergifosse, C. Bannwarth, and S. Grimme, J. Phys. Chem. A 123, 5815 (2019). 20M. de Wergifosse, J. Seibert, B. Champagne, and S. Grimme, J. Phys. Chem. A 123, 9828 (2019). 21M. de Wergifosse and S. Grimme, J. Chem. Theory Comput. 16, 7709 (2020). 22R. Rüger, E. van Lenthe, T. Heine, and L. Visscher, J. Chem. Phys. 144, 184103 (2016). 23N. Asadi-Aghbolaghi, R. Rüger, Z. Jamshidi, and L. Visscher, J. Phys. Chem. C 124, 7946 (2020). 24S. Havenridge, R. Rüger, and C. M. Aikens, J. Chem. Phys. 158, 224103 (2023). 25G. Giannone and F. Della Sala, J. Chem. Phys. 153, 084110 (2020). 26Z. Zhou, F. Della Sala, and S. M. Parker, J. Phys. Chem. Lett. 14, 1968 (2023). 27A.-S. Hehn, B. Sertcan, F. Belleflamme, S. K. Chulkov, M. B. Watkins, and J. Hutter, J. Chem. Theory Comput. 18, 4186 (2022). 28Y. Cho, S. J. Bintrim, and T. C. Berkelbach, J. Chem. Theory Comput. 18, 3438 (2022). 29M. de Wergifosse and S. Grimme, J. Phys. Chem. A 125, 3841 (2021). 30K. Nishimoto and N. Mataga, Z. Phys. Chem. 12, 335 (1957). 31K. Ohno, Theor. Chim. Acta 2, 219 (1964). 32G. Klopman, J. Am. Chem. Soc. 86, 4550 (1964). 33R. Pariser and R. G. Parr, J. Chem. Phys. 21, 466 (1953). 34L. Goerigk and S. Grimme, J. Chem. Phys. 132, 184103 (2010). 35N. Zadbuke, S. Shahi, A. Jadhav, and S. Borde, Int. J. Pharm. Pharm. Sci. 8, 234 (2016). 36O. S. Keen, I. Ferrer, E. Michael Thurman, and K. G. Linden, Chemosphere 117, 316 (2014). 37S. K. Fadjimata, A. Rabani, and J. IOSR, Appl. Chem 12, 51 (2019). 38O. Haggag, R. Baer, S. Ruhman, and A. I. Krylov, J. Chem. Phys. 160, 124111 (2024). 39The STDA code can be obtained from: https://github.com/grimme-lab/stda/. 40M. Beerepoot, D. Friese, N. List, J. Kongsted, and K. Ruud, Phys. Chem. Chem. Phys. 17, 19306 (2015). 41M. de Wergifosse and B. Champagne, J. Chem. Phys. 134, 074113 (2011). 42A. Zangwill and P. Soven, Phys. Rev. A 21, 1561 (1980). 43S. Hirata and M. Head-Gordon, Chem. Phys. Lett. 314, 291 (1999). 44F. Wang, C. Y. Yam, and G. Chen, J. Chem. Phys. 126, 244102 (2007). 45A. J. Thorvaldsen, K. Ruud, K. Kristensen, P. Jørgensen, and S. Coriani, J. Chem. Phys. 129, 214108 (2008). 46F. Zahariev and M. S. Gordon, J. Chem. Phys. 140, 18A523 (2014). 47S. M. Parker, D. Rappoport, and F. Furche, J. Chem. Theory Comput. 14, 807 (2018). 48T. Husch, A. C. Vaucher, and M. Reiher, Int. J. Quantum Chem. 118, e25799 (2018). 49E. Epifanovsky, A. T. B. Gilbert, X. Feng, J. Lee, Y. Mao, N. Mardirossian, P. Pokhilko, A. F. White, M. P. Coons, A. L. Dempwolff, Z. Gan, D. Hait, P. R. Horn, L. D. Jacobson, I. Kaliman, J. Kussmann, A. W. Lange, K. U. Lao, D. S. Levine, J. Liu, S. C. McKenzie, A. F. Morrison, K. D. Nanda, F. Plasser, D. R. Rehn, M. L. Vidal, Z.-Q. You, Y. Zhu, B. Alam, B. J. Albrecht, A. Aldossary, E. Alguire, J. H. Andersen, V. Athavale, D. Barton, K. Begam, A. Behn, N. Bellonzi, Y. A. Bernard, E. J. Berquist, H. G. A. Burton, A. Carreras, K. Carter-Fenk, R. Chakraborty, A. D. Chien, K. D. Closser, V. Cofer-Shabica, S. Dasgupta, M. de Wergifosse, J. Deng, M. Diedenhofen, H. Do, S. Ehlert, P.-T. Fang, S. Fatehi, Q. Feng, T. Friedhoff, J. Gayvert, Q. Ge, G. Gidofalvi, M. Goldey, J. Gomes, C. E. Gonzalez-Espinoza, S. Gulania, A. O. Gunina, M. W. D. Hanson-Heine, P. H. P. Harbach, A. Hauser, M. F. Herbst, M. Hernandez Vera, M. Hodecker, Z. C. Holden, S. Houck, X. Huang, K. Hui, B. C. Huynh, M. Ivanov, A. Jasz, H. Ji, H. Jiang, B. Kaduk, S. Kahler, K. Khistyaev, J. Kim, G. Kis, P. Klunzinger, Z. Koczor-Benda, J. H. Koh, D. Kosenkov, L. Koulias, T. Kowalczyk, C. M. Krauter, K. Kue, A. Kunitsa, T. Kus, I. Ladjanszki, A. Landau, K. V. Lawler, D. Lefrancois, S. Lehtola, R. R. Li, Y.-P. Li, J. Liang, M. Liebenthal, H.-H. Lin, Y.-S. Lin, F. Liu, K.-Y. Liu, M. Loipersberger, A. Luenser, A. Manjanath, P. Manohar, E. Mansoor, S. F. Manzer, S.-P. Mao, A. V. Marenich, T. Markovich, S. Mason, S. A. Maurer, P. F. McLaughlin, M. F. S. J. J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-16 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp https://doi.org/10.1021/acs.chemrev.8b00193 https://doi.org/10.1016/0009-2614(96)00440-x https://doi.org/10.1063/1.1353585 https://doi.org/10.1039/c2cs35394f https://doi.org/10.1063/1.4811331 https://doi.org/10.1016/j.comptc.2014.02.023 https://doi.org/10.1039/c3cp54517b https://doi.org/10.1063/1.4959605 https://doi.org/10.1080/00268976.2018.1510141 https://doi.org/10.1063/1.5037665 https://doi.org/10.1063/5.0020543 https://doi.org/10.1063/1.5080199 https://doi.org/10.1021/acs.jpca.2c02395 https://doi.org/10.1021/acs.jpca.9b03176 https://doi.org/10.1021/acs.jpca.9b08474 https://doi.org/10.1021/acs.jctc.0c00990 https://doi.org/10.1063/1.4948647 https://doi.org/10.1021/acs.jpcc.0c00979 https://doi.org/10.1063/5.0142240 https://doi.org/10.1063/5.0020545 https://doi.org/10.1021/acs.jpclett.2c03698 https://doi.org/10.1021/acs.jctc.2c00144 https://doi.org/10.1021/acs.jctc.2c00087 https://doi.org/10.1021/acs.jpca.1c02362 https://doi.org/10.1524/zpch.1957.12.5_6.335 https://doi.org/10.1007/bf00528281 https://doi.org/10.1021/ja01075a008 https://doi.org/10.1063/1.1698929 https://doi.org/10.1063/1.3418614 https://doi.org/10.1016/j.chemosphere.2014.07.085 https://doi.org/10.9790/5736-1212015158 https://doi.org/10.1063/5.0196641 https://github.com/grimme-lab/stda/ https://doi.org/10.1039/c5cp03241e https://doi.org/10.1039/c5cp03241e https://doi.org/10.1063/1.3549814 https://doi.org/10.1103/physreva.21.1561 https://doi.org/10.1016/s0009-2614(99)01149-5 https://doi.org/10.1063/1.2746034 https://doi.org/10.1063/1.2996351 https://doi.org/10.1063/1.2996351 https://doi.org/10.1063/1.4867271 https://doi.org/10.1021/acs.jctc.7b01008 https://doi.org/10.1002/qua.25799 The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp Menger, J.-M. Mewes, S. A. Mewes, P. Morgante, J. W. Mullinax, K. J. Oosterbaan, G. Paran, A. C. Paul, S. K. Paul, F. Pavosevic, Z. Pei, S. Prager, E. I. Proynov, A. Rak, E. Ramos-Cordoba, B. Rana, A. E. Rask, A. Rettig, R. M. Richard, F. Rob, E. Rossomme, T. Scheele, M. Scheurer, M. Schneider, N. Sergueev, S. M. Sharada, W. Skomorowski, D. W. Small, C. J. Stein, Y.-C. Su, E. J. Sundstrom, Z. Tao, J. Thirman, G. J. Tornai, T. Tsuchimochi, N. M. Tubman, S. P. Veccham, O. Vydrov, J. Wenzel, J. Witte, A. Yamada, K. Yao, S. Yeganeh, S. R. Yost, A. Zech, I. Y. Zhang, X. Zhang, Y. Zhang, D. Zuev, A. Aspuru-Guzik, A. T. Bell, N. A. Besley, K. B. Bravaya, B. R. Brooks, D. Casanova, J.-D. Chai, S. Coriani, C. J. Cramer, G. Cserey, A. E. DePrince, R. A. DiStasio, A. Dreuw, B. D. Dunietz, T. R. Furlani, W. A. Goddard, S. Hammes-Schiffer, T. Head-Gordon, W. J. Hehre, C.-P. Hsu, T.-C. Jagau, Y. Jung, A. Klamt, J. Kong, D. S. Lambrecht, W. Liang, N. J. Mayhall, C. W. McCurdy, J. B. Neaton, C. Ochsenfeld, J. A. Parkhill, R. Peverati, V. A. Rassolov, Y. Shao, L. V. Slipchenko, T. Stauch, R. P. Steele, J. E. Subotnik, A. J. W. Thom, A. Tkatchenko, D. G. Truhlar, T. Van Voorhis, T. A. Wesolowski, K. B. Whaley, H. L. Woodcock, P. M. Zimmerman, S. Faraji, P. M. W. Gill, M. Head-Gordon, J. M. Herbert, and A. I. Krylov, J. Chem. Phys. 155, 084801 (2021). 50Q. Sun, J. Comput. Chem. 36, 1664 (2015). 51T. Kosaka and T. Wakabayashi, Heterocycles 41, 477 (1995). 52K. D. Nanda and A. I. Krylov, J. Chem. Phys. 142, 064118 (2015). 53M. de Wergifosse, E. Botek, E. De Meulenaere, K. Clays, and B. Champagne, J. Phys. Chem. B 122, 4993 (2018). 54M. Drobizhev, N. S. Makarov, S. E. Tillo, T. E. Hughes, and A. Rebane, Nat. Methods 8, 393 (2011). 55K. Aidas, C. Angeli, K. L. Bak, V. Bakken, R. Bast, L. Boman, O. Christiansen, R. Cimiraglia, S. Coriani, P. Dahle, E. K. Dalskov, U. Ekström, T. Enevoldsen, J. J. Eriksen, P. Ettenhuber, B. Fernández, L. Ferrighi, H. Fliegl, L. Frediani, K. Hald, A. Halkier, C. Hättig, H. Heiberg, T. Helgaker, A. C. Hennum, H. Hettema, E. Hjertenæs, S. Høst, I. Høyvik, M. F. Iozzi, B. Jansík, H. J. A. Jensen, D. Jonsson, P. Jørgensen, J. Kauczor, S. Kirpekar, T. Kjærgaard, W. Klopper, S. Knecht, R. Kobayashi, H. Koch, J. Kongsted, A. Krapp, K. Kristensen, A. Ligabue, O. B. Lut- næs, J. I. Melo, K. V. Mikkelsen, R. H. Myhre, C. Neiss, C. B. Nielsen, P. Norman, J. Olsen, J. M. H. Olsen, A. Osted, M. J. Packer, F. Pawlowski, T. B. Pedersen, P. F. Provasi, S. Reine, Z. Rinkevicius, T. A. Ruden, K. Ruud, V. V. Rybkin, P. Sałek, C. C. M. Samson, A. S. de Merás, T. Saue, S. P. A. Sauer, B. Schimmelpfen- nig, K. Sneskov, A. H. Steindal, K. O. Sylvester-Hvid, P. R. Taylor, A. M. Teale, E. I. Tellgren, D. P. Tew, A. J. Thorvaldsen, L. Thøgersen, O. Vahtras, M. A. Watson, D. J. D. Wilson, M. Ziolkowski, and H. Ågren, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 4, 269 (2014). 56Dalton, a molecular electronic structure program, release dalton2018.0 (2018), see http://daltonprogram.org. 57TURBOMOLE V7.6 2021, a development of university of karlsruhe and forschungszentrum karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007, available from https://www.turbomole.org. 58A. Willetts, J. E. Rice, D. M. Burland, and D. P. Shelton, J. Chem. Phys. 97, 7590 (1992). 59J. A. Pople, D. L. Beveridge, and P. A. Dobosh, J. Chem. Phys. 47, 2026 (1967). 60J. A. Pople, Trans. Faraday Soc. 49, 1375 (1953). 61D. Jacquemin, V. Wathelet, E. A. Perpète, and C. Adamo, J. Chem. Theory Comput. 5, 2420 (2009). 62D. J. Tozer, J. Chem. Phys. 119, 12697 (2003). 63A. J. Cohen, P. Mori-Sánchez, and W. Yang, Chem. Rev. 112, 289 (2012). 64L. Kronik, T. Stein, S. Refaely-Abramson, and R. Baer, J. Chem. Theory Comput. 8, 1515 (2012). 65S. Löffelsender, P. Beaujean, and M. de Wergifosse, “Simplified quantum chem- istry methods to evaluate non-linear optical properties of large systems,” WIREs Comput. Mol. Sci. 14, e1695 (2023). 66M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, and J. A. Montgomery, J. Comput. Chem. 14, 1347 (1993). 67M. S. Gordon and M. W. Schmidt, in Theory and Applications of Computational Chemistry, edited by C. E. Dykstra, G. Frenking, K. S. Kim, and G. E. Scuseria (Elsevier, Amsterdam, 2005), pp. 1167–1189. 68M. Chołuj, M. M. Alam, M. T. P. Beerepoot, S. P. Sitkiewicz, E. Matito, K. Ruud, and R. Zaleśny, J. Chem. Theory Comput. 18, 1046 (2022). 69M. Véril, A. Scemama, M. Caffarel, F. Lipparini, M. Boggio-Pasqua, D. Jacquemin, and P.-F. Loos, WIREs Comput. Mol. Sci. 11, e1517 (2021). 70P.-F. Loos, M. Comin, X. Blase, and D. Jacquemin, J. Chem. Theory Comput. 17, 3666 (2021). 71P.-F. Loos and D. Jacquemin, J. Phys. Chem. A 125, 10174 (2021). J. Chem. Phys. 160, 204110 (2024); doi: 10.1063/5.0206380 160, 204110-17 Published under an exclusive license by AIP Publishing 29 M ay 2024 06:59:37 https://pubs.aip.org/aip/jcp https://doi.org/10.1063/5.0055522 https://doi.org/10.1002/jcc.23981 https://doi.org/10.3987/com-94-6954 https://doi.org/10.1063/1.4907715 https://doi.org/10.1021/acs.jpcb.8b01430 https://doi.org/10.1038/nmeth.1596 https://doi.org/10.1038/nmeth.1596 https://doi.org/10.1002/wcms.1172 https://doi.org/10.1002/wcms.1172 http://daltonprogram.org https://www.turbomole.org https://doi.org/10.1063/1.463479 https://doi.org/10.1063/1.1712233 https://doi.org/10.1039/tf9534901375 https://doi.org/10.1021/ct900298e https://doi.org/10.1021/ct900298e https://doi.org/10.1063/1.1633756 https://doi.org/10.1021/cr200107z https://doi.org/10.1021/ct2009363 https://doi.org/10.1002/wcms.1695 https://doi.org/10.1002/wcms.1695 https://doi.org/10.1002/jcc.540141112 https://doi.org/10.1021/acs.jctc.1c01056 https://doi.org/10.1002/wcms.1517 https://doi.org/10.1021/acs.jctc.1c00226 https://doi.org/10.1021/acs.jpca.1c08524