A quantum chemical investigation of the second hyperpolarizability of p-nitroaniline

Because of their potential applications in the fields of optical communications, optical switching and computing, optical data storage,^1–4 and bioimaging (cells and tissues),^5–11 molecules and materials bearing large nonlinear optical (NLO) properties have gained an increased interest. Applications are based on the second-order⁵ and third-order⁶ NLO responses, but the materials for the latter ones have been less developed. At the molecular scale, the third-order NLO responses are described by the second hyperpolarizabilities, γ. To target applications and design new material exhibiting large γ values, both experimental and computational studies have been carried out since several decades. Hence, families of organic and organometallic compounds have been targeted, owing to the high-tunability of their structures using synthetic approaches. Experimentally, γ is measured using the electric field-induced second harmonic generation (EFISHG),^12–14 third harmonic generation (THG),¹⁵ and the third-harmonic scattering (THS)^16,17 techniques. On the other hand, computational-chemistry investigations, based on quantum chemistry (QC) calculations, are useful for rationalizing experimental data and for a detailed understanding of the structure-γ relationships. Hence, QC determinations of γ have served as guidelines for experimentalists to design new third-order NLO compounds.^18–20 

Nevertheless, predicting the second hyperpolarizabilities remains challenging for QC because several aspects must be finely described, including the electron correlation,^21–25 the frequency dispersion,²⁶ the surrounding,^27–29 and the vibrational contributions.^30–32 In addition, the choice of large basis sets including diffuse functions is generally essential, particularly for small molecules.^22,33,34 Many QC methods have been developed to evaluate and interpret the second hyperpolarizabilities from time-(in)dependent perturbation theory, leading to sum-over-states (SOS) expressions of γ,^35–37 to derivatives and response theory approaches.³⁸ 

Among the derivative approaches (in which the hyperpolarizabilities are expressed as the successive derivatives of the energies or of lower-order properties like the dipole moment and the linear polarizability), several methods proceed via analytical derivatives, which allows calculating the (nonlinear) optical responses with respect to dynamic electric fields.^39,40 These techniques have their equivalent in response theory.^38,41,42 Alternatively, the derivatives can be performed completely or in part by adopting the finite field (FF) method.⁴³ In the case of these numerical derivatives, only responses to static electric fields can be evaluated. The FF method has the advantage that it can easily be combined with a broad range of approximations, allowing the comparison of their performance. This has been enacted using Møller–Plesset perturbation theory (MPPT) and Coupled-Cluster (CC) schemes, in which the treatment of electron correlation can be systematically improved. However, these schemes are very computationally demanding and cannot be applied to large systems. In this context, density functional theory (DFT)^44,45 and time-dependent DFT (TDDFT)⁴⁰ methods, which can account for electron correlation effects at lower computational costs, constitute relevant alternatives. However, their accuracy depends on the approximations used in the exchange–correlation functionals (XCFs), and several studies have highlighted the difficulty of defining an appropriate XCF for calculating the molecular responses to external electric fields,^46–52 since these are inherently nonlocal. Hence, conventional XCFs within the local density approximation (LDA) or the generalized gradient approximation (GGA) have been reported to drastically fail for computing the NLO properties. Due to their short-sightedness, they lead to over-delocalization of the responses to external fields, which has been related to the self-interaction error and the integer discontinuity of the XC potential.^47,49,50 Several solutions have been proposed, including the introduction of exact Hartree–Fock (HF) exchange within global or range-separated hybrid (RSH) XCFs. In the latter case, the system-specific tuning of the range separation parameter was shown to minimize the delocalization error (DE),⁴⁹ and thus, it can help in improving the accuracy of the calculated hyperpolarizabilities.^25,53,54

Recent experimental measurements of γ using THS^16,17 have revived interest in calculating molecular second hyperpolarizabilities. Some of us recently addressed small molecules,^34,55,56 and in this work, we move toward more complex systems by studying p-nitroaniline (pNA) (Fig. 1) as a prototypical push-pull π-conjugated molecule. The first hyperpolarizability of pNA has been calculated using a broad range of methods,^57–66 but the second hyperpolarizability has been much less studied.^67–69 The current study focuses on the effects of electron correlation on γ by using wavefunction methods (MPPT and CC) and addresses the performance of a broad range of XCFs. This work also investigates numerical aspects associated with the FF procedure and with the basis set convergence as well as the role of frequency dispersion in the third-order NLO responses of pNA, namely, the optical Kerr effect (OKE), EFISHG, and THG.

This paper is organized as follows: Sec. II defines the NLO properties, summarizes the numerical expressions in the context of the FF approach, and presents the theoretical methods that have been used. Then, Sec. III gives the computational details and the optimal parameters for numerical differentiations. The results are presented and discussed in Sec. IV before conclusions are drawn in Sec. V.

