Time-Dependent Multiconfigurational Short-Range Density Functional Theory with Generalized Valence Bond Wave Functions

TL;DR

This paper introduces a new computational method for calculating molecular excitation energies and spin couplings using a combination of density functional theory and generalized valence bond wave functions.

Contribution

The novel contribution is the development and implementation of TD-GVB-srDFT, which improves excitation energy accuracy compared to prior methods.

Findings

TD-GVB-srDFT achieves excitation energy accuracy comparable to CAS-srDFT with deviations of 0.2 eV.

The generalized Tamm-Dancoff approximation is essential for accurate triplet excitation calculations.

GVB-srDFT accurately computes spin–spin coupling constants in organic molecules and transition metal complexes.

Abstract

We present a theory and an efficient implementation of
TD-GVB-srDFT,
a time-dependent multiconfigurational range-separated density functional
theory based on generalized valence bond perfect-pairing (GVB-PP)
wave functions. In GVB-srDFT, dynamic correlation effects are incorporated
via range-separation of the Coulomb potential, using tailored Kohn–Sham
functionals of the density. The present implementation builds on our
earlier work on TD-GVB [Hapka et al. J. Chem. Phys.
2022, 156, 174102], which employs
direct Hessian techniques for both wave function optimization and
linear response. We benchmark the performance of TD-GVB-srDFT for
singlet and triplet excitation energies, as well as indirect spin–spin
coupling constants (SSCCs). Compared to the underlying GVB-PP model,
the method significantly improves excitation energies and achieves
accuracy comparable to the complete active space variant, CAS-srDFT,
with mean absolute deviations of 0.2 eV. The use of the generalized
Tamm-Dancoff approximation (gTDA) is mandatory for reliable treatment
of triplet excitations. For organic molecules, SSCCs computed with
GVB-srDFT closely match those from CAS-srDFT and HF-srDFT results,
whereas pure GVB-PP performs markedly worse than CASSCF for all coupling
terms. Both GVB-srDFT and CAS-srDFT accurately reproduce fluorine–metal
couplings in transition metal complexes, provided that gTDA is applied
to singlet contributions.

Figures

MAD (mean absolute deviation) and MSD (mean signed deviation) for singlet excitation energies obtained with DFT, srDFT, GVB-PP, and CASSCF models. The HF-, CAS-, and GVB-srDFT calculations used μ = 0.4 bohr–1.

Deviations of GVB-srLDA (top) and GVB-srPBE (bottom) triplet excitations from the CC3 benchmark , using full linear response (“full E [2], left) and gTDA (right). All GVB-srDFT calculations performed with μ = 0.4 bohr–1. The energy unit is eV.

MAD (mean absolute deviation) and MSD (mean signed deviation) for triplet excitation energies obtained with DFT, srDFT, GVB-PP, and CASSCF models. HF-, CAS-, and GVB-srDFT calculations used the μ = 0.4 bohr–1 value. The darker color indicates that gTDA was applied.

Errors relative to CC3 in the FC contributions to the SSCCs obtained in the “full E[2]” and gTDA approaches, sorted from the most negative to most positive error in CASSCF/CASsrDFT calculations. Lines are drawn to guide to eye.

Calculated SSCCs for the VF6 –1 complex. All the srDFT calculations are with μ = 1.0 bohr–1. The black dashed lines represent the reference coupling ±5 Hz. For each method, four results are reported: full E[2], full gTDA, only gTDA in singlet response, and only gTDA in triplet response.

Tables

1. Summary of Error Statistics (in eV) for Singlet Excitations: Mean Signed Deviation (MSD), Mean Absolute Deviation (MAD), Standard Deviation, and Maximum Absolute Error (MAX)

method	MSD	MAD	std. dev.	MAX
LDA	–0.67	0.72	0.55	2.31
PBE	–0.69	0.74	0.58	2.49
CASSCF	0.48	0.55	0.53	2.46
GVB-PP	0.66	0.72	0.62	2.43
HF-srLDA	0.00	0.23	0.33	1.50
HF-srPBE	0.00	0.24	0.34	1.49
CAS-srLDA	–0.02	0.21	0.28	0.80
CAS-srPBE	–0.04	0.19	0.24	0.76
GVB-srLDA	–0.03	0.19	0.24	0.68
GVB-srPBE	–0.03	0.19	0.24	0.67
GVB-srPBE*	–0.04	0.19	0.24	0.67

2. Summary of Error Statistics (in eV) for Triplet Excitations

method	MSD	MAD	std. dev.	MAX
LDA	–0.37	0.44	0.45	1.37
PBE	–0.50	0.51	0.40	1.64
CASSCF	0.33	0.48	0.66	3.15
GVB-PP	0.62	0.64	0.64	3.30
HF-srLDA	–0.09	0.22	0.25	0.86
HF-srPBE	–0.16	0.27	0.27	0.77
CAS-srLDA	–0.12	0.21	0.22	0.74
CAS-srPBE	–0.19	0.26	0.24	0.73
GVB-srLDA	–0.05	0.19	0.23	0.77
GVB-srLDA*	–0.06	0.19	0.23	0.76
GVB-srPBE	–0.13	0.24	0.26	0.77
GVB-srPBE*	–0.15	0.25	0.26	0.74

Full text

Introduction

1

