Efficient Kohn–Sham density-functional theory implementation of isotropic spectroscopic observables associated with quadratic response functions

For general exchange–correlation functionals with a dependence on the local spin densities and spin-density gradients, we provide computationally tractable expressions for the tensor-averaged quadratic response functions pertinent to the experimental observables in second-harmonic generation (SHG). We demonstrate how the tensor-averaged quantities can be implemented with reference to a derived minimal number of first- and second-order perturbed Fock matrices. Our consideration has the capability of treating a situation of resonance enhancement as it is based on damped response theory and allows for the evaluation of tensor-averaged resonant-convergent quadratic response functions using only ∼25% (one-photon off-resonance regions) and ∼50% (one-photon resonance regions) of the number of auxiliary Fock matrices required when explicitly calculating all the needed individual tensor components. Numerical examples of SHG intensities in the one-photon off-resonance region are provided for a sample of makaluvamine derivatives recognized for their large nonlinear optical responses as well as a benchmark set of small- and medium-sized organic molecules.


Introduction
Discovered in 1961 [1], secondary harmonic generation (SHG) is a scattering process in which light is observed at twice the frequency of the incident radiation. Since then and based on this process, spectroscopic methods and applied technologies have been developed. The SHG process enables the distinction between surface-bound molecule and molecules that are randomly oriented in the bulk phase as the latter do not produce a second harmonic signal due to destructive interference. SHG spectroscopy has been used to study the surface molecular orientations at interfaces [2][3][4] and in microscopy to visualize cell and tissue structure and function [5][6][7][8]. By relating spectroscopic response signals to molecular structure and organization, theoretical simulations play an important role in this work and also in the design of materials that exhibit strong second harmonic scattering.
Response theory is an established framework for calculating nonlinear optical properties of molecular materials with modern formulations and associated implementations based on either the Ehrenfest [9,10] or the quasi-energy approach [11,12]. These developments include situations when the external electromagnetic fields are in resonance with transition frequencies in the molecular system, and with the introduction of phenomenological damping parameters associated with the inverse lifetime of excited states, the resulting response functions are physically sound in the entire spectral region and they are in general complex-valued. However, the accurate determination of frequency-dependent molecular hyperpolarizabilities depends critically on the treatment of electron correlation [13]. As to provide a reasonable compromise between computational cost and accuracy, response theory at the level of density functional theory (DFT) has been developed [14,15], and the efficient and numerically stable solution of the

EFISHG and hyperpolarizabilities
Molecular polarizabilites and hyperpolarizabilites can be derived using response theory as corrections to the expectation value of the electric-dipole moment operator in the presence of an external electric field [19]. In order to describe nonlinear optical phenomenon, one assumes that the first moment of the charge distribution, i.e. the dipole moment, can be expressed in a power series expansion in the electric field strength. Under this assumption, the time-dependent dipole moment can have Fourier amplitudes oscillating at frequencies differing from that of the incident light. The expansion coefficients in this power series of the time-dependent dipole moment defines the response functions: The coefficient linear in the field strength, the polarizability α, determines the linear response and corresponds to phenomenon such as one-photon absorption, the coefficients for the squared field strength, the first-order nonlinear hyperpolarizability β, is related to the quadratic response functions, which is of main importance for this work as it relates to the process of second-harmonic generation. For an EFISHG experiment, a static electric field is applied along with a time varying field, subsequently when two quanta of light interact with the molecules, scattering of photons of double the frequency may be observed. If one assumes that the sample is isotropic, the measured second-harmonic intensity is related to Γ as: ⟨γ(−2ω; ω, ω, 0)⟩ = 1 15 x,y,z α,β (γ ααββ + γ αββα + γ αβαβ ) .
The EFISHG experiment measures the isotropic second-order dipole hyperpolarizability given by ⟨γ⟩. In dipolar samples, there is also a temperature dependent term that depends on the vectorial component of the first-order hyperpolarizability β tensor given as [20][21][22]:

Hyperpolarizabilities and response functions
Within the quasi-energy formulation, the equation that governs the responses of the wave function parameters with respect to the perturbation is the time-dependent variational principle for the time-averaged quasi-energy [23]: One can show that the quasi-energy satisfies the time-dependent Hellman-Feynman theorem, where the time-averaged quasi-energy is defined in terms of the time-independent HamiltonianĤ 0 , and the time-dependent perturbation operatorV(t) [11]: In the presence of the field, the quasi-energy is written as a power series expansion of Fourier amplitudes of the field: Combining the time-dependent Hellman-Feynman theorem, equation (5), and the expansion of the dipole moment, equation (1), it can be shown that [12]: When the time averaging is carried out one arrives at the conclusion that ω 1 must equal minus the sum of the other optical frequencies in order to get a nonzero result.