The above equalities are valid for exact wavefunctions. For approximate wavefunctions, when the induced dipole moment results from the minimization, with respect to all variational coefficients, of the total energy of the system interacting with the electric field, the responses obtained from the energy derivatives are also equivalent to those obtained from derivatives of the properties.

These derivatives can be evaluated numerically or analytically or by combining both, i.e., by numerically differentiating lower-order responses obtained within an analytical procedure. The analytical schemes encompass linear, quadratic, and cubic response theories at various levels of approximation (ranging from HF to DFT and CC schemes)^38–40 or sum-over-states methods.^70–72 These analytical approaches allow the calculation of both the static and the frequency-dependent responses. They are usually computationally more demanding (with difficulties when applied to large systems), and they are available for a narrower range of approximations. This is particularly the case for the second hyperpolarizability, which explains the success of (partial) numerical derivative approaches.

In this section, we briefly introduce the computational schemes employed for evaluating molecular second hyperpolarizabilities. We first present the quantum-chemical methods used to determine field-dependent energies, from wavefunction-based to DFT methods. Next, we briefly present response theory approaches, which can be combined with wavefunction or DFT approaches to calculate static and dynamic molecular properties.

The Hartree–Fock method defines a kind of zero-order method, with its well-known drawback of neglecting electron correlation. Yet, this approximate mean-field approach is the starting point of higher-order methods. Electron correlation can be accounted for by using Møller–Plesset (MP) perturbation theory. In this work, it was applied at second (MP2), third (MP3), and fourth (MP4) orders, which are standard levels available in many computational packages. The second- and third-order corrections originate only from doubly excited configuration state functions (doubles), while at fourth-order, singles, doubles, triples, and quadruples contribute. This allows us to address the different contributions to the MP4 correlation energy, among these, the MP4D method (only the contributions from the doubles are included), the MP4DQ method (the singles and triples are not included), and the MP4SDQ (only the triples are neglected). The MPn schemes are non-variational and size-consistent at any order (n).

An efficient way to improve the convergence of the electron correlation as a function of the nature of the excitations consists in employing the exponential coupled cluster (CC) ansatz on the reference HF wavefunction (Ψ^HF), $ΨCC=eT̂ΨHF$⁠, where $T̂$ is the cluster operator defined as the sum of single, double, and higher-order excitations. To get a hierarchy of methods, the natural choice is to truncate the cluster expansion to the doubles (CCD), then to the singles and doubles (CCSD), to the Coupled-Cluster Singles and Doubles with perturbative treatment of the Triples (CCSDT), and so on. Here, the CCSD scheme was adopted as well as the gold standard CCSD(T) method, where the $T̂$ expansion is truncated after the singles and doubles, while the triples are treated, a posteriori, in a perturbative manner. CCSD(T) is considered here as the reference method, used to assess the performance of the other approaches.^25,77,78

An efficient, computationally advantageous, alternative to these highly correlated wavefunction methods consists in adopting Kohn–Sham DFT, where the central quantity is the electron density $ρr$ so that the energy of the system is described as a functional of $ρr$⁠.^44,45 The exact XC functional is, however, unknown. To work out better and better approximations, the field of DFT has been very active over the last 40 years with the developments of families of improved XCFs. These have been obtained either with the strategy to fulfill mathematical and physical rules and/or to reproduce broader and broader set of data via their improved parameterizations. Yet, many approximate XCFs are not satisfactory to describe the linear and nonlinear optical responses of molecules due to their long-range nature. This has been evidenced for LDA XCFs but also for their evolutions using the gradient of the density (GGA) or its kinetic energy density (meta-generalized gradient approximation, mGGA), which lead to responses that can be orders of magnitude larger than reference values, especially for large π-conjugated systems.^46,47 This incorrect behavior has been attributed to the short-sightedness of these functionals or, in other words, to the incorrect field dependence of the response part of the exchange functional, which lacks a linear term counteracting the external electric field.⁷⁹ These drawbacks have been related to the self-interaction error or over-delocalization of the response to external electric fields.⁶⁵ 

An alternative to improve the XCF consists in replacing a fraction of the DFT correlation by second-order perturbation theory (PT2) contribution, evaluated from a preliminary DFT calculation, leading to double hybrid XCFs.^90,91 Their reliability to evaluate the first hyperpolarizability of push-pull π-conjugated systems was demonstrated in Ref. 48, though MP2 remains superior. On the other hand, much less is known about their performance for computing the second hyperpolarizabilities.