The geminal wave function ansatz, expressed as an antisymmetrized product of two-electron functions, provides a natural extension of a single-determinant orbital description of electronic structure by explicitly accounting for electron pairing. ?,? Geminal models that satisfy the strong orthogonality condition, known as strictly localized geminal (SLG) wave functions, ?−? ? occupy a special place in quantum chemistry due to their relative simplicity, variational formulation, and ability to partially capture static correlation effects, giving a correct description of single bond dissociation. This work centers on one such approach: the generalized valence bond perfect pairing (GVB-PP) ?,? method, a special case of the more general antisymmetrized product of strongly orthogonal geminals (APSG) ?,? ansatz.

A fundamental limitation of SLG wave functions is their inability to describe correlation between electrons assigned to different geminals. Several strategies have been developed to recover the missing intergeminal correlation energy as a correction to the variational geminal energy. These include second-order perturbation theory, ?−? ? ? ? ? ? as well as coupled cluster (CC)-based methods, such as valence bond CC models, ?−? ? block-correlated formulations, ?,? a linearized multireference CC? theory, and a GVB-based extension of ring CCD. ?,? Another class of methods ?,? employs extended random phase approximation (ERPA)? to recover the missing part of the two-electron reduced density matrix. This approach can be derived in the framework of multireference adiabatic connection theory,? and was later extended using embedding techniques to capture dispersion interactions ?,? and to mitigate symmetry breaking in aromatic systems.?

Next to wave function methods, density functional theory (DFT) has also been employed to improve the accuracy of SLG models. The straightforward approach is to supplement the SLG energy with an additional term obtained from a DFT correlation functional evaluated using the SLG density. This idea was first pursued by Kraka? who combined the LDA functional? with GVB-PP. Two decades later, Cagg and Rassolov? proposed the APSG-DFT method, which relied on the PBE? correlation energy scaled by an empirical, basis-set-dependent factor to mitigate double counting. Odoh et al.? proposed an extension of GVB-PP to separated-pair approximation and combined it with multiconfigurational pair-density functional theory.? Similar hybrid strategies were also explored in geminal theories with lifted strong orthogonality constraint, specifically using pair coupled cluster doubles (pCCD) wave functions.? An alternative way of merging geminals with DFT was explored by Filatov et al.? who developed a variant of the spin-restricted ensemble Kohn–Sham (KS) approach based on the GVB-PP ansatz.

One of the directions in the development of geminal-based methods is their extension to the time-dependent (TD) domain. The advantage of APSG linear response is the ability to access singly excited states and certain types of doubly excited states in multireference systems at low computational cost. Pioneering works in the field by Pernal and co-workers ?,? examined singlet excitation energies for several model molecules both in weakly and strongly correlated regimes. More recently, our study? of singlet excited states on a broader set of small and medium-sized systems confirmed that, due to the absence of dynamic correlation, TD-GVB-PP provides only modest improvements over the Hartree–Fock treatment. Moreover, while the method captures certain doubly excited states, their accuracy is lower than obtained with linear-response complete active space self-consistent field (LR-CASSCF), and some double excitations are entirely absent from the spectrum. To address the missing static correlation between different electron pairs, Li and co-workers ?,? developed an equation-of-motion (EOM) formulation of their block-correlated CC method. When up to three-block correlations are included,? the approach achieves excellent agreement with the density matrix renormalization group (DMRG) in systems that demand large active spaces.

In this work, we present an efficient implementation of a variational DFT-based approach that incorporates dynamic correlation into both the GVB-PP wave function and its linear response to time-dependent perturbations. This is achieved by combining GVB-PP with DFT via range separation of the Coulomb potential. ?,? In range-separated multiconfigurational DFT (MC-srDFT), ?,? the DFT functional is responsible for efficient description of short-range electron correlation, while the long-range wave function captures static correlation without requiring spin-symmetry breaking. Splitting the Coulomb operator not only avoids double counting of correlation in a rigorous manner, but also removes the electron–electron cusp in the long-range wave function. The immediate benefits are shorter wave function expansion and faster convergence with the basis set size. ?,?,?

In geminal theory, range-separated APSG model and its time-dependent extension were first derived and implemented by Pernal et al.? On the example of water and formaldehyde molecules, the authors demonstrated that TD-APSG-srDFT significantly reduces the errors with respect to the pure TD-APSG results in the vicinity of the equilibrium. However, TD-APSG-srDFT fails to describe dissociation curves. In the dissociation limit, the accuracy of the method is inferior to pure TD-APSG due to inadequate description of the ground state by GVB-srDFT originating from self-interaction and static correlation errors in the short-range functional. ?,? Garza et al.? applied range-separation together with the pCCD ansatz. By incorporating the on-top pair density into short-range functionals via auxiliary spin densities, ?,? pCCD-srDFT predicts accurate potential energy curves for single-bond breaking, including symmetric dissociation of water.

Our aim is to provide a systematic assessment of the time-dependent GVB-srDFT model. To enable excitation energy and magnetic response property calculations, we formulate TD-GVB-srDFT for perturbations of both singlet and triplet symmetry. We compare the performance of the method against both the uncorrected GVB-PP linear response and time-dependent complete active space short-range DFT (TD-CAS-srDFT) approach. Singlet and triplet excitation energies are computed for a representative set of small- and medium-sized molecules. In addition, we examine the ability of GVB-srDFT to describe indirect spin–spin coupling constants. Beyond benchmarking on standard data sets, we present a case study of the challenging VF_6_ ^–^ complex.? Special attention is given to the accuracy of GVB-PP relative to CASSCF for triplet-state response, which to the best of our knowledge has not been previously analyzed in the literature.

In Section, we present the details of the implementation, focusing on aspects specific to GVB-PP. The technical details of our calculations are described in Section, followed by a discussion of the results in Section. We conclude in Section with a summary of the key results and observations.