Kohn-Sham reference state parameterization
The time evolution of the Kohn-Sham determinant is parameterized in terms of the time-dependent orbital rotations: The anti-Hermitian orbital-rotation operatorκ can be written as an inner product: where the excitation and de-excitation operators are given by: With indices i, j, . . ., a, b, . . ., and p, q, . . . we denote occupied, unoccupied, and general spin orbitals or spatial orbitals when the operators are accompanied with a spin index, respectively. We also introduce compound indices, n, m, . . ., to denote pairs of orbital indices {a, i}. Furthermore, we use Greek letters for Cartesian coordinates in the molecular frame of reference. Using equation (10), we define the state-rotation parameter vector in two blocks as: Using the Baker-Campbell-Hausdorff expansion one can expand the density in terms of the orbital-rotation parameters as [24]: (13) where the density and density gradient operator are defined as [25]:

Response functions and orbital responses
In the presence of the perturbation we write the order-corrections to the orbital-rotation parameters as, The first and second-order corrections to the κ-vectors can then be constructed in terms of the Fourier amplitudes of the field as: κ (2) Using the parametrization of equation (9), we can get working expressions for the order-corrections to the quasi-energy of equation (7) with respect to orbital-rotation vectors as: The evaluation of the order-corrections to the quasi-energy requires the order-corrections to the orbital rotation parameters, which are themselves solutions of the response equations. With use of the time-dependent variational principle equation (4), we get the linear and quadratic response equations: Solving for the second-order response vector in equations (19) and (20) where we have symmetrized the response vectors in frequency variables ω 1 , ω 2 we get: The quadratic response function being part of the third-order correction to the time-averaged quasi-energy can by the 2n + 1 rule be written in terms of the first-order response of the wave function parameters. If we substitute equation (22) into equation (18) and symmetrize with respect to ω 2 and ω 3 we get the working formula for the quadratic response function as: where again ω σ = (ω 2 + ω 3 ) and where we have made the substitution: Some important features to note here are that the evaluation of the quadratic response function requires the contraction of the tensor ∂ 3 QT ∂κ −ωσ ∂κ ω 1 ∂κ ω 2 which contains the electronic cubic-Hessian, ∂κ −ωσ ∂κ ω 1 ∂κ ω 2 , see equations (A.10) and (A.11). The two-electron part of the electronic Hessian has been treated elsewhere, here we are only interested in the exchange-correlation contribution to the electronic Hessian. In summary, by solving the response equations equations (19) and (20) we obtain the corrections to the orbital-rotation parameters in the field. In turn, we can use these to find the corrections to the density: Using the chain rule and the orbital rotation parametrization of the density of equation (13) and the first-order correction to the orbital-rotation parameters, we can get the Fourier amplitudes of the first-order correction to the density as: Differentiating the density expansion of equation (13) we get the derivatives: Combining equations (26) and (27), we get that the perturbed densities can be written as: where we have defined the operator:κ

Contraction of the second-order Hessian with a response vector
Owing to the large size of the electronic Hessian, the generalized gradient vector formed from the contraction of the electronic Hessian with a response vector is derived analytically. Here we derive the transformed exchange-correlation contribution to the generalized gradient vector required for the exchange-correlation part of equation (19): Using the expression for the exchange-correlation contribution to the electronic Hessian of equation (A.9), for a functional of the form of equation (A.1), the lower part of the resultant vector on the right hand-side can be written as: where the perturbed density is given by equation (28). The perturbed exchange-correlation matrix of equation (31) can be seen as composed of two parts, one where the κ operator of equation (29) is commuted with the zeroth-order exchange-correlation operator of equation (A.7) and one term which only depends on the first-order exchange-correlation operator of equation (33): The first term,v ω1 xc;σ , in equation (32) is what we call the first-order exchange-correlation matrix operator and is defined by combining equations (28), (31) and (32): where an element of the tensor v (1) xc;pq,rs;σ,σ ′ in equation (33) is given by:

Contraction of the third-order Hessian with two response vectors
In the direct approach the cubic electronic-Hessian, a rank three tensor, is not computed explicitly, instead we here derive the analytical expression for the generalized gradient formed from contracting the electronic cubic-Hessian with two response vectors: where each element of the double transformed exchange-correlation matrices are given in terms of the permutation operator P ω1,ω2 as: Using the expression for the cubic-Hessian of equation (A.11) with the expression for the density expansion of equation (13), we then get that equation (36) can be written as: where we have that the second-order perturbed density is given by the contraction of the second-order derivative of the density expansion of equation (13) with respect to the orbital rotation parameters with the corrections to the orbital rotation parameters in the presence of the perturbation: We can write equation (37) compactly using the expressions for the first and second-order exchange-correlation operators as, The two-time transformed exchange-correlation matrix of equation (39) can be used to define the second-order exchange-correlation operator: where v (1) xc,pq,ai;σ,σ ′ , is given by equation (34), and we define v (2) xc,pq,tu,rs;σ,σ ′ ,σ ′ ′ as: v (2) xc;pq,tu,rs;σ,σ ′ ,σ ′ ′ =ˆ where g pqia represents an element of the two-electron integral tensor and D refers to perturbed density matrices. Sums of transformed Fock matrices can be computed by first summing the transformed density matrices: The linear transformation of equation (44) allows for the evaluation of sums of Fock matrices at the same computational cost of a single Fock matrix. Within the DFT frame-work, transformed Fock matrices, can with use of equations (34) and (41) be constructed using v (1) xc,pq,ai;σ,σ ′ and v (2) xc,pq,tu,ai;σ,σ ′ ,σ ′ ′ which are quantities that are in-dependent of the perturbed densities: Using the expressions for the exchange-correlation kernels of equations (34) and (41), and the expression for the DFT Fock matrix, equation (46), it is clear that sums of two-time transformed Fock matrices at the DFT level of theory are linear transformations of the perturbed densities: One can thus in principle reduce the number of kernel integrations required in order to compute quadratic response functions at the DFT level of theory by computing sums of Fock matrices instead of computing them individually.

Compounded Fock matrices for the vector components of first-order hyperpolarizability
In this section, we will derive the compounded Fock matrices required for computing the vector components of the β tensor. With use of equation (47) one can reduce the number of kernel integrations if one can reformulate the vector components of β in terms of sums of Fock matrices. Combining the expression for the isotropic average, equation (3) with the last term of the expression for the quadratic response function equation (23) and the contraction of the cubic-Hessian, equation (35), we get that the total contraction of the cubic-Hessian can be written as: where we have defined the compounded two-time transformed Fock matrices: Recalling that the compound index for the elements of a response vector refers to a pair of virtual (v) and occupied (o) molecular orbitals, see equation (10), we scatter the elements of said vector into the ov-and vo-blocks of a matrix according to: such that Ω λ αβ;σ and Ω π α;σ are given by: Using the expression for the two-time transformed Fock matrices, equation (46), we can write equations (49) and (50) in terms of sums or compounded perturbed densities as: where when using equation (51), the compounded densities are given by: D ω pq;β;σ = κ ω β;σ , D σ pq .

Reduced expression
We can further reduce the number of Fock matrices required to perform the E [3] tensor contractions needed to calculateβ(−2ω; ω, ω). The E [3] contributions to the real and imaginary parts of β α in equation (3) are given by: Hence, one needs both the real F σ,R α and imaginary F σ,I α components of the two-time transformed Fock matrices to form the real and imaginary component ofβ. In one-photon off-resonance regions however, the response vectors with frequency ω are approximately real while the 2ω response vectors might have an imaginary component due to two-photon resonances. This means that the perturbed densities of equations (56)-(58) are all approximately real, hence: One can therefore approximate the expressions of equations (59) and (60) as:

Minimal set of Fock matrices
In this section, a comparison is made between the number of auxiliary Fock matrices required for the calculation ofβ, equation (3) In table 1, the tensor-average algorithm in its full and reduced forms refer to the calculation ofβ using equations (59), (60) and equations (62), (63), respectively. In the tensor-average approach, one does not obtain information about individual β tensor elements but instead directly obtains the elements of the β vector. Using the subspace extraction method of Ahmadzadeh et al [17], we can write the first-order Fock matrices F ω as linear combinations of the electronic Hessian contractions with the vectors that span the subspace and, hence, circumvent the need to recalculate them. As for the compounded two-time transformed Fock matrices of equations (54) and (55) found in equation (48), there are six unique two-time transformed Fock matrices. There are three real and three imaginary Fock matrices of the form F π α , one for each spatial component α. For the F λ αβ two-time transformed Fock matrix, α ̸ = β and the compounded density D λ αβ of equation (56) have permutation symmetry with respect to interchange of α and β. Furthermore, as seen from equation (37), interchanging the perturbed densities for the second-order kernel integration has no effect of the net result, that is p,q,t,u v (2) xc,pq,tu,ai;σ,σ ′ ,σ ′ ′ D ω pq;α;σ ′ D ω tu;β;σ ′ ′ = p,q,t,u v (2) xc,pq,tu,ai;σ,σ ′ ,σ ′ ′ D ω pq;β;σ ′ D ω tu;α;σ ′ ′ , hence F λ αβ = F λ βα and there are thus only three unique F λ αβ Fock matrices. As seen in table 1, the full tensor-average algorithm offers a reduction in the number of kernel integrations and two-electron integral contractions with the perturbed densities by 50%. In regions far from one-photon resonance, where the only linear response vectors used to construct the Fock matrices, are nearly real, the reduced tensor average algorithm allows for the calculation ofβ with 75% fewer Fock matrix constructions as compared to the tensor component algorithm. We will illustrate this method with numerical examples in the following sections.

Computational details
We adopt optimized structures for three makaluvamines (A-C), see figure 1, at the B3LYP [26,27] level of theory with the polarized split-valence def2-SVP [28] basis set. All response calculations were performed using our tensor-average quadratic response implementation in the VeloxChem quantum chemistry software [29] using B3LYP with the def2-SVPD basis set. For all response calculations, a damping parameter of ℏγ = 62.0 meV was employed.

Makaluvamine derivatives A-C
The smallest makaluvamine derivative used, molecule A, which is an isolated cationic pyrroloiminoquinone core moiety, has a strong two-photon resonance at λ = 652 nm corresponding to the S 4 ← S 0 transition and a staticβ 0 = −42.0 a.u. The two remaining molecules, B and C possess an aromatic side-chain that is conjugated with the cationic pyrroloiminoquinone core moiety. The side-chain extended, methyl substituted makaluvamine derivative, molecule B, has a two-photon resonance at λ = 1358 nm corresponding to the S 1 ← S 0 transition and a staticβ 0 = −15 630 a.u. The hydrogen substituted side-chain extended makaluvamine derivative, molecule C, has a two-photon resonance at λ = 1355 nm corresponding to the S 1 ← S 0 transition and a staticβ 0 = −15 170 a.u. The EFISHG spectrum in figure 2 refers to the one-photon off-resonance region, and it is calculated based on the full, equations (59) and (60), as well as the reduced, equations (62) and (63), expressions for the isotropic second-order hyperpolarizability.
As seen in figure 2, results obtained with the full-and reduced-form expressions are in perfect quantitative agreement in the entire spectral regions under consideration. This provides a strong indication that the assumptions underlying the reduced-form expressions in equations (62) and (63) are sound and valid in off-resonance as well as low-lying two-photon resonance regions of the spectrum.

Computational efficiency in quadratic response SHG calculation
The main focus of the present work is to maximize the efficiency of the construction ofβ by reducing the number of auxiliary Fock matrix constructions required for the evaluation of the β-vector components of equation (3).
In this section we discuss the calculations gains obtained from using our algorithm for constructing the EFISHG spectra in figure 2 and compare the full and reduced algorithms. The most time-demanding term for the evaluation of the β-vector components is the contraction of the third-order electronic Hessian of equation (23). Equation (48) represents the total cubic-Hessian contractions required to get the α-component of the β-vector and its real and imaginary components are given by equations (59) and (60).

Calculation of orbital-rotation corrections to first-order
In order to compute the quadratic response function, we require the first-order corrections to the orbital-rotation parameters. For the 35 frequency points used for the construction of the spectra in figure 2, we need to solve response equations of the form equations (21) and (24) to get the first-order corrections to κ: ∂κ ω /∂F ω β and ∂κ −2ω /∂F −2ω α . These equations are solved using a preconditioned iterative subspace approach with symmetrized trial vectors as outlined in the work of Kauczor et al [30,31]. An extension of this algorithm is implemented in our work where multiple frequencies and multiple right-hand sides are solved simultaneously using a common subspace. For the 35 optical frequencies, we need a total of 210 linear response equations out of which 105 are nearly real-valued in one-photon off-resonance regions and 105 are typically complex-valued as they can be close to two-photon resonances.