They have been derived in a general frame, where the time evolution of the expectation value of a Hermitian operator (⟨A⟩) is described for a system undergoing a time-dependent perturbation, such as an oscillating electric field.⁹⁷ The expansion terms are the successive linear, quadratic, cubic, … response functions, which can be static/dynamic.⁹⁸ Hence, the first-order dipole moment at pulsation ω induced by an electric field at the same pulsation defines the frequency-dependent polarizability, $−⟨⟨μi;μj⟩⟩ω=αij−ω;ω$⁠. Similarly, when applying electric fields at pulsations ω₁ and ω₂, the second-order induced dipole, oscillating at the sum frequency ω_σ = ω₁ + ω₂, defines the dynamic first hyperpolarizability, $−⟨⟨μi;μj,μk⟩⟩ω1,ω2=βijk−ωσ;ω1,ω2$ and so on for $−⟨⟨μi;μj,μk,μl⟩⟩ω1,ω2,ω3=γijkl−ωσ;ω1,ω2,ω3$⁠, where $ωσ=∑i1,2,…ωi$ is the frequency of the generated light. These response functions can be evaluated at different levels of approximation.³⁸ At the HF level, it leads to the time-dependent Hartree–Fock (TDHF)³⁹ scheme, while its DFT counterpart is known as the time-dependent DFT (TDDFT)⁴⁰ approach. In the static limit, one recovers the coupled-perturbed HF (CPHF) and coupled-perturbed Kohn–Sham (CPKS) scheme, respectively. Since the HF method is variational, the analytical CPHF scheme provides the same molecular properties as the FF method employing HF energies. The same equivalence holds for CPKS and FF/KS, provided there is no approximation in the XC kernels. Besides TDHF and TDDFT, the response theory approaches have also been worked out at higher levels of approximation, including CC levels. This has defined a hierarchy of approximations: CCS, CC2, CCSD, CC3, ….^99,100

The geometry of pNA was fully optimized at the MP2/6-311G(d) level in gas phase, with thresholds on the forces and displacements of 1.5 × 10⁻⁵ hartree bohr⁻¹ and 6.0 × 10⁻⁵ Å, respectively. All harmonic vibrational frequencies calculated for the optimized structure are real, confirming that it corresponds to a minimum on the potential energy surface. The nitro-amino charge-transfer (CT) direction defines the x-axis, while the phenyl ring determines the xy-plane (Fig. 1). The molecular energies were obtained using the SCF procedure (HF and DFT) with a threshold of 10⁻¹¹ a.u. on the energy. The same convergence threshold was applied in the CCSD iterative procedure. Field-dependent energies were calculated with a broad set of wavefunction methods, namely, HF, MP2, MP3, MP4D, MP4DQ, MP4SDQ, CCSD, and CCSD(T). In addition, field-dependent polarizabilities were calculated at the HF and MP2 level. DFT was used to evaluate the field-dependent energies as well as the field-dependent polarizabilities and first hyperpolarizabilities. The selected XCFs are given in Table I together with their key characteristics. Static responses have been evaluated using the CPKS scheme, while dynamic responses [$α−ω;ω,β−ω;ω,0,andβ−2ω;ω,ω$] were calculated at the TDDFT level. The photon energies (ℏω) range from 0.0 to 0.1 a.u., corresponding to ℏω = 2.72 eV (λ = 456 nm). In those CPHF/TDHF and CPKS/TDDFT calculations, the convergence threshold on the response of the density matrix was set to 10⁻¹⁰ a.u. These calculations were carried out using the Gaussian16 program.¹⁰³ Moreover, cubic response functions (CRFs) have been calculated at the CCSD level using the Dalton program¹⁰⁴ using the unrelaxed field-independent orbitals.

The analyses concentrate on several γ quantities, i.e., (i) γ_xxxx, the diagonal component along the CT axis, (ii) γ_//, which is probed by EFISHG,^12–14 and (iii) the $γTHS=⟨γZZZZ2⟩+⟨γZXXX2⟩1/2$ and $DR=⟨γZZZZ2⟩/⟨γZXXX2⟩$ quantities of THS.^16,17,34 All second hyperpolarizabilities are reported in atomic units (1 a.u. of γ = 6.235 377 × 10⁻⁶⁵ C⁴ m⁴ J⁻³ = 5.0367 × 10⁻⁴⁰ esu), within the T convention, as defined in Eq. (1).¹⁰⁵ 

Owing to the high numerical precision that is required to perform the numerical differentiations (see Subsection II B), a key computational parameter in DFT calculations is the density of the integration grid. These numerical integrations are performed atom by atom, in spherical coordinates centered at the nuclei, with the three-dimensional integration decomposed into a one-dimensional radial and a two-dimensional angular integration.^106–108 The default grid in Gaussian16 calculations is called ultrafine (99 radial shells and 590 angular points per shell), but it is of lower quality than the superfine grid (175 radial shells and 974 angular points per shell for H atoms and 250 radial shells and 974 angular points per shell for the other atoms). To gauge the impact of the choice of integration grid, static γ quantities were calculated using both grids and a selection of XCFs (Table S1). Both grids provide similar γ_xxxx, γ_//, and γ_THS values with differences of a few tens of a.u., which is small in comparison to the amplitude of the whole responses (10⁴–10⁵ a.u.). Though the differences are larger for CAM-B3LYP and LC-BLYP in comparison to B3LYP, M06, M06-HF, and M06-2X, these remain small. The effect on depolarization ratio (DR) is negligible. Note that the computational needs do not increase much when using the denser grid, with central processing unit (CPU) time on average 40%–45% larger. Therefore, the superfine integration grid was used for all DFT/TDDFT calculations.

The basis set effects were investigated by considering the static γ responses, as evaluated with the hybrid numerical derivative scheme where field-dependent polarizability values are calculated at the CPKS level with the CAM-B3LYP XCF (Table II and Fig. 2). The list of basis set encompasses Pople’s basis sets¹⁰⁹ [6-31G(d), 6-31(d,p), 6-311G(d), 6-311G(d,p), 6-31+G(d), 6-31+G(d,p), 6-311+G(d), and 6-311+G(d,p)] and Dunning’s correlation consistent basis sets¹¹⁰ [cc-pVXZ (X = D, T, Q), aug-cc-pVXZ (X = D, T, Q), and d-aug-cc-pVXZ (X = D, T, Q)]. The results obtained with the augmented or doubly augmented Dunning’s basis sets show a nice convergence with the number of valence functions, which allows concluding that (i) in the complete basis set limit, γ_xxxx ∼ 110 × 10³ a.u., γ_// ∼ 32 × 10³ a.u., γ_THS ∼ 43 × 10³ a.u., and DR ∼ 11.0 and that (ii) the d-aug-cc-pVQZ basis set can be considered as reference. As observed for other systems,^33,111 the lack of diffuse functions in both Pople and Dunning basis set families is detrimental to accurate second hyperpolarizability values and leads to large underestimations even if the cc-pVXZ series displays the good trend. Comparing cc-pVDZ to d-aug-cc-pVDZ, γ_xxxx is underestimated by 37%, while γ_THS is underestimated by 44% and γ_// by 56%, highlighting that the component along the CT axis is less impacted than the perpendicular and out-of-plane ones. Among Pople’s basis sets, the 6-311+G(d) basis set performs well with underestimations of γ_xxxx, γ_THS, and γ_// by 1.5%, 9.4%, and 5.3%, respectively. On the other hand, the number of basis functions [238] is more than five times less than in d-aug-cc-pVQZ [1426]. Adding p polarization functions in Pople basis sets has a marginal impact. In both Pople’s or Dunning’s series, the effects of adding diffuse functions are similar for double-ζ, triple-ζ, and quadruple-ζ basis sets. These basis set effects were confirmed at the MP2 level for a selected case [6-311+G(d) vs d-aug-cc-pVDZ] (Table S2). From these calculations, the 6-311+G(d) basis set was identified as a good compromise between cost and accuracy and used in the rest of this study.

As discussed above, both pure numerical and hybrid differentiation schemes were used to evaluate the second hyperpolarizabilities. The amplitudes of the electric field were generated using a geometric progression $Ek,n=2knE0$⁠, with E⁰ = 0.0003 a.u., n = 1, and k = 0, 1, 2, …. In order to obtain all tensor components, these fields were applied in each direction of the molecular frame. The Romberg iteration formula and the subsequent selection of the best numerical derivative estimate and its error bar were enacted with the NACHOShomemade code.^56,112

To determine the convergence region in Romberg’s tables, electric field amplitudes from 0.0003 to 0.0192 a.u. were employed for a selection of calculations. A typical table is provided in Table S3 for the CAM-B3LYP/6-311+G(d) calculation of γ_xxxx, highlighting that the (k, m) = (2, 2) value is in a stability domain, proposing that γ_xxxx ∼ (10 858 ± 2) × 10¹ a.u. The same approach has been followed when adopting the first-order derivative of analytically evaluated β_xxx or second-order derivative of analytically evaluated α_xx (Table S4). As expected, it turns out that for α and β, the stability domain is broader and that γ_xxxx can be estimated with a higher precision of 1.0 a.u. (γ_xxxx = 108 582 ± 1 a.u.). These comparisons have been extended to all the γ quantities defined above and to different levels of approximation (Table S5).

These results confirm that when analytical derivatives are available to calculate lower-order properties, they should be used for a matter of precision and efficiency. Furthermore, they also provide the precision of the numerical estimates obtained from fourth-order derivatives. For γ, it amounts to about 10 a.u., which is largely sufficient for the purpose of this work.

Using the 6-311+G(d) basis set and the optimal parameters for the calculation of numerical derivatives, the static γ values were evaluated using wavefunction methods and compared to reference CCSD(T) results (Table III). All methods underestimate the γ values, and the largest underestimation is found with the HF method (50%), highlighting the crucial role played by electron correlation. The MP2 approximation is a good compromise between accuracy and computational cost, with underestimations of the order of 10%. It performs better than MP3 and also better than the incomplete MP4 levels. Moreover, the MP2 error is reduced by a factor of 3 upon including the third- and whole fourth-order corrections, demonstrating that (i) the performance of the MP2 scheme results from errors compensations, i.e., from partial cancellations between higher-order terms (singles, doubles, and triples bring positive contributions, which is partially counterbalanced by the quadrupoles) and that (ii) the contribution of the triples at fourth-order is important and amounts to 18%. This latter contribution is also responsible for the 15% underestimation of CCSD(FF) and the fact that it is less accurate than MP4 (and also than MP2). Using the CRF approach with unrelaxed MOs, the CCSD values are very close to the reference CCSD(T) results. This is another evidence of cancellation between two contributions: the inclusion of perturbative triples and the relaxation of MOs. These results on pNA are in line with other theoretical studies where it was found that in the case of the static electrical properties, the inclusion of orbital relaxation does not necessarily provide superior results.^66,113–115 In addition, the effects of the methods on DR are smaller, with an overestimation by 10% at the HF level and small underestimations with the correlated methods.

In this section, we address the performance of various DFT XCFs, ranging from LDA to double hybrids (Table IV). The results were obtained with the 6-311+G(d) basis set and the superfine integration grid. Except for CCSD(T), where γ’s are evaluated as 4th-order derivatives of the energy, the values have been obtained from the hybrid numerical differentiation scheme i.e., from the second-order derivatives of the polarizability. As observed with the wavefunction methods (though there is no obvious link), all XCFs underestimate the second hyperpolarizabilities. Apart from the mPW2PLYP and B2-PLYP double-hybrid functionals, the errors are larger than 10%, by 15%–39% for γ_//, and by 22%–41% for γ_THS. These results show that calculating γ of pNA is challenging for all types of functionals, whether conventional, global hybrids, or RSHs.

Isolating systematic trends between the nature of the XCF and the γ values is not straightforward. In the Minnesota series that includes M06-L (0% of HF exchange), M06 (27%), M06-2X (54%), and M06-HF (100%), γ first increases up to M06 and then decreases. This highlights to a given point the role played by the amount of HF exchange, usually associated with a decrease in the linear and nonlinear responses. Figure 3 illustrates this trend for a larger set of XCFs. This exact exchange effect is also observed when decreasing the μ value in RSHs, which corresponds to lowering the amount of HF exchange at short-range. This is illustrated by the increase in the γ values from ωB97 (μ = 0.4 bohr⁻¹), to ωB97X (μ = 0.3 bohr⁻¹), and ωB97X-D (μ = 0.2 bohr⁻¹), although such a comparison is partially biased since all parameters, not only μ, differ in this series of functionals. Globally, the γ values obtained with RSHs using their default μ parameters do not differ much from those obtained with the corresponding global hybrids. The impact of tuning μ is investigated in more detail in Sec. IV D.

The depolarization ratio (DR) of the THS response ranges from 9.0 to 12.4 (Table IV). LDA and GGA XCFs clearly perform poorly, in comparison to the hybrids and double hybrids. Obviously, the accuracy of the values raises with the amount of exact exchange, with the best results obtained with some of the range-separated hybrids (M11, LC-ωPBE, and LC-BLYP) and double hybrids.

Assuming Kleinman’s symmetry conditions (which is the case for static values), the γ tensor can be decomposed into its isotropic (γ_J=0), quadrupolar (γ_J=2), and hexadecapolar (γ_J=4) spherical tensor components.¹¹⁶ As defined in Refs. 17 and 34, the anisotropy parameters $ρ0/2=γJ=0/γJ=2$ and $ρ4/2=γJ=4/γJ=2$ were evaluated (Tables V and S6). The reference CCSD(T) calculations predict that the quadrupolar term is slightly dominant and that the isotropic contribution is larger than the hexadecapolar one. These trends are generally well reproduced with global hybrids, range-separated hybrids, and double hybrids. Again, the B2-PLYP double-hybrid functional provides accurate estimates of these components compared to the CCSD(T) values.

Among the XCFs, the double-hybrid functionals provide the most accurate results, especially B2-PLYP and mPW2PLYP. Nevertheless, as illustrated in Fig. 4, their performance depends on the amount of PT2 correlation and HF exchange. On the one hand, while keeping the % of HF exchange at its default 53% value, for 0% of PT2 correlation, γ_THS amounts to 37 000 a.u. and it increases linearly with the amount of PT2 correlation (the sum of PT2 and DFT correlations is equal to 1) to reach a value of 111 000 a.u. for 100% of PT2 correlation (Table S7). The best results are obtained using PT2 percentages of 27%–30%. On the other hand, when increasing the amount of HF exchange but keeping the % of PT2 correlation to its 27% default value, γ_THS decreases (Fig. 4), consistently with the results of Fig. 3. Contrary to the PT2 case, the dependence of γ_THS on the percentage of HF exchange is no longer linear since the HF exchange term intervenes directly in the orbital optimization i.e., in the SCF procedure while the PT2 contribution is an a posteriori correction. γ_THS vs % (HF) curve crosses the CCSD(T) line at 53%, which is exactly the default value in the B2-PLYP functional. Among all the functionals, the original B2-PLYP that contains 27% of PT2 correlation and 53% of HF exchange shows the best accuracy and performs also better than MP2. B2-PLYP is thus recommended for the prediction of the static second hyperpolarizabilities of the pNA molecule.

As shown in the literature, tuning the range-separating parameter (μ) offers the possibility to better describe the molecular responses to electric fields.^25,89 Using LC-BLYP, CAM-B3LYP and ωB97 XCFs, μ was varied to monitor the $ΔIPμ$ function [Eq. (13)] and its components [Fig. 5(a) and Table S8]. $ΔIPμ$ decreases monotonically from μ = 0.05 bohr⁻¹ until it cancels at 0.28 bohr⁻¹ (0.26 bohr⁻¹) for LC-BLYP (ωB97), which corresponds therefore to the value where Koopmans’ theorem is satisfied. Then, $ΔIPμ$ gets more and more negative. For CAM-B3LYP, $ΔIPμ$ decreases slowly but remains positive for all μ values and with an asymptotic value of about 0.21 eV (Fig. S1). For that range of μ values, γ was calculated [Fig. 5(b) and Table S9]. With both LC-BLYP and ωB97, γ increases with μ until reaching a maximum at μ = 0.12 bohr⁻¹ and then decreases almost linearly. As already discussed above, the decrease in γ with μ is expected as it corresponds to increasing the amount of HF exchange at short-range (the amount of HF exchange at long-range always tends to 100% for LC-BLYP and ωB97 and to 65% for CAM-B3LYP). The presence of a maximum is, however, less intuitive. For the CAM-B3LYP functional, the γ value decreases monotonically with μ. Note that whatever the value of μ, γ remains underestimated with respect to CCSD(T). At μ = 0.12 bohr⁻¹, the underestimation amounts to 19% (24%) for γ_THS using LC-BLYP (ωB97) [and to 11% (17%) for γ_//], while for the “optimal” μ value, γ_THS is underestimated by about 25% (29%) [and 21% (24%) for γ_//]. Compared to results obtained with the default μ value, the error is slightly reduced. The μ value could also be chosen from the tuning scheme proposed by Besalú-Sala et al.⁸³ dedicated specifically to the calculation of second hyperpolarizabilities. Owing to the value of the LC-BLYP/6-311+G* longitudinal polarizability of pNA, α_xx = 135.15 a.u.., the size extensive descriptor of Ref. 83, amounts to 0.274, which corresponds to μ_Tα = 0.30 bohr⁻¹. If considering the average polarizability of pNA instead of its longitudinal component, the μ_Tα value becomes 0.33 bohr⁻¹. Both μ_Tα values are thus larger than the optimal one obtained using Eq. (13) and, therefore, do not improve the γ value of pNA. This is probably because no π-conjugated push-pull systems were included in the molecular dataset considered in Ref. 83.

The performances of XCFs to calculate γ, as presented in Subsections IV B and IV D, are now discussed in relation with their (de)localization errors. Figure 6 and Figs. S2–S4 plot the $ΔEδ$ quantity as a function of δ as well as its deviations with respect to the exact linear behavior described by $ΔEδ$ [Eq. (14)]. Quantitative estimates of these deviations are given by their curvatures in both the −1 ≤ δ ≤ 0 and 0 ≤ δ ≤ 1 intervals (see the values provided in the parenthesis in Fig. 6). Positive and negative curvatures (corresponding to negative and positive ΔΔɛ values) are associated with delocalization and localization errors, respectively.^93,117

Although none of the functionals display a linear behavior of $ΔEδ$ in any of the intervals, the deviations from linearity vary considerably between the XCFs. The HF method displays the largest localization error and then come XCFs with a large amount of HF exchange: M06-HF, LC-ωPBE, ωB97, ωB97X, and LC-BLYP (μ = 0.47 bohr⁻¹ and μ = 0.40 bohr⁻¹). Those are the XCFs that most underestimate the γ responses. The next functionals, in order of decreasing localization, are ωB97X-D, M11, and LC-BLYP (μ = 0.28 bohr⁻¹) with negative or positive curvatures of small amplitudes (smaller than 0.2 eV). All other XCFs present a delocalization error (DE). DE is the largest for functionals without HF exchange [SVWN, BLYP, Perdew–Burke–Ernzerhof (PBE), B97-D, and M06-L] or with small amounts of HF exchange (B3LYP, PBE0, M06, and MN15). CAM-B3LYP with “only” 65% HF exchange in the long-range and M06-2X with 54% HF exchange have also large delocalization errors, as evidenced by the curvature amplitudes. Among the Minnesota XCFs, M11 is behaving the best. LC-BLYP with the optimal μ = 0.28 bohr⁻¹ displays the smallest curvatures, while LC-BLYP with smaller or larger μ values perform more poorly. However, small DE values (i.e., small curvature coefficients) are not associated with more accurate γ values (Figs. 7 and S5 and Table S10). On the other hand, the good behavior of global hybrids, such as B3LYP and PBE0, for predicting the γ values over the tuned LC-BLYP and other XC functionals results in the compensation of their large DE by the dynamic electron correlation as it has been reported for other π-conjugated systems.^25,54 This is confirmed by the B2-PLYP results where 27% of dynamic PT2 correlation perfectly compensates the DE and, thus, provides the most accurate second hyperpolarizabilities. Including more than 27% of PT2 correlation unbalances the compensation, and thus, the second hyperpolarizabilities increase drastically.

Before concluding on the performance of wavefunctions and DFT methods to evaluate the second hyperpolarizability of pNA, it is interesting (i) to relate the performance of a given level of approximation with the required CPU times and (ii) to assess, on an equal foot, how these methods perform for estimating its first hyperpolarizability. On the one hand, Tables S12 and S13 list indicative CPU times for calculating the static γ tensor from using a single field amplitude [Eqs. (7) and (8) with E_i = 0.0012 a.u.] with wavefunction and DFT methods. Energy calculations, even with a tight threshold, are fast, of the order of the minute of CPU time. The CPU time is multiplied by a factor of about two between HF and MP2 and again by about two between MP2 and MP3. At fourth order, calculating the contributions from the triples becomes substantially more demanding (about 2 h of CPU time). Each CCSD energy calculation amounts also to about two hours so that a CCSD(T) single point calculation takes about twice longer. When evaluating the γ tensor from field derivatives of the polarizability, in addition to the SCF cycle, three CPHF/CPKS cycles need to be carried out (one for each Cartesian direction of the molecular frame). At the HF level, the CPU time is about twice larger than what is needed for the pure numerical differentiation scheme. This larger value is, however, compensated by the fact that the hybrid numerical differentiation scheme requires less field amplitudes (smaller number of k values) than the pure numerical differentiation scheme to achieve convergence in the Romberg table. At the MP2 level, the balance is, however, clearly in favor of the scheme based on energy derivatives. At the DFT level, each polarizability calculation takes, on average, between 21 and 104 min, with the longest CPU times for double hybrids and the smallest for LDA. Besides double hybrids, the advantage of DFT with respect to MP2 in terms of CPU is small, which is, in part, attributed to the use of superfine integration grids in the former.

On the other hand, since the β response is a lower-order property than γ, the β values are provided at no extra cost. Among the wavefunction methods, the MP2 and MP4 levels reproduce closely the CCSD(T) results on β_// and β_SHS (Table S14). Yet, like for the γ quantities, (i) the other levels of approximation underestimate the β quantities, (ii) the accuracy of the MP2 values originates from error cancelations between the doubles and quadruples on one side and the singles and triples on the other side, (iii) unrelaxed quadratic response function (QRF) CCSD values are more accurate than the relaxed CCSD results (from FF calculations), and (iv) globally, the errors with respect to CCSD(T) are smaller for the β quantities than for the γ quantities. At the DFT level, the results are contrasted with respect to those observed for the γ quantities (Table S15): (i) in the absence of HF exchange, the β_// and β_SHS values are overestimated (11%–25%, from M06-L to BLYP), (ii) adding more and more HF exchange decreases the β values so that, for instance, M06 overestimates β_SHS by 5%, whereas M06-2X underestimates it by 13%, (iii) similarly, the β values decrease from BLYP (overestimation of β_SHS by 25%) to B3LYP (overestimation of β_SHS by 12%) and to LC-BLYP (underestimation of β_SHS by 19%), (iv) the β values also increase when decreasing the range-separating parameter (from ωB97 to ωB97X-D), and (v) the double hybrids perform well but not necessarily better than some global or range-separated hybrids. In summary, the PBE0, M06, MN15, CAM-B3LYP, PBE0DH, and MPW2PLYP XCFs give β quantities within 10% of the reference CCSD(T) value, in comparison to the fact that all XCFs underestimate the γ quantities. Additional calculations performed with the modified B2-PLYP double hybrid XCF show the variations in the β quantities as a function of the percentages of HF exchange or PT2 correlation while keeping the other contribution fixed at its default value (Table S16 and Fig. S6). Here, the closest agreement with CCSD(T) occurs for smaller amount of PT2 correlation (15%) but larger percentage of HF exchange (63%), again showing the interplay between both contributions. Then, with RSHs and their default μ values, β values are underestimated, but decreasing μ (decreasing the amount of HF exchange at short range) leads to an increase of β, and there exists an optimal μ value for reproducing the CCSD(T) results (Table S17 and Fig. S7). This optimal μ amounts to 0.19 bohr⁻¹ for LC-BLYP and CAM-B3LYP and to 0.14 bohr⁻¹ for ωB97, which differ from the μ values that satisfy $ΔIPμ=0$⁠. At those best μ values for β, the γ quantities are underestimated by about 20% (Table S9).

In summary, owing to the revival of interest on the second hyperpolarizability, γ, associated with recent measurements of the third harmonic scattering responses,^16,17 this paper has addressed the performance of wavefunction and DFT methods to evaluate γ of the pNA prototypical push-pull π-conjugated molecule. The performance of these methods has been assessed in comparison to the CCSD(T) reference results. Among wavefunction-based methods, MP2 offers the best accuracy/cost ratio for computing the static γ, although its good performance results from compensation of errors. At the DFT level, γ remains challenging to compute. Conventional, global hybrid, or RSH XCFs underestimate static γ values by at least 15%. Then, γ is very sensitive to the dynamic correlation, as estimated by tuning the PT2 contribution of double hybrids. Among all functionals, the original double-hybrid B2-PLYP functional, which benefits from the 27% of PT2 correlation and 53% HF exchange, provides accurate estimates of static γ values. Unfortunately, the best performing methods for γ are not necessarily reliable for the β quantities, except MP2. In fact, the β quantities of pNA could be predicted, with a good accuracy, with a range of hybrid XCFs, but these systematically underestimate γ.

Finally, frequency-dependent OKE and EFISHG second hyperpolarizabilities have been calculated with the CAM-B3LYP XCF for different incident photon wavelengths. Polynomial expressions linking the static and frequency-dependent hyperpolarizabilities have been deduced and allowed to describe frequency dispersion effects occurring in NLO processes implying only dynamic incident electric fields, such as THG.

Further work will consider other π-conjugated organic molecules searching for QC methods or DFT XCFs that can be generalized to the second hyperpolarizabilities calculation and that could be used in parallel with experimental measurements to address structure–NLO property relationships. These investigations will also aim at assessing how the DR is impacted by the molecular structure, the size of the π-conjugated backbone, and the number and nature of the donor/acceptor substituents.

See the supplementary material for (i) the effect of the integration grid and the basis set on the second hyperpolarizability of pNA, (ii) Romberg tables, (iii) comparisons between first, second, and fourth-order numerical derivative schemes, (iv) the decomposition of the second hyperpolarizability into its isotropic invariants, (v) the variations of the second hyperpolarizability as a function of the percentages of HF exchange and of PT2 correlation in modified B2-PLYP XCF, (vi) the relationships between the (de)localization error, the second hyperpolarizability, and the range-separating parameter, (vii) frequency dispersion of the second hyperpolarizability, and (viii) CPU times for calculating the tensor with wavefunction and DFT methods as well as (ix) the selected results on the related first hyperpolarizability of pNA.

The calculations were performed on the computers of the Consortium des équipements de Calcul Intensif (CÉCI) (http://www.ceci-hpc.be), including those of the UNamur Technological Platform of High-Performance Computing (PTCI) (http://www.ptci.unamur.be) and those of the Tier-1 supercomputer of the Fédération Wallonie-Bruxelles, for which we gratefully acknowledge the financial support from the FNRS-FRFC, the Walloon Region, and the University of Namur (Convention Nos. 2.5020.11, GEQ U.G006.15, U.G018.19, U.G011.22, RW1610468, RW/GEQ2016, 1117545, and RW2110213).

The authors have no conflicts to disclose.

Komlanvi Sèvi Kaka: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Pierre Beaujean: Formal analysis (supporting); Methodology (supporting); Validation (equal); Writing – review & editing (equal). Frédéric Castet: Formal analysis (supporting); Validation (equal); Writing – review & editing (equal). Benoît Champagne: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.