Theory

2

GVB-srDFT Wave Function Optimization

2.1

The central concept in the MC-srDFT theory is the separation of the Coulomb operator into long-range (lr) and short-range (sr) contributions

$$[eqn]$$

satisfying asymptotic conditions

$$[eqn]$$

The separation enables the use of a density functional to describe the short-range electron interactions, while the wave function is employed to handle the long-ranged electron interactions. The ground state energy is obtained by minimizing the range-separated density functional and is expressed as the sum of wave function and functional contributions

$$[eqn]$$

where the long-range Hamiltonian is the sum of kinetic and nuclei-electron operators, T̂ and V̂ ne, respectively, the nuclear repulsion term, V nn, and the long-range electron interaction part

$$[eqn]$$

where V̂ ^lr^ = ∑_ j>i _ v ^lr^(r _ ij _). We follow the typical choice of representing the long-range interaction via an error function

$$[eqn]$$

with μ denoting the range-separation parameter. The E H ^sr^ and E xc ^sr^ terms in eq are the short-range Hartree and exchange-correlation (xc) functionals, respectively. For generalized gradient approximation (GGA) models considered here, the latter depends on charge and spin densities, denoted ρ_ C (r, λ) and ρ S (r, λ), respectively, and their gradients, γ XY (r, λ) = ∇ρ X (r, λ)∇ρ Y (r, λ), where X, Y = C, S; the dependence ξ = {ρ C , ρ S , γ CC , γ SS , γ CS _}. The λ vector collects the variational wave function parameters which we discuss below.

In this work, we consider long-range wave functions built from singlet coupled geminals

$$[eqn]$$

where x̅ = (r, σ) is a combined spatial and spin coordinate, {φ_ p } and {c _ p } are natural orbitals and expansion coefficients, respectively, and P is the geminal index. For each geminal the coefficients are normalized ∀ P _ ∑ p∈P _ c _ p _ ^2^ = 1, and relate to natural spin–orbital occupation numbers, n _ p _ = c _ p _ ^2^, where n _ p _ ∈ [0, 1]. The geminals are strongly orthogonal which means that for each pair it occurs that

$$[eqn]$$

We focus on the special case of GVB-PP wave functions where each geminal is represented by a pair of orbitals

$$[eqn]$$

and the full wave function is an antisymmetrized product of N G = N/2 geminals, where N is the even number of electrons. Note, though, that our implementation can also treat APSG wave functions, if desired.

The variational parameters of a GVB-PP wave function are geminal expansion coefficients and orbital-rotation variables. In ref ? we introduced a geminal-norm conserving parametrization of a GVB-PP wave function

$$[eqn]$$

where the parameters x _ p _ are unconstrained and ∀_ P : ∑ p∈P _(c _ p _ ^0^)^2^ = 1. To describe the orbital rotations we use the antisymmetric real singlet orbital-rotation operator

$$[eqn]$$

where κ is an antisymmetric matrix and Ê _ pq _ is the singlet excitation operator, Ê _ pq _ = â _ pα_ ^†^ â _ qα_ + â _ pβ_ ^†^ â _ qβ_.

The GVB-srDFT wave function optimization is carried out using the restricted-step second-order optimization algorithm, ?−? ? as implemented in Dalton.? The algorithm is based on second-order Taylor expansion of the electronic energy in the λ parameters, around λ = 0

$$[eqn]$$

where g and H are the electronic gradient and electronic Hessian, respectively, and λ = (** x **, κ) gathers the geminal coefficients and orbital rotation variables. Differentiating the GVB-srDFT energy with respect to coefficients and orbital parameters leads to a block structure of both the gradient and Hessian

$$[eqn]$$

The Hessian is obtained as H = P K P, where K is the augmented Hessian for the real-valued GVB-srDFT wave function, and P is the projector operator that removes redundant variables in the coefficient space

$$[eqn]$$

where for each geminal Q we have introduced the redundancy vector p ^ Q ^

$$[eqn]$$ $$[eqn]$$

where I _ r _ denotes the geminal that orbital r belongs to. Similar to the gradient, the Hessian contains both wave function and srDFT components

$$[eqn]$$

The algorithm avoids explicit construction of the full Hessian matrix by directly computing contributions to the Hessian using linear transformations with trial vectors b _ n _

$$[eqn]$$

The trial vectors are chosen to obey P b _ n _ = b _ n _, to facilitate handling of the redundancies.

The GVB-PP long-range contributions to the gradient and σ _ n _ elements are identical to the regular, i.e., full-range, terms described in ref ?, except that the full-range two-electron integrals are replaced with the long-range ones. The short-range exchange-correlation contributions to the gradient and Hessian are obtained analogously to the general MC-srDFT framework.? Below we summarize the key points relevant to our GVB-srDFT implementation.

The orbital contribution to the sr exchange-correlation gradient depends only on the reduced one-electron charge density matrix [c.f., eqs 23–25 in ref ?] which for GVB-srDFT wave function in the NO representation takes the diagonal form

$$[eqn]$$

The geminal contribution to the gradient is the same as in the GVB-PP case,? except it is evaluated using long-range two-electron integrals and a one-electron Hamiltonian modified by the addition of the short-range xc potential, which reads

$$[eqn]$$

where e xc ^sr^(ξ) is the short-range xc energy density, and Ω_ pq (r) = φ p _ ^*^(r)­φ_ q (r). Contributions to the gradient from the Hartree term E H ^sr^[ρ C _] are analogous, with the short-range Hartree potential in place of the short-range xc potential.

Evaluation of the srDFT components of the Hessian, $[eqn]$ , involves two different types of terms

$$[eqn]$$

which follows from the nonlinearity of the srDFT functional. The direct Hessian technique, eq, requires contracting eq with the geminal and orbital trial vectors

$$[eqn]$$

In the implementation, we restrict ourselves to closed-shell singlet GVB-srDFT reference wave functions, cf. eq, so that derivatives with respect to the spin density variables are not considered. Equation can be cast in a convenient form using transition density matrices, $[eqn]$ , of two types: the one-index transformed ?,? density matrix for orbital trial vectors, which depends solely on the density matrix

$$[eqn]$$

and the transition density matrix for geminal trial vectors, which for the specific case of a GVB wave function takes the form

$$[eqn]$$

where we have introduced state vectors for the current geminal trial vector

$$[eqn]$$

These quantities are used to introduce effective potential operators which take the form

$$[eqn]$$

where $[eqn]$ are matrix elements of the charge component of the short-range exchange-correlation kernel [see, e.g., ref ? for standard GGA form and eq 40 in ref ? for srDFT]. With these effective potentials, one may define all necessary Hessian σ-vectors for wave function optimization, cf. Equation 42 in ref ?. The geminal part of σ-vectors before projection takes a simple form

$$[eqn]$$ $$[eqn]$$

where Ṽ srxc ^ C ^ is the one-index transformed ?,? form of the gradient operator in eq, given as

$$[eqn]$$

The wave function parametrization, eq, requires that the σ^ ′ ^-vectors defined in eqs and ? have to be projected onto the x-representation using the P projector introduced in eq: σ = P σ ^ ′ ^. Expressions for the orbital part of the Hessian σ-vectors are identical as in the MC-srDFT case [cf. eqs 42b and 42d in ref ?], except for the evaluation of density matrices and transition density matrices, cf. eqs and ?, respectively.

Although our primary focus is on the linear response models, it is worth noticing that for the standard value of μ = 0.4 the GVB-srDFT optimization yields geminals with a more diffuse character compared to GVB-PP. In particular, the π orbitals of aromatic systems may become delocalized, as illustrated for the paradigmatic benzene molecule in Figures S1–S5 of the Supporting Information.

Linear Response for GVB-srDFT

2.2

The GVB-srDFT linear response equations can be written in a general matrix form ?,?

$$[eqn]$$

where X = C or X = S for singlet and triplet properties, respectively, and V ^[1],X ^ denotes external perturbation vector. Excitation energies correspond to poles of the response functions, i.e., they follow by setting the external perturbation to zero. The nonperturbed problem can be cast in a matrix form

$$[eqn]$$

where Λ_ i _ ^ X ^ denotes the eigenvector, and the ω_ i _ eigenvalues are approximations to excitation energies.

In MC-srDFT, the Hessian separates as a sum of long-range and short-range contributions, E ^[2]^ = E ^[2],lr^ + E ^[2],sr^. The long-range contributions are identical to those given in ref ?, with the lr two-electron integrals instead of the full-range ones. Expressions for the srDFT contributions to the Hessian follow from eqs and ?. In practice, linear response equations are solved using direct iterative algorithms described in ref ?.

To handle both real and imaginary orbital rotations for the solutions of response equations, a more general parametrization of orbital rotation operators is required?

$$[eqn]$$ $$[eqn]$$

where T̂ _ pq _ = â _ pα_ ^†^ â _ qα_ – â _ pβ_ ^†^ â _ qβ_ is the triplet one-electron operator. With these general operators, the one-index transformed matrices take the form

$$[eqn]$$

for both X = C and X = S. Note that eqs and ? are recovered for κ_ pq _ = −κ_ qp _ and b _ pq _ ^ o ^ = −b _ qp _ ^ o ^, i.e., when restricting κ̂ to real orbital rotations.

Since triplet response equations are solved here assuming a real singlet reference wave function, all quantities which are of an odd order dependency on the spin density variable S do not contribute to the srDFT part of the Hessian, as defined in eq [see also eq 30 in ref ?]:

$$[eqn]$$

The srDFT components of the Hessian for triplet response are conveniently expressed via effective triplet potential operators

$$[eqn]$$

where $[eqn]$ are matrix elements of the spin-component of the short-range xc kernel,? and the spin transition density matrices $[eqn]$ are defined with respect to the spin-density matrix D _ pq _ ^ S ^ = ⟨Ψ^GVB^|T̂ _ pq _|Ψ^GVB^⟩. Note that for a singlet-reference wave function D ^ S ^ = 0, but D ^ S,(1λ)^≠0. Since our implementation of the GVB-PP wave function does not include triplet components, i.e., geminals in eq are singlet coupled, the geminal part of the effective triplet potential is equal to zero, V srxc ^ S, [1g]^ = 0. In consequence, the only nonvanishing contributions to the Hessian σ-vectors come from the orbital trial vectors. These read

$$[eqn]$$ $$[eqn]$$

where H̃̂ ^lr^ is the one-index transformed lr Hamiltonian, and V̂ srxc ^ S, [1o]^ = ∑_ pq _ V̂ _srxc, pq _ ^ S, [1o]^ T̂ _ pq _. Note that eq does not contain contribution from the one-index transformed sr-xc gradient, since $[eqn]$ . Equations and ? are computed as in the MCSCF case, see, e.g., eq 53 in ref ?. For completeness, in the Supporting Information we also provide explicit expressions for contributions to the σ-vectors in the singlet linear response equations.

Indirect Nuclear Spin–Spin Coupling

Constants

2.3