Contraction of electronic-Hessian
The contraction of the cubic-Hessian is in our double-direct implementation written as transformed Fock matrices as illustrated in equation (48). The terms resulting from the contraction of the cubic Hessian require the construction of one-time time transformed Fock matrices, equation (45) and two-time transformed Fock matrices equation (46). The one-time transformed Fock matrices are found in equations (52) and (53). One-time transformed Fock matrices are constructed once each iteration of the subspace procedure for obtaining the first-order corrections to the orbital-rotation parameters, that is when solving the linear response equations of equations (21) and (24). We have previously shown [17] that the one-time transformed Fock matrices in equations (52)

Computational timings
The spectrum calculations for molecules A-C required the construction of less than 300 one-time transformed Fock matrices to obtain the first-order corrections to the orbital-rotation parameters. For the full method, the E [3] contractions amount to 63%, 61%, and 60% of the computation time for molecules A, B, and C, respectively. For the reduced method that requires only half the number of two-time transformed Fock matrices, the E [3] contractions were done faster than solving the linear-response equation solver and amounted to 47%, 45%, and 44% of the computation time for molecules A, B, and C, respectively, as seen in table 2. Separate timings for the two-electron Coulomb and DFT exchange-correlation parts of the E [3] contractions are presented in table 3. It is seen that the kernel integrations comprise a mere 34%, 36%, 17%, 19%, 21%, and 21% of the total calculation time for molecules A-C with the full and the reduced methods, respectively. Given the algebraic complexity and length of the kernels in quadratic response that are presented in detail in the supplementary information, this suggests that our DFT implementation is highly efficient.
That said, however, the most important message conveyed in this section and specifically in table 2 is the low number of auxiliary Fock matrix constructions that has been reached with the presented tensor-average approach.

Assessment of the reduced form expression
In order to further assess the accuracy of the reduced form expression of the EFISHG observable, we have performed a benchmark investigation involving a set of six molecules. This set includes small-and medium-sized molecules that are of biochemical interest calculated in the one-photon off-resonance region, as seen in figure 3.
The first-one photon resonances, ω 10   molecules in the benchmark set, we have plottedβ(−2ω, ω, ω) in the spectral region ranging up to ω 10 − 2ℏγ, where 2ℏγ = 124 meV corresponds to the full width at half maximum. As expected from equations (59), (60) and (62), (63) and as seen in figure 3, the full and reduced formulations are displaying some discrepancies the closer one gets to a one-photon resonance point. Most clearly this is seen for caffeine, which has the strongest S 1 ← S 0 one-photon resonance point of all the molecules in this benchmark set. Less distinct but visible, there is also a discrepancy seen for p-nitroaniline that has a strong S 2 ← S 0 transition at 4.14 eV. However, based on the results presented in figure 3, it can be safely concluded that the reduced formulation gives quantitatively correct results for the complex first-order hyperpolarizabilityβ(−2ω, ω, ω) in one-photon off-resonance spectral regions.

Summary and conclusions
Computationally tractable expressions for the calculation of resonance-enhanced SHG spectra from the real part of the first-order hyperpolarizability have been derived and implemented using the complex polarization propagator approach to quadratic order in perturbation theory. The present work adopts the density functional theory approximation with a Fock-matrix driven handling of the contractions of secondand third-order generalized electronic Hessian matrices with trial and first-order response vectors.
We have demonstrated a highly efficient algorithm for obtaining the vector components of the quadratic response function for SHG spectra for molecules in an EFISHG experiment where the vector components of the β-tensor is put in computational focus without the explicit reference to individual tensor components. The presented tensor-average algorithm introduces compounded Fock matrices without the introduction of approximations, SHG intensities are determined with a mere 50% of the number of unique Fock matrices required in the tensor-component approach. With the introduction of physically well motivated approximations in one-photon off-resonance regions of the spectrum, the number of Fock matrices in the calculation has been further reduced by a factor of two. The error in the SHG intensities as introduced by these approximations are shown to be negligible in example calculations on a set of three makaluvamines as well as for a wider benchmark study of small-and medium-sized molecules that are of biochemical interest.
The complex polarization propagator approach presented here is perfectly applicable to systems with a high density of states. Combined with the fact that the presented novel reduced-form algorithm largely removes CPU and memory bottlenecks associated with quadratic-response based SHG calculations, we argue that our work will be important so as to enable SHG spectrum simulations for large-scale systems of real technical interest.

Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Taking the vector derivative of the scalar function equation (A.13) we get the vector: Likewise, the second-order vector derivative of the scalar function equation (A.13) we get the rank-two tensor: