Recent measurements of the third harmonic scattering responses of molecules have given a new impetus for computing molecular second hyperpolarizabilities (*γ*) and for deducing structure–property relationships. This paper has employed a variety of wavefunction and density functional theory methods to evaluate the second hyperpolarizability of the *p*-nitroaniline prototypical push-pull *π*-conjugated molecule, addressing also numerical aspects, such as the selection of an integration grid and the impact of the order of differentiation vs the achievable accuracy by using the Romberg quadrature. The reliability of the different methods has been assessed by comparison to reference Coupled-Cluster Singles and Doubles with perturbative treatment of the Triples results. On the one hand, among wavefunction methods, the MP2 scheme offers the best accuracy/cost ratio for computing the static *γ*. On the other hand, using density functional theory, *γ* remains a challenging property to compute because all conventional, global hybrid or range-separated hybrid, exchange–correlation functionals underestimate static *γ* values by at least 15%. Even tuning the range-separating parameter to minimize the delocalization errors does not enable to improve the *γ* values. Nevertheless, the original double-hybrid B2-PLYP functional, which benefits from 27% of PT2 correlation and 53% Hartree–Fock exchange, provides accurate estimates of static *γ* values. Unfortunately, the best performing exchange–correlation functionals for *γ* are not necessarily reliable for the first hyperpolarizability, *β*, and vice versa. In fact, the *β* of *p*-nitroaniline (*p*NA) could be predicted, with a good accuracy, with several hybrid exchange–correlation functionals (including by tuning the range-separating parameter), but these systematically underestimate *γ*. As for *γ*, the MP2 wavefunction method remains the best compromise to evaluate the first hyperpolarizability of *p*NA at low computational cost.

## I. INTRODUCTION

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^{5} and third-order^{6} 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),^{15} 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,^{26} 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.^{38}

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.^{43} 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)^{40} 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),^{49} 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 (*p*NA) (Fig. 1) as a prototypical push-pull *π*-conjugated molecule. The first hyperpolarizability of *p*NA 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 *p*NA, 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.

## II. THEORY AND METHODOLOGY

### A. Definitions of the static polarizabilities and hyperpolarizabilities

*i*,

*j*,

*k*, … directions are described by the successive contributions to the induced dipole moment. For the

*i*th Cartesian component of

*μ*, it reads

*α*is the first-order molecular response or molecular polarizability, while

*β*and

*γ*are the second- and third-order responses, the first and second hyperpolarizabilities, respectively. These successive quantities are tensors of rank 1, 2, 3, …, which is also evidenced by their number of indices. In that event, the total electronic energy can be written as a function of the electric fields and the field-dependent dipole moment,

- Polarizability:(3)$\alpha ij=\u2212\u22022E\u2202Ei\u2202EjE\u20d7\u21920=\u2202\mu i\u2202EjE\u20d7\u21920.$
- First hyperpolarizability:(4)$\beta ijk=\u2212\u22023E\u2202Ei\u2202Ej\u2202EkE\u20d7\u21920=\u22022\mu i\u2202Ej\u2202EkE\u20d7\u21920=\u2202\alpha ij\u2202EkE\u20d7\u21920.$
- Second hyperpolarizability:(5)$\gamma ijkl=\u2212\u22024E\u2202Ei\u2202Ej\u2202Ek\u2202ElE\u20d7\u21920=\u22023\mu i\u2202Ej\u2202Ek\u2202ElE\u20d7\u21920=\u22022\alpha ij\u2202Ek\u2202ElE\u20d7\u21920=\u2202\beta ijk\u2202ElE\u20d7\u21920.$

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.

### B. Second hyperpolarizabilities from numerical derivatives

*i*axis) tensor components of the polarizability and second hyperpolarizability read

^{43}

*μ*and

*β*. Equations (6) and (7) demonstrate that the finite difference expression approximates

*α*

_{ii}and

*γ*

_{iiii}with contaminations that grow as the second, fourth, … powers of the field amplitude and that these contaminations involve the higher-order responses of the same parity (here, even-order responses). Obtaining $\gamma iiiiEi$ values is therefore quite straightforward, provided field-dependent energies can be evaluated. This can usually be performed with a broader range of methods than in the case of analytical derivative schemes. From now on, such a technique (i.e., the fourth-order derivative of the energy to get

*γ*) is named the “

*pure numerical differentiation scheme*.” Yet, the second hyperpolarizability can be evaluated from field-dependent polarizabilities or first hyperpolarizabilities, calculated using an analytical scheme. This procedure is referred to as the “

*hybrid numerical differentiation scheme*,” and typical finite difference equations read

*ɛ*

_{iiiiii}, …) so that obtaining precise values requires that the field amplitudes need to be as small as possible. Conversely, the smaller the field amplitude, the lower the precision of the required derivatives because the number of significant digits in the field-dependent energy differences decreases. A compromise has thus to be found, and to improve the numerical precision, the Richardson extrapolation procedure also known as the Romberg procedure

^{66,73–75}can be employed. Starting from Eq. (7), the combination of estimated $\gamma iiiiEi$ values obtained for electric field amplitudes of

*E*

_{i}and 2

*E*

_{i}provides an improved estimate, $\gamma iiii\u2032Ei$, where the

*ɛ*

_{iiiiii}contamination has been removed,

*E*

^{k}= 2

^{k}

*E*

^{0}relationship, at iteration

*m*, the

*γ*estimate (we have dropped the tensor indices for clarity) reads

*γ*

^{m,k}values obtained in successive iterations are then sorted in a table, which helps monitoring the convergence of the iterative process, and allows us to select the best estimate of the derivative, with its related precision.

### C. Methods to calculate the field-dependent energies and molecular properties

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.

#### 1. Wavefunction methods^{76}

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}), $\Psi CC=eT\u0302\Psi HF$, where $T\u0302$ 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\u0302$ 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}

#### 2. Density functional theory (DFT)

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 $\rho r$ so that the energy of the system is described as a functional of $\rho 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.^{79} These drawbacks have been related to the self-interaction error or over-delocalization of the response to external electric fields.^{65}

*r*

_{12}).

^{80}For the purpose of calculating responses to electric fields, these are

*long-range*corrected (LC) hybrids, where the amount of Hartree–Fock exchange increases with

*r*

_{12}. Hence, 1/

*r*

_{12}is split into short- and long-range parts by using a standard error function to allow a smooth evolution with

*r*

_{12},

^{81}

^{,}

*α*corresponds to the amount of HF exchange at zero interelectronic distance, while

*α*+

*β*corresponds to the amount in the asymptotic limit.

*μ*is the range-separating parameter that controls their

*r*

_{12}-dependence. RSH functionals have been shown to correctly describe the asymptotic behavior and are among the best functionals to predict the NLO properties.

^{25,48,50,53,82,83}Moreover, it has been demonstrated that another improvement can be achieved by optimizing the range-separating parameter

*μ*with fixed fraction of the HF exchange in the short-range and/or by applying a two-dimensional tuning procedure where two parameters (

*α*and

*μ*) in Eq. (12) are varied (with a constrain that

*β*= 1 −

*α*in the LC functional). Among these strategies, the most popular non-empirical tuning of RSH functionals consists in selecting

*μ*in such a way that Koopmans’ theorem is satisfied for the first ionization energy.

^{84,85}In other words, the following difference, $\Delta IP\mu $, should be as close as possible to zero:

^{86–88}These studies have shown that the optimum

*μ*value is system dependent. The advantage of using non-empirically tuned RSHs to evaluate hyperpolarizabilities (generally, this was achieved for

*β*) has, however, not been confirmed as a general trend, and some authors have concluded that the tuning scheme in Eq. (13) does not improve the accuracy of NLO properties,

^{53}especially second hyperpolarizabilities.

^{25,89}In this context, Besalú-Sala and coworkers

^{83}recently proposed a new tuned RSH functional, T

_{α}-LC-BLYP, in which the “optimal”

*μ*parameter is determined using an empirical correlation between the static polarizability (the diagonal component along the main inertial axis, as calculated with the original LC-BLYP functional) and the values of

*μ*that reproduces the CCSD(T) second hyperpolarizability for 60 medium-size organic molecules.

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.

^{92}The $EN$ curvature can be computed explicitly using electronic structure packages that are capable of treating fractional electron numbers, but it requires modifications in standard codes.

^{93,94}Alternatively, as done here, it can be interpolated by using the cubic spline interpolation formula,

^{94}

^{,}

*et al.*,

^{94}the quantity in the curly brackets, $\Delta \Delta E\delta $, is the deviation with respect to the linear behavior in

*δ*, as described by $\Delta E\delta $. This expression has already been employed to analyze the localization and delocalization errors along the use of the HF and DFT methods for calculating molecular properties.

^{49,54,95,96}

#### 3. Response theory methods

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.^{97} The expansion terms are the successive linear, quadratic, cubic, … response functions, which can be static/dynamic.^{98} Hence, the first-order dipole moment at pulsation *ω* induced by an electric field at the same pulsation defines the frequency-dependent polarizability, $\u2212\u27e8\u27e8\mu i;\mu j\u27e9\u27e9\omega =\alpha ij\u2212\omega ;\omega $. Similarly, when applying electric fields at pulsations *ω*_{1} and *ω*_{2}, the second-order induced dipole, oscillating at the sum frequency *ω*_{σ} = *ω*_{1} + *ω*_{2}, defines the dynamic first hyperpolarizability, $\u2212\u27e8\u27e8\mu i;\mu j,\mu k\u27e9\u27e9\omega 1,\omega 2=\beta ijk\u2212\omega \sigma ;\omega 1,\omega 2$ and so on for $\u2212\u27e8\u27e8\mu i;\mu j,\mu k,\mu l\u27e9\u27e9\omega 1,\omega 2,\omega 3=\gamma ijkl\u2212\omega \sigma ;\omega 1,\omega 2,\omega 3$, where $\omega \sigma =\u2211i1,2,\u2026\omega i$ is the frequency of the generated light. These response functions can be evaluated at different levels of approximation.^{38} At the HF level, it leads to the time-dependent Hartree–Fock (TDHF)^{39} scheme, while its DFT counterpart is known as the time-dependent DFT (TDDFT)^{40} 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}

### D. Frequency dispersion of the second hyperpolarizability

^{70}Bishop and co-workers

^{26,101,102}demonstrated relationships between the static and frequency-dependent hyperpolarizabilities, which are satisfied for the diagonal tensor components as well as for rotational averages. In the case of the second hyperpolarizability, such relationships read

*A*,

*B*, … are frequency independent expansion coefficients. Moreover, they proved that

*A*is the same for all third-order processes, though it depends on the molecule. In Eq. (15), $\omega Ln$ (

*n*= 2, 4, …) is an “effective” frequency at power “

*n*” and is defined as

*ω*

_{i}quantities encompass the incident and generated light pulsations. Therefore, for typical NLO processes where the electric fields are static or oscillate at the same frequency

*ω*, $\omega Ln=\u2113\omega n$, with

*ℓ*an integer. For instance, for the leading term (

*n*= 2),

*ℓ*= 2 and $\omega L2=2\omega 2$ in the case of the optical Kerr effect (OKE) process [$\gamma \u2212\omega ;\omega ,0,0$],

*ℓ*= 6 for electric-field induced second harmonic generation [$\gamma \u22122\omega ;\omega ,\omega ,0$],

*ℓ*= 4 for degenerate four-wave mixing [$\gamma \u2212\omega ;\omega ,\u2212\omega ,\omega $], and

*ℓ*= 12 for third harmonic generation/scattering [$\gamma \u22123\omega ;\omega ,\omega ,\omega $]. Hence, if the

*A*coefficient is determined from the frequency dependence of the OKE response, it can then be used to estimate, at first approximation (i.e., in the low frequency regime), the THG response.

## III. COMPUTATIONAL ASPECTS

The geometry of *p*NA was fully optimized at the MP2/6-311G(d) level in gas phase, with thresholds on the forces and displacements of 1.5 × 10^{−5} hartree bohr^{−1} and 6.0 × 10^{−5} Å, 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^{−11} 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 [$\alpha \u2212\omega ;\omega ,\beta \u2212\omega ;\omega ,0,and\beta \u22122\omega ;\omega ,\omega $] 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^{−10} a.u. These calculations were carried out using the Gaussian16 program.^{103} Moreover, cubic response functions (CRFs) have been calculated at the CCSD level using the Dalton program^{104} using the unrelaxed field-independent orbitals.

Type . | Acronym . | % HF . | μ (L)
. | % PT2 . |
---|---|---|---|---|

LDA | SVWN | ⋯ | ⋯ | ⋯ |

GGA | BLYP | ⋯ | ⋯ | ⋯ |

PBE | ⋯ | ⋯ | ⋯ | |

B97-D | ⋯ | ⋯ | ⋯ | |

mGGA | M06-L | ⋯ | ⋯ | ⋯ |

Global hybrid GGA | B3LYP | 20 | ⋯ | ⋯ |

PBE0 | 25 | ⋯ | ⋯ | |

Global hybrid mGGA | M06 | 27 | ⋯ | ⋯ |

M06-2X | 54 | ⋯ | ⋯ | |

M06-HF | 100 | ⋯ | ⋯ | |

MN15 | 44 | ⋯ | ⋯ | |

Range-separated hybrid GGA | ωB97 | 0–100 | 0.40 (1.323) | ⋯ |

ωB97X | 15.7–100 | 0.30 (1.764) | ⋯ | |

ωB97X-D | 22.2–100 | 0.20 (2.646) | ⋯ | |

LC-ωPBE | 0–100 | 0.40 (1.323) | ⋯ | |

CAM-B3LYP | 19–65 | 0.33 (1.604) | ⋯ | |

LC-BLYP | 0–100 | 0.47 (1.126) | ⋯ | |

T_{α}-LC-BLYP | 0–100 | 0.30 (1.764) | ⋯ | |

Range-separated hybrid mGGA | M11 | 42.8–100 | 0.25 (2.117) | ⋯ |

Double hybrid GGA | PBE0-DH | 50 | ⋯ | 12.5 |

mPW2-PLYP | 55 | ⋯ | 25 | |

B2-PLYP | 53 | ⋯ | 27 |

Type . | Acronym . | % HF . | μ (L)
. | % PT2 . |
---|---|---|---|---|

LDA | SVWN | ⋯ | ⋯ | ⋯ |

GGA | BLYP | ⋯ | ⋯ | ⋯ |

PBE | ⋯ | ⋯ | ⋯ | |

B97-D | ⋯ | ⋯ | ⋯ | |

mGGA | M06-L | ⋯ | ⋯ | ⋯ |

Global hybrid GGA | B3LYP | 20 | ⋯ | ⋯ |

PBE0 | 25 | ⋯ | ⋯ | |

Global hybrid mGGA | M06 | 27 | ⋯ | ⋯ |

M06-2X | 54 | ⋯ | ⋯ | |

M06-HF | 100 | ⋯ | ⋯ | |

MN15 | 44 | ⋯ | ⋯ | |

Range-separated hybrid GGA | ωB97 | 0–100 | 0.40 (1.323) | ⋯ |

ωB97X | 15.7–100 | 0.30 (1.764) | ⋯ | |

ωB97X-D | 22.2–100 | 0.20 (2.646) | ⋯ | |

LC-ωPBE | 0–100 | 0.40 (1.323) | ⋯ | |

CAM-B3LYP | 19–65 | 0.33 (1.604) | ⋯ | |

LC-BLYP | 0–100 | 0.47 (1.126) | ⋯ | |

T_{α}-LC-BLYP | 0–100 | 0.30 (1.764) | ⋯ | |

Range-separated hybrid mGGA | M11 | 42.8–100 | 0.25 (2.117) | ⋯ |

Double hybrid GGA | PBE0-DH | 50 | ⋯ | 12.5 |

mPW2-PLYP | 55 | ⋯ | 25 | |

B2-PLYP | 53 | ⋯ | 27 |

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 $\gamma THS=\u27e8\gamma ZZZZ2\u27e9+\u27e8\gamma ZXXX2\u27e91/2$ and $DR=\u27e8\gamma ZZZZ2\u27e9/\u27e8\gamma ZXXX2\u27e9$ quantities of THS.^{16,17,34} All second hyperpolarizabilities are reported in atomic units (1 a.u. of γ = 6.235 377 × 10^{−65} C^{4} m^{4} J^{−3} = 5.0367 × 10^{−40} esu), within the T convention, as defined in Eq. (1).^{105}