The nonrelativistic indirect spin–spin coupling tensor describes the correction to the interaction between naked nuclear dipoles due to the presence of the electron cloud. It can be represented as a sum of four contributions ?,?

$$[eqn]$$

where the diamagnetic spin–orbit (DSO) and paramagnetic spin–orbit (PSO) terms arise from the coupling of nuclear magnetic moments to the orbital motion of the electrons, whereas the spin-dipole (SD) and Fermi-contact (FC) terms result from the coupling of the nuclear magnetic moments to the spin of the electron. The DSO term is evaluated as an expectation value. The remaining contributions are recovered in second-order perturbation theory from linear response functions at zero frequency: the PSO term requires solving singlet response equations; both FC and SD hyperfine tensors follow from triplet linear response functions. For the specific form of the corresponding property gradient vectors, see, e.g., ref ?.

From the isotropic part of the coupling tensor, $[eqn]$ , one recovers the indirect spin–spin coupling constant, which we focus on in this work. We have integrated our GVB-srDFT linear response codes with the SSCC implementation available in Dalton, originally developed for MCSCF wave functions.?

Computational Details

3

The second-order GVB-srDFT wave function optimization and TD-GVB-sr-DFT linear response equations were implemented in a development version of the Dalton program. ?,? For all systems, initial geminals were generated using our black-box, two-step scheme, described in ref ?. The starting strongly occupied orbitals were obtained by localizing relevant SCF orbitals via the Foster-Boys? scheme. GVB-srDFT calculations were performed with short-range LDA (srLDA)? and short-range PBE (srPBE)? exchange-correlation functionals. The computational cost of GVB-srDFT wave function optimization is identical to GVB-PP scaling for the long-range GVB-PP part, plus DFT scaling for the srDFT part, as the total energy is the sum of two contributions. Analogously, the scaling of TD-GVB-srDFT response equations corresponds to the combined scaling of the TD-GVB-PP and TD-DFT parts.

The data set of reference excitation energies is identical to that used in ref ?, comprising 33 molecules with geometries taken either from the Thiel’s set ?,? (16 molecules) or from the set of Loos et al.? (17 molecules). The data set includes 110 singlet and 110 triplet excitation entries, with reference values taken from CC3 calculations as reported in refs ? and ?. The LDA, PBE, CASSCF and CAS-srDFT data for triplet excitations were taken from ref ? (for the active space selection, see the Supporting Information of ref ?). For singlet excitations, we employed the TZVP basis set? for systems of the Thiel’s data set and the aug-cc-pVTZ basis set? for the remaining molecules. For all triplet excitations, the aug-cc-pVTZ basis set was selected. In our previous work, we established that GVB-PP fails to recover certain singlet doubly excited states from the Thiel setspecifically, the 2^1^ A g states of butadiene, hexatriene and octatetraene, as well as the 2^1^ A 1 state of cyclopentadiene.? Since these will also be absent at the GVB-srDFT level, we exclude them from the benchmark set across all methods to ensure a consistent comparison.

All GVB-srDFT calculations were performed employing the full GVB-PP valence geminal space. In the Supporting Information, we also provide excitation energies obtained with a smaller active space selected based on MP2-srDFT natural occupation numbers, analogously to the CAS-srDFT active space defined in ref ? [see Table S1 for the GVB-srDFT­(μ = 0.4) active space composition].

The data set of reference SSCCs by Faber et al.? analyzed in Section comprises 45 data points for 13 small single-reference molecules computed at the CC3/aug-ccJ-pVTZ level of theory. The KS-DFT, CASSCF and CAS-srDFT data for SSCCs were taken from ref ?. All GVB-PP and GVB-srDFT calculations of SSCCs were performed in the aug-ccJ-pVTZ basis set.?

Results

4

Singlet and Triplet Excitations

4.1

Table and Figure present error statistics for singlet excitations. As reported in ref ?, GVB-PP systematically underperforms compared to CASSCF, with both mean absolute deviations (MADs) and mean signed deviations (MSDs) larger by about 0.2 eV. When dynamic correlation is included via a short-range DFT functional, the MAD decreases from 0.7 eV for GVB-PP to 0.2 eV for GVB-srDFT, matching the accuracy of CAS-srDFT. The srLDA and srPBE models perform essentially the same when combined with either GVB-PP or HF; for CASSCF, srPBE offers slightly better accuracy. Since molecules in the set do not feature strong correlation effects, the results for HF-, GVB- and CAS-srDFT are generally similar. Still, both GVB-srDFT and CAS-srDFT reduce the error spread relative to HF-srDFT, with the standard deviation dropping from 0.3 to 0.2 eV. On average, GVB-srDFT exhibits identical mean errors as CAS-srDFT. However, for states with a more pronounced doubly excited character, as indicated by the T 1 diagnostic, GVB-srDFT resembles HF-srDFT and is 0.3–0.4 eV less accurate than CAS-srDFT. In our data set, this occurs for the 2^1^ A 1 state of furan and the 1^1^ B 2u state of s-tetrazine.

1: Summary of Error Statistics (in eV) for Singlet Excitations: Mean Signed Deviation (MSD), Mean Absolute Deviation (MAD), Standard Deviation, and Maximum Absolute Error (MAX)

MAD (mean absolute deviation) and MSD (mean signed deviation) for singlet excitation energies obtained with DFT, srDFT, GVB-PP, and CASSCF models. The HF-, CAS-, and GVB-srDFT calculations used μ = 0.4 bohr–1.

In Figures and ?, we compare error statistics for triplet excitations obtained using the full response (“full E[2]”) and the Tamm-Dancoff approximation, using CC3 results as the benchmark. The gTDA treatment proves beneficial for all methods. The largest improvement occurs for HF- and GVB-srDFT, where gTDA reduces the MAD by 0.2 eV. The effect is most pronounced for excitations below 3 eV, which are expected to suffer from near-triplet instabilities due to the presence of small elements in the Hessian matrix. Notably, the accuracy of GVB-srDFT triplet excitation energies using gTDA remains robust in the entire range studied, up to 10 eV (see Figure S6 in the Supporting Information).

Deviations of GVB-srLDA (top) and GVB-srPBE (bottom) triplet excitations from the CC3 benchmark , using full linear response (“full E [2], left) and gTDA (right). All GVB-srDFT calculations performed with μ = 0.4 bohr–1. The energy unit is eV.

MAD (mean absolute deviation) and MSD (mean signed deviation) for triplet excitation energies obtained with DFT, srDFT, GVB-PP, and CASSCF models. HF-, CAS-, and GVB-srDFT calculations used the μ = 0.4 bohr–1 value. The darker color indicates that gTDA was applied.

Similar to the singlet response, GVB-PP systematically overestimates triplet excitations relative to the CC3 reference and is less accurate than CASSCF, with MAD values of 0.6 and 0.5 eV, respectively (see Table). Pure LDA and PBE underestimate the excitation energies to a comparable extent, the respective MSDs amounting to −0.4 and −0.6 eV. In MC-srDFT, the srLDA functional is slightly more accurate compared to srPBE, regardless of the long-range wave function. Range-separation effectively remedies the deficiencies of the GVB-PP model: GVB-srDFT performs on a par with CAS-srDFT, achieving MAD of 0.2 eV and a standard deviation of similar magnitude (see Table). Overall, GVB-srLDA and CAS-srLDA are the two best-performing methods. The 3^3^ B 2 state of pyridine is the outlier in the set, as both GVB- and CAS-srLDA overestimate the CC3 value by as much as 0.7–0.8 eV. Here, GVB-srLDA­(gTDA) and CAS-srLDA­(gTDA) excitation energies of 8.06 and 8.03 eV, respectively, are much closer to the CASPT2 result? in the same aug-cc-pVTZ basis (7.6 eV) than to the CC3 reference (7.29 eV).

2: Summary of Error Statistics (in eV) for Triplet Excitations

A detailed inspection reveals that the largest deviations between GVB-srDFT and CAS-srDFT triplet excitations, ranging from 0.3 to 0.5 eV, occur for five systems characterized by increased static correlation:? benzene (the 1^3^ B 1u and 1^3^ B 2u states), naphthalene (the 2^3^ B 1g and 3^3^ A g states), all-E-octatetraene (the 1^3^ B u state), pyridine (the 3^3^ A 1 state), and s-tetrazine (the 1^3^ B 1u state). For these excitations, GVB-srDFT is closer to HF-srDFT. Considering the subset of 36 triplet excitations associated with these molecules, the differences between the srDFT approaches are reflected mainly in the error spread: GVB-srDFT improves over HF-srDFT, but remains slightly less accurate than CAS-srDFT, with standard deviations of 0.26, 0.30, and 0.22 eV, respectively (see Table S2 in the Supporting Information). For singlet excitations in the same subset, the differences between multi- and single-reference srDFT methods are even more pronounced. At the srPBE level, both GVB- and CAS-srPBE reduce the MAD by roughly 0.2 eV relative to HF-srPBE, while the standard deviation decreases from 0.5 eV (HF-srPBE) to 0.2 and 0.3 eV for GVB-srPBE and CAS-srPBE, respectively.

In general, MC-srDFT methods require smaller active spaces than their μ = ∞ limit, provided the range-separation parameter is sufficiently small to allow short-range correlation effects to be captured by the srDFT functional.? For CAS-srDFT, the value of μ = 0.4 has been shown to meet this criterion. ?,?,? To verify if the same choice of μ allows for smaller active spaces for geminal-based wave functions, we compared excitation energies from GVB-srDFT (μ = 0.4) calculations performed either with the full-valence GVB-PP active space or with reduced number of geminals. The reduced space was selected based on MP2-srDFT occupation numbers. The differences between individual excitation energies from the two active spaces do not exceed 0.1 eV, and the corresponding error statistics are virtually identical (see GVB-srDFT* entries in Tables and ?, as well as Table S3 in the Supporting Information). Thus, using the full-valence active space offers no clear advantage, and restricting geminals to the strongly correlated orbitals is recommended, consistent with the typical approach in CAS-srDFT.

SSCCs for Single-Reference Molecules

4.2

In this section, we assess the performance of GVB-PP and GVB-srDFT methods in computing indirect spin–spin coupling constants. Since our reference data set consists of small molecules of a single-reference character, we report errors relative to the CC3 benchmark values.

Overall, couplings from GVB-PP are less accurate than the CASSCF ones, with largest deviations between the two models occurring for couplings involving fluorine. Figure (top panels) compares the errors of GVB-PP and CASSCF for the FC contributions to the coupling tensor (results for PSO and SD terms are given in Supporting Information, see Figures S9 and S10). The FC contributions are of a particularly poor quality compared to the remaining contributions computed at the GVB-PP level: GVB-PP deviates from CASSCF not only for fluorine, but also for J CC (FCCH, FCCF) and ^1^ J HC (HCCH, HFCO, FCCH) couplings. Invoking the gTDA approximation worsens the agreement between GVB-PP and CASSCF even further. A boxplot representation of errors reveals that gTDA reduces the GVB-PP error for the triplet SD term, yet the quality of the FC and PSO contributions deteriorates (see Figure S11 in the Supporting Information). In fact, the failure of the gTDA model for the FC contribution persist across all methods, including KS-DFT and MC-srDFT approaches, consistently with earlier findings in the literature. ?,?

Errors relative to CC3 in the FC contributions to the SSCCs obtained in the “full E[2]” and gTDA approaches, sorted from the most negative to most positive error in CASSCF/CASsrDFT calculations. Lines are drawn to guide to eye.

Adding dynamic correlation via the short-range density functional significantly improves the accuracy of GVB-based couplings. As shown in the bottom panels of Figure, the GVB-srLDA and CAS-srLDA results are nearly identical. In the “full E^[2]^” approach, the deviations do not exceed 2 Hz. When gTDA is employed, the differences between CAS-srLDA and GVB-srLDA become more pronounced: for three couplings involving fluorine the methods deviate by 10–20 Hz [^1^ J NF (FNO) and ^1^ J OF (OHF, OF_2_)]. In these cases, GVB-srLDA is closer to HF-srLDA. Overall, gTDA improves SD and PSO terms irrespective of the long-range wave function in the MC-srDFT model. For FC components the picture is different: gTDA is beneficial only for FC couplings involving fluorine, while the remaining types are more accurate when gTDA is not applied (see also Figures S11 and S12 in the Supporting Information). Irrespective of gTDA, the most problematic FC contributions are of the ^2^ J FF and ^1^ J CF types; their exclusion markedly improves MC-srDFT accuracy (cf. Figure S12 in the Supporting Information). Finally, since srLDA and srPBE perform similarly, the cheaper srLDA model is recommended.

Transition Metal MF6

–z Complexes

4.3

In ref ?, Kjellgren and Jensen analyzed SSCCs of the ^1^ J MF and ^1^ J MC types in transition metal complexes stabilized by either carbon monoxide [M­(CO)_ y _ ^–z ^] or fluorine [MF_6_ ^–z ^] ligands. In this work, we focus on the most challenging system from that test set: VF_6_ ^–^. This is the only complex in the group where both FC and PSO terms contribute to the coupling by comparable magnitude. Since these terms have opposite signs, getting the correct cancellation between them becomes particularly challenging. As reported in ref ?, satisfactory accuracy for carbonyl complexes is already achieved at the Kohn–Sham DFT level, and introducing range separation brings little or no improvement. Therefore, these complexes are not included in the present study. Results for the remaining fluorine-containing systems from ref ? (ScF_6_ ^–3^ and TiF_6_ ^–2^) are provided in the Supporting Information. Reference data are taken from ref ?. Each MC-srDFT method is applied with four variants for solving the linear response problem: full response, gTDA, or gTDA applied to only the singlet or only the triplet part. Following ref ?, the range separation parameter μ is set to 1.0 bohr^–1^ for all transition-metal systems. All CAS-srDFT results were obtained using a CAS­(10,10) active space, which includes the full d-shell. Similar to the single-reference systems discussed in the previous section, srLDA and srPBE functionals yield SSCCs of comparable quality (see also Table S4 in the Supporting Information), and only srLDA results are discussed below.

In Figure we show the ^1^ J VF couplings in the VF_6_ ^–^ complex. The application of gTDA has a large impact on the couplings, with a significant effect both in the triplet (FC, SD) and singlet (PSO) contributions. This stays in contrast to ScF_6_ ^–3^ and TiF_6_ ^–2^, where the FC term dominates, so that invoking gTDA in the singlet response equations has little effect (see Figures S14 and S15 in the Supporting Information). As shown in Figure, MC-srDFT methods clearly outperform both DFT and pure wave function (GVB-PP and CASSCF) approaches. In agreement with ref ?, only one variant of MC-srDFT response provides accurate results, i.e., using gTDA for singlet components and solving the full response for the triplet ones. This approach yields similar accuracy regardless of the long-range wave function type, when combined with the srLDA functional. The relatively small influence of nondynamical correlation is consistent with the large natural occupation numbers, all exceeding 1.98 at the GVB-PP level. For ScF_6_ ^–3^ and TiF_6_ ^–2^, MC-srDFT again clearly outperforms DFT and markedly improves upon the pure CASSCF and GVB-PP methods, with best results obtained using either full response or gTDA in the singlet part (see Figures S14 and S15 in the Supporting Information).

Calculated SSCCs for the VF6 –1 complex. All the srDFT calculations are with μ = 1.0 bohr–1. The black dashed lines represent the reference coupling ±5 Hz. For each method, four results are reported: full E[2], full gTDA, only gTDA in singlet response, and only gTDA in triplet response.

Overall, the largest improvement of MC-srDFT over DFT comes from mitigating self-interaction error through range separation. The effect of accounting for multireference character is considerably smaller, but still noticeable: GVB-srDFT and CAS-srDFT offer similar accuracy and are systematically better than HF-srDFT.

Conclusions

5

In this work we presented a method for treating dynamic correlation in strongly orthogonal geminal wave function theory, applicable to both ground and excited states. The approach is formulated in the framework of multiconfiguration range-separated DFT theory ?,? and its time-dependent extension. ?,? It combines long-range wave function of the GVB-PP type with short-range functionals from the LDA (srLDA?) and GGA (srPBE?) rungs of the Jacob’s ladder. GVB-srDFT wave function optimization is performed employing the restricted step second-order algorithm previously developed for CAS-srDFT? as an extension of the original MCSCF implementation. ?,?,? For the initial geminal guess we use the black-box protocol first introduced for GVB-PP and APSG wave functions.? Excited states of singlet and triplet character are accessed via linear-response GVB-srDFT theory, derived for singlet reference. Linear response equations are solved using a direct iterative algorithm ?,? as implemented in the Dalton ?,? program.