### A. Integration grid in DFT calculations

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^{4}–10^{5} 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.

### B. Basis set effects

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^{109} [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^{110} [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^{3} a.u., *γ*_{//} ∼ 32 × 10^{3} a.u., *γ*_{THS} ∼ 43 × 10^{3} 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.

Basis sets . | γ_{xxxx}
. | γ_{//}
. | γ_{THS}
. | DR . |
---|---|---|---|---|

6-31G(d) [162] | 71 459 (−35) | 13 886 (−56) | 24 968 (−42) | 6.6 (−40) |

6-31G(d,p) [180] | 71 404 (−35) | 13 867 (−56) | 24 943 (−42) | 6.6 (−40) |

6-311G(d) [198] | 76 421 (−31) | 15 688 (−51) | 26 921 (−37) | 7.0 (−36) |

6-311G(d,p) [216] | 76 140 (−31) | 15 616 (−51) | 26 813 (−37) | 7.0 (−36) |

cc-pVDZ [170] | 68 950 (−37) | 13 949 (−56) | 24 237 (−43) | 6.9 (−37) |

cc-pVTZ [384] | 77 375 (−30) | 16 307 (−49) | 27 403 (−36) | 7.2 (−34) |

cc-pVQZ [730] | 84 397 (−23) | 18 934 (−40) | 30 288 (−29) | 7.7 (−30) |

6-31+G(d) [202] | 108 456 (−2) | 29 271 (−8) | 40 762 (−5) | 9.8 (−11) |

6-31+G(d,p) [220] | 108 418 (−2) | 29 282 (−8) | 40 756 (−5) | 9.8 (−11) |

6-311+G(d) [238] | 108 582 (−1) | 28 817 (−9) | 40 648 (−5) | 9.5 (−13) |

6-311+G(d,p) [256] | 108 101 (−2) | 28 690 (−10) | 40 467 (−6) | 9.5 (−13) |

aug-cc-pVDZ [234] | 105 619 (−4) | 29 318 (−10) | 40 389 (−6) | 10.3 (−6) |

aug-cc-pVTZ [518] | 107 436 (−3) | 30 485 (−4) | 41 493 (−3) | 10.7 (−3) |

aug-cc-pVQZ [1076] | 108 725 (−1) | 31 129 (−2) | 42 168 (−2) | 10.9 (−1) |

d-aug-cc-pVDZ [398] | 111 708 (1) | 32 190 (1) | 43 469 (1) | 10.9 (−1) |

d-aug-cc-pVTZ [812] | 110 273 (−0.2) | 31 787 (0.02) | 42 913 (−0.1) | 11.0 (0) |

d-aug-cc-pVQZ [1422] | 110 455 | 31 782 | 42 954 | 11.0 |

Basis sets . | γ_{xxxx}
. | γ_{//}
. | γ_{THS}
. | DR . |
---|---|---|---|---|

6-31G(d) [162] | 71 459 (−35) | 13 886 (−56) | 24 968 (−42) | 6.6 (−40) |

6-31G(d,p) [180] | 71 404 (−35) | 13 867 (−56) | 24 943 (−42) | 6.6 (−40) |

6-311G(d) [198] | 76 421 (−31) | 15 688 (−51) | 26 921 (−37) | 7.0 (−36) |

6-311G(d,p) [216] | 76 140 (−31) | 15 616 (−51) | 26 813 (−37) | 7.0 (−36) |

cc-pVDZ [170] | 68 950 (−37) | 13 949 (−56) | 24 237 (−43) | 6.9 (−37) |

cc-pVTZ [384] | 77 375 (−30) | 16 307 (−49) | 27 403 (−36) | 7.2 (−34) |

cc-pVQZ [730] | 84 397 (−23) | 18 934 (−40) | 30 288 (−29) | 7.7 (−30) |

6-31+G(d) [202] | 108 456 (−2) | 29 271 (−8) | 40 762 (−5) | 9.8 (−11) |

6-31+G(d,p) [220] | 108 418 (−2) | 29 282 (−8) | 40 756 (−5) | 9.8 (−11) |

6-311+G(d) [238] | 108 582 (−1) | 28 817 (−9) | 40 648 (−5) | 9.5 (−13) |

6-311+G(d,p) [256] | 108 101 (−2) | 28 690 (−10) | 40 467 (−6) | 9.5 (−13) |

aug-cc-pVDZ [234] | 105 619 (−4) | 29 318 (−10) | 40 389 (−6) | 10.3 (−6) |

aug-cc-pVTZ [518] | 107 436 (−3) | 30 485 (−4) | 41 493 (−3) | 10.7 (−3) |

aug-cc-pVQZ [1076] | 108 725 (−1) | 31 129 (−2) | 42 168 (−2) | 10.9 (−1) |

d-aug-cc-pVDZ [398] | 111 708 (1) | 32 190 (1) | 43 469 (1) | 10.9 (−1) |

d-aug-cc-pVTZ [812] | 110 273 (−0.2) | 31 787 (0.02) | 42 913 (−0.1) | 11.0 (0) |

d-aug-cc-pVQZ [1422] | 110 455 | 31 782 | 42 954 | 11.0 |

### C. Romberg’s scheme and numerical precision

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} = 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^{1} 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.

## IV. RESULTS AND DISCUSSION

### A. Electron correlation effects with wavefunction methods

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 *p*NA 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.

Methods . | γ_{//}
. | γ_{THS}
. | DR . |
---|---|---|---|

HF | 19 534 (−49) | 26 280 (−52) | 10.2 (10) |

MP2 | 33 679 (−11) | 49 649 (−10) | 8.7 (−5) |

MP3 | 27 677 (−27) | 40 409 (−27) | 8.8 (−4) |

MP4 | 36 778 (−3) | 53 475 (−3) | 9.1 (−1) |

MP4D | 30 072 (−21) | 44 606 (−19) | 8.5 (−8) |

MP4DQ | 26 643 (−30) | 38 581 (−30) | 8.9 (−3) |

MP4SDQ | 30 175 (−21) | 43 305 (−22) | 9.2 (0) |

CCSD(FF) | 32 214 (−15) | 46 457 (−16) | 9.2 (0) |

CCSD(CRF) | 36 959 (−3) | 54 057 (−2) | 8.3 (−10) |

CCSD(T) | 37 986 | 55 163 | 9.2 |

Methods . | γ_{//}
. | γ_{THS}
. | DR . |
---|---|---|---|

HF | 19 534 (−49) | 26 280 (−52) | 10.2 (10) |

MP2 | 33 679 (−11) | 49 649 (−10) | 8.7 (−5) |

MP3 | 27 677 (−27) | 40 409 (−27) | 8.8 (−4) |

MP4 | 36 778 (−3) | 53 475 (−3) | 9.1 (−1) |

MP4D | 30 072 (−21) | 44 606 (−19) | 8.5 (−8) |

MP4DQ | 26 643 (−30) | 38 581 (−30) | 8.9 (−3) |

MP4SDQ | 30 175 (−21) | 43 305 (−22) | 9.2 (0) |

CCSD(FF) | 32 214 (−15) | 46 457 (−16) | 9.2 (0) |

CCSD(CRF) | 36 959 (−3) | 54 057 (−2) | 8.3 (−10) |

CCSD(T) | 37 986 | 55 163 | 9.2 |

### B. Performance of DFT functionals

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 *p*NA is challenging for all types of functionals, whether conventional, global hybrids, or RSHs.

XCFs . | γ_{//}
. | γ_{THS}
. | DR . |
---|---|---|---|

SVWN | 28 757 (−24) | 37 149 (−33) | 12.3 (34) |

BLYP | 32 129 (−15) | 41 369 (−25) | 12.4 (35) |

PBE | 31 136 (−18) | 40 248 (−27) | 12.3 (34) |

B97-D | 30 584 (−19) | 39 589 (−28) | 12.1 (31) |

M06-L | 25 192 (−34) | 35 136 (−36) | 9.8 (6) |

B3LYP | 31 368 (−17) | 43 286 (−22) | 10.1 (10) |

PBE0 | 29 811 (−22) | 41 572 (−25) | 9.8 (6) |

M06 | 28 589 (−25) | 40 148 (−27) | 9.6 (4) |

M06-2X | 26 583 (−30) | 37 883 (−31) | 9.3 (1) |

M06-HF | 23 380 (−39) | 32 799 (−41) | 9.6 (4) |

M11 | 24 563 (−35) | 35 585 (−35) | 9.0 (−2) |

MN15 | 29 709 (−22) | 41 840 (−24) | 9.6 (4) |

ωB97 | 25 342 (−33) | 35 697 (−35) | 9.5 (3) |

ωB97X | 26 580 (−30) | 37 003 (−33) | 9.8 (6) |

ωB97X-D | 27 789 (−27) | 38 648 (−30) | 9.8 (6) |

LC-ωPBE | 24 624 (−35) | 35 254 (−36) | 9.2 (0) |

CAM-B3LYP | 28 856 (−24) | 40 675 (−26) | 9.5 (3) |

LC-BLYP | 24 965 (−34) | 35 892 (−35) | 9.1 (−1) |

PBE0DH | 30 886 (−19) | 44 596 (−19) | 9.1 (−1) |

mPW2PLYP | 35 642 (−6) | 51 823 (−6) | 9.0 (−2) |

B2-PLYP | 37 272 (−2) | 54 121 (−2) | 9.1 (−1) |

CCSD(T) | 37 986 | 55 163 | 9.2 |

XCFs . | γ_{//}
. | γ_{THS}
. | DR . |
---|---|---|---|

SVWN | 28 757 (−24) | 37 149 (−33) | 12.3 (34) |

BLYP | 32 129 (−15) | 41 369 (−25) | 12.4 (35) |

PBE | 31 136 (−18) | 40 248 (−27) | 12.3 (34) |

B97-D | 30 584 (−19) | 39 589 (−28) | 12.1 (31) |

M06-L | 25 192 (−34) | 35 136 (−36) | 9.8 (6) |

B3LYP | 31 368 (−17) | 43 286 (−22) | 10.1 (10) |

PBE0 | 29 811 (−22) | 41 572 (−25) | 9.8 (6) |

M06 | 28 589 (−25) | 40 148 (−27) | 9.6 (4) |

M06-2X | 26 583 (−30) | 37 883 (−31) | 9.3 (1) |

M06-HF | 23 380 (−39) | 32 799 (−41) | 9.6 (4) |

M11 | 24 563 (−35) | 35 585 (−35) | 9.0 (−2) |

MN15 | 29 709 (−22) | 41 840 (−24) | 9.6 (4) |

ωB97 | 25 342 (−33) | 35 697 (−35) | 9.5 (3) |

ωB97X | 26 580 (−30) | 37 003 (−33) | 9.8 (6) |

ωB97X-D | 27 789 (−27) | 38 648 (−30) | 9.8 (6) |

LC-ωPBE | 24 624 (−35) | 35 254 (−36) | 9.2 (0) |

CAM-B3LYP | 28 856 (−24) | 40 675 (−26) | 9.5 (3) |

LC-BLYP | 24 965 (−34) | 35 892 (−35) | 9.1 (−1) |

PBE0DH | 30 886 (−19) | 44 596 (−19) | 9.1 (−1) |

mPW2PLYP | 35 642 (−6) | 51 823 (−6) | 9.0 (−2) |

B2-PLYP | 37 272 (−2) | 54 121 (−2) | 9.1 (−1) |

CCSD(T) | 37 986 | 55 163 | 9.2 |

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^{−1}), to ωB97X (*μ* = 0.3 bohr^{−1}), and ωB97X-D (*μ* = 0.2 bohr^{−1}), 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.^{116} As defined in Refs. 17 and 34, the anisotropy parameters $\rho 0/2=\gamma J=0/\gamma J=2$ and $\rho 4/2=\gamma J=4/\gamma 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.

Methods . | ρ_{0/2}
. | ρ_{4/2}
. |
---|---|---|

SVWN | 1.08 | 0.71 |

B3LYP | 0.93 | 0.73 |

M06-2X | 0.88 | 0.75 |

CAM-B3LYP | 0.90 | 0.75 |

LC-BLYP | 0.86 | 0.76 |

B2-PLYP | 0.84 | 0.71 |

HF | 1.01 | 0.86 |

CCSD(T) | 0.84 | 0.69 |

Methods . | ρ_{0/2}
. | ρ_{4/2}
. |
---|---|---|

SVWN | 1.08 | 0.71 |

B3LYP | 0.93 | 0.73 |

M06-2X | 0.88 | 0.75 |

CAM-B3LYP | 0.90 | 0.75 |

LC-BLYP | 0.86 | 0.76 |

B2-PLYP | 0.84 | 0.71 |

HF | 1.01 | 0.86 |

CCSD(T) | 0.84 | 0.69 |

### C. Tuning the percentages of HF exchange and PT2 correlation in B2-PLYP

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 *p*NA molecule.

### D. Tuning the range-separating parameter

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 $\Delta IP\mu $ function [Eq. (13)] and its components [Fig. 5(a) and Table S8]. $\Delta IP\mu $ decreases monotonically from *μ* = 0.05 bohr^{−1} until it cancels at 0.28 bohr^{−1} (0.26 bohr^{−1}) for LC-BLYP (ωB97), which corresponds therefore to the value where Koopmans’ theorem is satisfied. Then, $\Delta IP\mu $ gets more and more negative. For CAM-B3LYP, $\Delta IP\mu $ 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^{−1} 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^{−1}, 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.*^{83} dedicated specifically to the calculation of second hyperpolarizabilities. Owing to the value of the LC-BLYP/6-311+G* longitudinal polarizability of *p*NA, *α*_{xx} = 135.15 a.u.., the size extensive descriptor of Ref. 83, amounts to 0.274, which corresponds to *μ*_{Tα} = 0.30 bohr^{−1}. If considering the average polarizability of *p*NA instead of its longitudinal component, the *μ*_{Tα} value becomes 0.33 bohr^{−1}. Both *μ*_{Tα} values are thus larger than the optimal one obtained using Eq. (13) and, therefore, do not improve the *γ* value of *p*NA. This is probably because no *π*-conjugated push-pull systems were included in the molecular dataset considered in Ref. 83.

### E. Delocalization error of XCFs

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 $\Delta E\delta $ quantity as a function of *δ* as well as its deviations with respect to the exact linear behavior described by $\Delta E\delta $ [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 $\Delta E\delta $ 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^{−1} and *μ* = 0.40 bohr^{−1}). 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^{−1}) 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^{−1} 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.

### F. Frequency dispersion effects

*ℏω*= 0.0–1.9 eV, $\gamma //\u22122\omega ;\omega ,\omega ,0$ increases by about one order of magnitude, while the enhancement of $\gamma //\u2212\omega ;\omega ,0,0$ amounts only to 50%. At

*ℏω*= 0.08 a.u. = 2.18 eV, and again at

*ℏω*= 0.09 a.u. = 2.45 eV, $\gamma //\u22122\omega ;\omega ,\omega ,0$ changes sign (see Table S11). These sign changes are consistent with the values of the two first excitation energies of

*p*NA, evaluated at 3.96 and 4.38 eV at the same level of approximation. Indeed, $\gamma //\u22122\omega ;\omega ,\omega ,0$ increases with the frequency and its frequency dispersion becomes stronger when approaching $\u210f\omega =4.382eV=2.19eV$ owing to the two-photon resonance, before it changes sign. Plotting

*γ*

_{//}as a function of $\omega L2$ [Fig. 8(b)] shows that the curves are superimposed at low values of $\omega L2$, as predicted in Refs. 101 and 102. Equation (15) was then fitted to each set of

*γ*

_{//}values (before resonance), leading to the following expression for the frequency dispersion of

*γ*

_{//}at small frequencies:

## V. FURTHER DISCUSSIONS, CONCLUSIONS, AND OUTLOOK

Before concluding on the performance of wavefunctions and DFT methods to evaluate the second hyperpolarizability of *p*NA, 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^{−1} for LC-BLYP and CAM-B3LYP and to 0.14 bohr^{−1} for ωB97, which differ from the *μ* values that satisfy $\Delta IP\mu =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 *p*NA 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 *p*NA 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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

*Second-Order Nonlinear Optical Characterization Techniques*

*Advances in Chemical Physics*

*Density-Functional Theory of Atoms and Molecules*

*Molecular Electronic‐Structure Theory*

*Propagators in Quantum Chemistry*