We focused on determining whether GVB-srDFT adequately describes excited states through linear response. To this end, we assessed its accuracy for both singlet and triplet excitation energies on a diverse set of molecules from benchmark data sets of Loos et al.? and Schreiber et al.? for which CC3 reference results are available. These small and medium-sized systems also allow direct comparison with a more accurate member of the MC-srDFT family, namely the CAS-srDFT method. GVB-srDFT offers a substantial improvement over the pure GVB-PP model: for singlet excitation energies, the MAD is reduced from 0.7 eV (GVB-PP) to 0.2 eV (GVB-srDFT). For triplet excitations, the corresponding MADs amount to 0.6 and 0.2 eV, respectively, when computed using the generalized Tamm-Dancoff approximation. When gTDA is not invoked, the performance deteriorates, with MADs increasing to 0.4 and 0.5 eV at the GVB-srLDA and GVB-srPBE levels, respectively. The more expensive srPBE functional offers no advantage over srLDA; in fact, the latter is systematically more accurate for triplet response and is the recommended choice.?

Importantly, GVB-srDFT closely matches the accuracy of CAS-srDFT, which is similar to that of CC2 and ADC(2) methods. ?,? On a subset of systems with enhanced static correlation,? GVB-srDFT and CAS-srDFT perform on par for singlet excitations, both reducing the MAD by 0.2 eV relative to HF-srDFT. However, for triplet excitations in the same subset, GVB-srDFT remains closer to HF-srDFT, with a larger spread of errors than CAS-srDFT, as indicated by standard deviations of 0.3 eV (GVB-srDFT) and 0.2 eV (CAS-srDFT). A key limitation of both GVB-PP and GVB-srDFT is their inability to describe certain doubly excited states (see the discussion in refs ?, ? ). This leads to larger errors for excitations with significant contributions from doubly excited determinantsfor example, the 2 ^1^ A 1 state of furan, where GVB-srLDA is approximately 0.3 eV less accurate than CAS-srLDA.

Indirect spin–spin couplings serve as a sensitive probe of the quality of linear response functions, as they depend on both singlet and triplet response components. We verified the performance of GVB-srDFT for SSCCs by comparison with CC3 benchmarks using a data set of small single-reference molecules,? as well as the challenging transition metal complex VF_6_ ^–^.? In the molecular data set, couplings computed with pure GVB-PP are clearly less accurate than those from CASSCF. In contrast, GVB-srDFT and CAS-srDFT results match closely. Since nondynamical correlation effects are negligible in this set, HF-srDFT is as accurate as both MC-srDFT models. For couplings not involving fluorine, range separated approaches improve upon both the pure KS-DFT and wave function limits, offering a quantitative description. However, for the most problematic cases, that is, couplings involving fluorine, range separation alone does not compensate for the deficiencies of the short-range functionals. In consequence, none of the MC-srDFT variants proves satisfactory for couplings with fluorine, with ^1^ J FF, ^1^ J OF, and ^1^ J CF presenting the greatest challenge.

Our study supports the recommendation that, within the MC-srDFT framework, gTDA should be applied only to the singlet response when calculating SSCCs. As observed in ref ?, gTDA systematically improves the singlet PSO contributions, but is detrimental to FC triplet contributions. This indicates that low-lying triplet excitations, which are better described with gTDA, are not critical for getting the correct FC terms. While gTDA also benefits the SD terms, these are generally small in magnitude.

For the vanadium–fluorine coupling in the VF_6_ ^–^ complex, use of TDA/gTDA is crucial to achieving the correct cancellation between the negative PSO and positive FC components. Across three isoelectronic MF_6_ ^–z ^ complexes (ScF_6_ ^–3^, TiF_6_ ^–2^, and VF_6_ ^–^), GVB-srLDA consistently outperforms HF-srLDA and is comparable to, or slightly better than, CAS-srLDA.

One of the advantages of MC-srDFT methods is that they provide accurate results with smaller active spaces than required in full wave function approaches. We confirmed that this also applies to GVB-srDFT: singlet and triplet excitation energies computed with a reduced number of geminals deviate by no more than 0.1 eV from those obtained using a full-valence active space, and the associated error statistics are virtually identical. To define the reduced active space, we selected geminals based on MP2-srDFT natural occupation numbers. Finally, we recommend using the standard value of μ = 0.4 bohr^–1^ for GVB-srDFT, as it provides excellent accuracy for excitation energies and SSCCs in organic molecules. In contrast, reliable prediction of SSCCs in transition metal complexes requires a larger value, μ = 1.0 bohr^–1^, likely due to medium-range nondynamical correlation effects involving the metal d-shell.?

To summarize, our results demonstrate that TD-GVB-srDFT provides a reliable description of excited states. The method is especially promising for systems in which the active space required for CAS-srDFT becomes prohibitively large. Future extensions employing on-top pair densities in the short-range functionals should enable a correct description of spin degeneracies and significant reduction of self-interaction error. ?,? Moreover, GVB-srDFT linear response can be adapted to describe noncovalent interactions in systems with significant nondynamical correlation, within the framework of multiconfigurational symmetry-adapted perturbation theory, SAPT­(MC). ?,? Ongoing work in our groups is focused on developing this approach.

Supplementary Material