Matrix elements between nonorthogonal Slater determinants represent an essential component of many emerging electronic structure methods. However, evaluating nonorthogonal matrix elements is conceptually and computationally harder than their orthogonal counterparts. While several different approaches have been developed, these are predominantly derived from the first-quantized generalized Slater–Condon rules and usually require biorthogonal occupied orbitals to be computed for each matrix element. For coupling terms between nonorthogonal excited configurations, a second-quantized approach such as the nonorthogonal Wick’s theorem is more desirable, but this fails when the two reference determinants have a zero many-body overlap. In this contribution, we derive an entirely generalized extension to the nonorthogonal Wick’s theorem that is applicable to all pairs of determinants with nonorthogonal orbitals. Our approach creates a universal methodology for evaluating any nonorthogonal matrix element and allows Wick’s theorem and the generalized Slater–Condon rules to be unified for the first time. Furthermore, we present a simple well-defined protocol for deriving arbitrary coupling terms between nonorthogonal excited configurations. In the case of overlap and one-body operators, this protocol recovers efficient formulas with reduced scaling, promising significant computational acceleration for methods that rely on such terms.

## I. INTRODUCTION

Matrix elements between nonorthogonal Slater determinants are increasingly common in emerging electronic structure methods. For example, capturing strong correlation using a linear combination of nonorthogonal Slater determinants is a relatively old idea^{1–4} that has seen a renaissance in the past decade.^{5–14} Similarly, nonorthogonal matrix elements arise in projected Hartree–Fock methods,^{15,16} while the combination of geminal-based nonorthogonal functions is an area of ongoing research.^{17} In each method, the nonorthogonality of different determinants can capture strong static correlation effects by breaking and restoring symmetries of the Hamiltonian,^{15} or it can provide a quasi-diabatic representation of dominant electronic configurations.^{5,9} Alternatively, multiple wave functions built from different orbitals arise in orbital-optimized excited states identified through methods such as ΔSCF,^{18–21} excited-state mean-field theory,^{22,23} or the complete active space self-consistent field (SCF).^{24,25} In these cases, orbital optimization can significantly improve predictions of charge transfer excitations, but nonorthogonal matrix elements are required for inter-state coupling terms such as oscillator strengths.

A variety of different approaches have been developed for the efficient evaluation of nonorthogonal matrix elements,^{26–30} which are predominantly derived from Löwdin’s general formula.^{31} The most popular framework in quantum chemistry is the generalized Slater–Condon rules,^{27,32} where biorthogonal occupied orbitals are constructed^{33,34} and a modified form of the Slater–Condon rules^{35} is applied depending on the number of zero-overlap orbital pairs in the biorthogonal basis. This approach is applicable to any pair of determinants but requires the diagonalization of the occupied orbital overlap matrix each time. In contrast, the development of many-body correlation methods using orthogonal determinants has greatly benefited from the second-quantized Wick’s theorem.^{36} While a nonorthogonal variant of Wick’s theorem exists, it is limited to determinants that have a non-zero many-body overlap and is not applicable if there are zero-overlap orbital pairs in the biorthogonal basis.^{37,38} This limitation arises because the Thouless theorem,^{39} used to relate two nonorthogonal determinants via an exponential transformation, breaks down when the two determinants have a zero many-body overlap. As a result, the nonorthogonal Wick’s theorem has seen only limited use in quantum chemistry.^{16,40,41}

Computationally efficient nonorthogonal matrix elements become increasingly important in methods that use orthogonally excited configurations from nonorthogonal reference determinants. For example, including post-NOCI dynamic correlation in methods such as perturbative NOCI-MP2^{42–44} and NOCI-PT2,^{14} or nonorthogonal multireference CI,^{3,16,41,45} requires overlap, one-body, or two-body coupling terms between excitations from nonorthogonal determinants. The number of nonorthogonal matrix elements therefore grows rapidly, and repeated biorthogonalization of the occupied orbitals becomes prohibitively expensive. In principle, the nonorthogonal Wick’s theorem could allow these matrix elements to be evaluated using only biorthogonal reference orbitals, but until now, this requires the reference determinants to have a strictly non-zero overlap.

In this contribution, we derive an entirely generalized nonorthogonal form of Wick’s theorem that applies to any pair of determinants with nonorthogonal orbitals, even if the overall determinants have a zero overlap. This new framework, which we call the “extended nonorthogonal Wick’s theorem,” provides the most general approach for deriving matrix elements using second-quantization, allowing Wick’s theorem and the generalized Slater–Condon rules to be unified for the first time.

One particular advantage of our approach is that it allows all matrix elements between excited configurations from a pair of nonorthogonal determinants to be derived using a single well-defined protocol. While some of the resulting expressions have previously been derived for various bespoke applications (see, e.g., Refs. 6, 42, and 45), these have often relied on the properties of matrix determinants to account for orbital excitations. Instead, we present a unifying theory that can recover all of these results and is automatically applicable in cases where the reference determinants have a zero overlap. Furthermore, we show how evaluating intermediates for a given pair of determinants can reduce the scaling of overlap and one-body coupling terms between excited configurations to $O1$. These particular non-orthogonal matrix elements then become almost as straightforward as the orthogonal Slater–Condon rules or Wick’s theorem, promising considerable acceleration for methods that rely on such terms.

To derive our generalized nonorthogonal matrix elements, we first define our notation in Sec. II. In Sec. III, we extend Thouless theorem^{39} to the case where the two determinants have nonorthogonal orbitals and a zero many-body overlap. Section IV combines this extended Thouless transformation with Wick’s theorem to create a generalized protocol for evaluating nonorthogonal matrix elements using second-quantization. We then illustrate the application of this approach by re-deriving the generalized Slater–Condon rules for one- and two-body operators in Sec. V. Finally, in Sec. VI, we extend our framework to the matrix elements between excited configurations and show how $O1$ scaling can be achieved for overlap and one-body operators.

## II. NOTATION

We will consider matrix elements between the two determinants $\Phi w$ and $\Phi x$. Each determinant is constructed from a bespoke set of molecular orbitals (MOs), represented in terms of the atomic spin-orbital basis functions $\chi \mu $ as

Here, we employ the nonorthogonal tensor notation of Head-Gordon *et al.*^{46} to explicitly keep track of any required overlap matrices. Occupied MOs are indexed as *i*, *j*, *k*, etc., virtual MOs are indexed as *a*, *b*, *c*, etc., and any general MO is indexed as *p*, *q*, *r*, etc. We emphasize that the MOs are orthogonal *within* each Slater determinant but are nonorthogonal *between* the different determinants.

In second-quantization, the *N*-electron determinant $\Phi w$ is defined as

where $\u2212$ is the physical vacuum and the molecular orbital creation operators $bi\u2020x$ satisfy the standard fermionic anticommutation rules.^{47} Using the expansion [Eq. (1)], the MO creation and annihilation operators can be represented in terms of the covariant atomic spin-orbital creation $a\mu \u2020$ and annihilation *a*_{μ} operators as

The covariant atomic spin-orbital operators have only one non-zero anticommutator,^{47}

where *g*_{μν} = ⟨*χ*_{μ}|*χ*_{ν}⟩ defines the corresponding covariant metric tensor (overlap matrix).^{46}

Throughout this paper, we will need to express the atomic spin-orbital creation and annihilation operators in terms of the MO creation and annihilation operators using the inverse of Eq. (3),

To avoid the introduction of overlap matrices throughout our expressions, we will often use the contravariant atomic spin-orbital operators $(a\mu )\u2020$ and *a*^{μ} defined as

where *g*^{μν} is the contravariant metric tensor corresponding to the inverse covariant overlap matrix,^{46} i.e.,

If the AO basis is overcomplete, this contravariant metric tensor becomes the pseudo-inverse of the covariant overlap matrix. Note that the anticommutator of these contravariant atomic spin-orbital operators is

## III. EXTENDED THOULESS TRANSFORMATION

### A. Conventional Thouless transformation

The conventional form of Thouless theorem allows two nonorthogonal determinants to be related by an exponential operator of single excitations as^{39}

To derive the single excitation operator $Z$, the occupied orbitals can be transformed to a biorthogonal basis using Löwdin pairing^{33,34} such that

The virtual–occupied and virtual–virtual blocks of the biorthogonalized overlap matrix become

and the transformed molecular orbital creation and annihilation operators are given as

The single excitation operator in Eq. (9) is then defined as

with the ^{xw}*Z*_{ai} matrix elements given by

A brief derivation of this result can be found in Appendix A.

Unfortunately, this exponential representation relies on the strict nonorthogonality of the two determinants ⟨^{x}Φ$|w$Φ⟩ ≠ 0; in other words, it is not applicable to a pair of determinants that are orthogonal but contain mutually nonorthogonal orbitals. Our first step is therefore a generalization of the Thouless transformation to the case where ⟨^{x}Φ$|w$Φ⟩ = 0.

### B. Introducing zero-overlap orbitals

We begin in the biorthogonal basis identified through Löwdin pairing, with orbital coefficients satisfying Eq. (10). For a general pair of nonorthogonal orbitals, it is possible for orbital pairs to have a zero-overlap in the biorthogonal basis, where $S\u0303ixw=0$. Taking the case with *m* zero overlaps between orbitals *k*_{1}, …, *k*_{m}, we construct “reduced” determinants by removing the electrons in these zero-overlap orbitals to give

These reduced determinants are strictly nonorthogonal with the non-zero reduced overlap defined as

Therefore, the Thouless transformation can now be applied to these reduced determinants to give

Here, we have introduced the reduced single excitation operator $Z\u0303$ that only contains excitations from occupied orbitals with a non-zero overlap as

The full *N*-electron determinants are then related through second-quantization as

Equation (19) can be further simplified by exploiting the commutativity relation $[Z\u0303,b\u0303kx]=0$ to shift the $exp(Z\u0303)$ operator to the far right-hand side, giving

We can then introduce single-electron excitation operators $z^k$ for the zero-overlap orbitals as

where we have exploited the biorthogonality and zero-overlap of the occupied orbitals such that $S\u0303ikxw=0$ for all *i*. The commutativity of these single excitation operators with the $b\u0303kx$ annihilation operators leads to the simplified relationship

Introducing the relationship $z^k=exp(z^k)\u22121$ allows Eq. (22) to be expanded as

Expanding the product of $(exp(z^k)\u22121)$ terms then leads to a sum of exponential transformations where every combination of the zero-overlap single excitation operators $z^k$ is either included or excluded with an appropriate phase factor, giving

Here, we have introduced the compound index $Cn$ to denote a particular combination of *n* out of *m* zero-overlap orbitals, while the superscript notation $Z\u0303Cn$ indicates which particular zero-overlap excitations are included in the corresponding operator, i.e.,

We refer to this transformation in Eq. (24) as the “extended Thouless transformation.” To explicitly illustrate its application, the case of two zero-overlaps in orbital pairs *k*_{1} and *k*_{2} leads to

where $Z\u0303(k1,k2)=Z\u0303+z^k1+z^k2$ and $Z\u0303(k1)=Z\u0303+z^k1$.

## IV. EXTENDED NONORTHOGONAL WICK’S THEOREM

### A. Conventional Wick’s theorem

Efficiently deriving matrix elements using the conventional Wick’s theorem requires the introduction of contractions, defined for two creation or annihilation operators ^{x}*b*_{p} and ^{x}*b*_{q} as

where {^{x}*b*_{p}^{x}*b*_{q}} represents a normal-ordered operator string with respect to the reference Fermi vacuum ⟨^{x}Φ|⋯|^{x}Φ⟩.^{36} The only non-zero contractions between creation and annihilation operators with respect to this symmetric Fermi vacuum are

Through Wick’s theorem, the Fermi vacuum expectation of an operator product is given by the sum over all fully contracted products of operators, e.g.,

### B. Zero-overlap transformed operators

Using the extended Thouless transformation, we can now extend the nonorthogonal Wick’s theorem^{37,38,48} to derive matrix elements between *any* pair of determinants with mutually nonorthogonal orbitals. In what follows, we will consider the contravariant atomic spin-orbital operators (see Sec. II) to avoid large numbers of overlap matrices in our expressions. The matrix elements for general operators expressed in the atomic spin-orbital basis requires the evaluation of terms containing a string of creation and annihilation operators, such as $\u27e8\Phi x|(a\mu )\u2020(a\nu )\u2020\cdots a\sigma a\tau |\Phi w\u27e9$. Applying the extended Thouless transformation leads to the linear combination

To evaluate each constituent matrix element for the combinations $Cn$, we follow the approach described in Refs. 38 and 40 and introduce a similarity-transformed set of spin-orbital creation and annihilation operators as

These operators clearly depend on the particular combination $Cn$ of included zero-overlap single excitation operators. Expanding the similarity transformation as

and similarly for $(a\mu )\u2020$, leads to the explicit forms

An explicit derivation of these relationships can be found in Appendix B. Exploiting the relationship

and the resolution of the identity

then allows the constituent matrix elements within Eq. (30) to be expressed as

### C. The fundamental contraction

The extended Thouless transformation essentially converts the nonorthogonal matrix element with an asymmetric Fermi vacuum ⟨^{x}Φ|⋯$|w$Φ⟩ to a transformed matrix element with respect to symmetric Fermi vacuum ⟨^{x}Φ|⋯|^{x}Φ⟩. Since the transformed operators $d[Cn]\mu $ and $(d[Cn]\mu )\u2020$ are expressed purely in terms of the $b\u0303px$ creation and $b\u0303p\u2020x$ annihilation operators, with respect to the ⟨^{x}Φ|⋯|^{x}Φ⟩ vacuum, their non-zero contractions with respect to ⟨^{x}Φ|⋯|^{x}Φ⟩ can be derived by combining Eqs. (28) and (33) to give

where we have introduced the general notation

These expressions closely resemble the co-density matrices in Ref. 5, and their derivation can be found in Appendix C.

The overall matrix element $\u27e8\Phi x|(a\mu )\u2020(a\nu )\u2020\cdots a\sigma a\tau |\Phi w\u27e9$ requires the derivation of contractions between the *a*^{μ} and $(a\mu )\u2020$ atomic spin-orbital operators rather than the transformed $d[Cn]\mu $ and $(d[Cn]\mu )\u2020$ operators. To show how a general string of operators can be evaluated using the contractions defined in Eqs. (37a) and (37b), we first demonstrate the derivation of the one-body co-density matrix element

Assuming that some form of Wick’s theorem can be derived, we expect this matrix element to be represented by the single contraction

Taking the most general case with *m* zero-overlap orbitals, the extended Thouless transformation leads to

Reversing the order of summation over *k* and $Cn$ yields

The first term in square brackets $\u2211Cn1$ can be recognized as the total number of ways to pick *n* orbitals from the *m* zero-overlap orbitals, given by $mn$. Similarly, the second term in square brackets $\u2211Cn\u0338\u220b\u2009k1$ is the total number of ways to pick *n* orbitals from the *m* − 1 zero-overlap orbitals that remain when orbital *k* is excluded, given by $m\u22121n$. These combinatorial expansions allow Eq. (42) to be expressed as

Here, we note that there are no ways to exclude an orbital *k* when all zero-orbital overlaps are included in the complete combination $Cm$. Exploiting the binomial expansion

then leads to the closed-form expression

The reduced overlap $S\u0303xw$ will be a prefactor for every matrix element between these nonorthogonal determinants. The remaining terms can then be used to define “fundamental contractions” for second-quantization operators with respect to the asymmetric Fermi vacuum ⟨^{x}Φ|⋯$|w$Φ⟩. The form of these contractions depends on the number of zero-overlap orbitals *m*, and we define the first fundamental contraction as

Similarly, the second fundamental contraction can be identified as

Crucially, we emphasize that these fundamental contractions are defined with respect to the asymmetric Fermi vacuum ⟨^{x}Φ|⋯$|w$Φ⟩.

### D. Combining several contractions

Next, we show how the fundamental contractions can be combined to derive matrix elements for longer products of creation and annihilation operators. As an example, consider the two-body reduced co-density matrix element, defined as

Applying Wick’s theorem, this matrix element should be given by the sum of the two contractions,

Note that the second term in this expression carries a phase of −1 from the intersection of the contraction lines, representing the fundamental parity of fermionic operators.^{36} Taking the first contraction as an example, we apply the extended Thouless transformation and the transformed contractions defined in Eqs. (37a) and (37b) to give

Once again, the reduced overlap $S\u0303xw$ appears as an overall prefactor. The order of summation over the *k*_{1}, *k*_{2} indices and $Cn$ can then be swapped to give

The third term in square brackets $\u2211Cn\u0338\u220bk1,k21$ is simply the number of ways to pick *n* orbitals from the *m* − 2 zero-overlap orbitals that remain when orbitals *k*_{1} and *k*_{2} are removed, given by $m\u22122n$. Applying the binomial expansion Eq. (44) in an analogous way to the single contraction leads to the closed form

A similar expression can be derived for the second contraction as

which, analogously with the orthogonal case, will carry a −1 phase factor from the intersection of the contraction lines. Combining Eqs. (52) and (53), with their associated phase factors, yields the full expression for the two-body reduced co-density matrix elements with *m* zero-overlap orbitals as

### E. General rules for constructing matrix elements

To simplify the derivation of even longer operator strings, we note that the two-body matrix elements in Eq. (52) can be factorized into the product of two fundamental contractions with individual *m*_{1} and *m*_{2} values under the constraint *m*_{1} + *m*_{2} = *m*, giving

This factorization can be extended for a general product of contractions with each *m*_{i} restricted to the values 0 or 1 for the overall product to be non-zero. We can therefore define an intuitive approach for extending Wick’s theorem to generalized nonorthogonal matrix elements as follows:

Construct all fully contracted combinations of the operator product with the associated phase factors.

For each term, sum every possible way to distribute

*m*zeros among the contractions such that ∑_{i}*m*_{i}=*m*.For every set of {

*m*_{i}} in each term, construct the relevant contribution as a product of fundamental contractions depending on whether each contraction has*m*_{i}= 0 or 1.Multiply the final combined expression with the reduced overlap $S\u0303xw$.

These rules and contractions therefore allow any matrix element to be evaluated with respect to the asymmetric Fermi vacuum ⟨^{x}Φ|⋯$|w$Φ⟩ for nonorthogonal orbitals, regardless of whether the many-body determinants are orthogonal or not. As a result, our formulation is the most flexible form of Wick’s theorem and reduces to the previous nonorthogonal^{37,38,48} or orthogonal^{36} variants under suitable restrictions on the MO coefficients. We refer to this approach as the “extended nonorthogonal Wick’s theorem.” In Secs. V and VI, we will show how these steps can be applied to recover the generalized Slater–Condon rules^{32} and to derive matrix elements between excited configurations with respect to the nonorthogonal reference determinants.

## V. GENERALIZED SLATER–CONDON RULES

The generalized Slater–Condon rules provide a first-quantized approach for evaluating matrix elements between nonorthogonal determinants. A detailed description of these rules, their derivation, and their application can be found in Ref. 32. However, to the best of our knowledge, the generalized Slater–Condon rules have never previously been derived though a fully second-quantized framework. In this section, we will show how these rules can be recovered using the extended nonorthogonal Wick’s theorem.

### A. One-body operators

Consider first the one-body operator

with the corresponding matrix element

Applying the extended nonorthogonal Wick’s theorem yields only one non-zero contraction to give

Substituting the fundamental contraction [Eq. (46)] and considering the possible values of *m* immediately yield the one-body generalized Slater–Condon rules^{32} in the form presented in Ref. 5 as

Here, we note that ^{xw}*M*^{νμ} = ^{xw}*W*^{νμ} when there are no zero-overlap orbitals, while $P\nu \mu xw=Pk\nu \mu xw$ for one-zero overlap in orbital *k*. The matrices ^{xw}*M* and ^{xw}*P*_{k} correspond, respectively, to the weighted and unweighted co-density matrices discussed in Ref. 5.

### B. Two-body operators

Next, consider a general two-body operator $v^$ defined in second-quantization as

with the two-electron integrals in the spin-orbital basis defined as

The corresponding matrix element is given by

Recognizing the $\u27e8\Phi x|(a\mu )\u2020(a\nu )\u2020a\tau a\sigma |\Phi w\u27e9$ term as the two-body reduced co-density matrix derived in Eq. (54), we immediately recover

Here, note that the *m* = 1 terms each contain two terms with the single zero-overlap in either the first or second contraction. Exploiting the symmetry $v$_{μντσ} = $v$_{νμστ} and the identity

then allows the two-body generalized Slater–Condon rules^{32} to be recovered as

Note that for the *m* = 1 case with one zero-overlap in orbital *k*, we have exploited the identities $Pk\tau \mu xw\u2009Pk\sigma \nu xx=Pk\sigma \mu xw\u2009Pk\tau \nu xx$ and $Pk\tau \mu xw\u2009Pk\sigma \nu xw=Pk\sigma \mu xw\u2009Pk\tau \nu xw$ to introduce the simplification ^{xw}*M* → ^{xw}*W*.

## VI. MATRIX ELEMENTS FOR EXCITED CONFIGURATIONS

While re-deriving the generalized Slater–Condon rules provides an important verification of the extended nonorthogonal Wick’s theorem, it does not provide any computational advantage over the original framework. On the contrary, the primary focus of our new framework involves deriving matrix elements between excited configurations from a pair of nonorthogonal determinants, e.g., $\u27e8\Phi ij\u2026ab\u2026x|\cdots |\Phi kl\u2026cd\u2026w\u27e9$. Terms of this form arise in perturbative corrections to NOCI,^{14,42–44} the NOCI–CIS approach for core excitations,^{8,49} and the evaluation of ⟨*S*^{2}⟩ coupling terms in NOCI expansions.^{6} Furthermore, these nonorthogonal matrix elements will be required to evaluate inter-state coupling elements between orbital-optimized excited-state wave functions identified using excited-state mean-field theory^{22} or state-specific complete active space SCF.^{24}

Until now, evaluating these matrix elements has required the direct application of the generalized Slater–Condon rules to each pair of excitations. This approach leads to significant computational costs associated with the biorthogonalization of the excited determinants, which scales as $ON3$ each time. In this section, we show how the extended nonorthogonal Wick’s theorem allows these matrix elements to be evaluated using only biorthogonalized reference determinants. When large numbers of coupling elements are required, avoiding the biorthogonalization of these excited configurations significantly reduces the overall computational cost. Furthermore, the cost of certain nonorthogonal matrix elements (including all one-body operators) becomes independent of the number of electrons or basis functions if additional intermediate matrix elements are computed for each pair of reference determinants.

In this section, we describe how the extended nonorthogonal Wick’s theorem can be applied to matrix elements between excited configurations. First, we discuss how the fundamental contractions described in Sec. IV C can be modified to an MO-based form that makes excited configurations easier to handle. We then illustrate the derivation of certain overlap, one-body, and two-body matrix elements between nonorthogonal excited configurations. The resulting expressions are entirely generalized for any pair of nonorthogonal reference determinants and can significantly reduce the computational scaling compared to a naïve application of the generalized Slater–Condon rules.

### A. Asymmetric representation

Evaluating matrix elements between two excited determinants will require the evaluation of the asymmetric contractions with respect to the asymmetric Fermi vacuum ⟨^{x}Φ|⋯|^{w}Φ⟩. Expanding the molecular orbital creation and annihilation operators using Eq. (3) yields

Introducing the fundamental contractions for and defined in Eqs. (46) and (47), respectively, then reduces these expressions to different forms depending on the number of zero-overlaps *m*_{k} associated with the contraction,

Here, we have defined the “screened” overlap terms

and

Although the indexing notation used in Eq. (68a) may seem counterintuitive, we find that it helps to keep track of the bra and ket orbital coefficients in the screened overlap terms.

Crucially, the orbital coefficients used to evaluate these contractions do not need to be the same as those used to evaluate the ^{xw}*M* and ^{xw}*P* matrices. This feature is particularly advantageous as the excited configurations can be represented in terms of the original orbital basis while the ^{xw}*M* and ^{xw}*P* matrices are evaluated in the biorthogonal basis. As a result, only the reference determinants need to be biorthogonalized, and the remaining matrix elements are evaluated in terms of these screened overlap terms. Furthermore, the screened overlap elements are themselves one-body matrix elements that can be computed once for a given pair of determinants and stored, before being combined to evaluate more complicated matrix elements. With *N*_{ref} determinants, the total cost of computing these intermediates therefore scales as $ONref2\u2009n3$.

To take full advantage of these asymmetric contractions, the one- and two-body operators can also be represented in terms of one set of molecular orbitals as

where we have defined the transformed matrix elements

and

Evaluating nonorthogonal matrix elements through the extended nonorthogonal Wick’s theorem then proceeds using similar steps outlined in Sec. IV E:

Assemble all fully contracted combinations of the asymmetric operator and excitation operator product and compute the corresponding phase factors.

For each term, sum every possible way to distribute

*m*zeros among the contractions such that ∑_{i}*m*_{i}=*m*.For every set of {

*m*_{i}} in each term, construct the relevant contribution as a product of fundamental contractions defined in Eqs. (67a) and (67b) depending on whether each contraction has*m*_{i}= 0 or 1.Multiply the combined expression with the reduced overlap $S\u0303xw$ of the reference determinants.

Note that the number of zero-overlap orbitals in these expressions corresponds to the biorthogonalized reference determinants, not the excited configurations.

In Secs. IV B–IV E, we will illustrate this process through a series of typical nonorthogonal matrix elements. As we shall see below, evaluating the sum of every combination of zero-overlap indices assigned to each contraction quickly leads to complicated equations. However, we can derive the general structure of a matrix element for the case with no zero-overlap orbitals *m* = 0 and then recover the forms for different values of *m* by distributing the zero-overlap indices over each fundamental contraction. We will therefore focus on the parent equation with *m* = 0, which we refer to as the “canonical form” and denote using the notation ⟨⋯⟩_{0}. Furthermore, any matrix element with more zero-overlap orbitals than the total number of contractions must be strictly zero.

### B. Overlap terms

First, we consider the overlap element between excited determinants. The simplest overlap matrix element involves only a single excitation $\u27e8\Phi x|\Phi iaw\u27e9$. For *m* zero-overlap orbitals in the reference determinants, this single excitation matrix element can be identified using one contraction as

Here, the canonical form is simply

and there is only one way to assign one zero-overlap orbital to the single contraction. Note that the reduced overlap between the *reference* determinants remains a prefactor for the overall matrix elements.

Next, the overlap of two single excitations $\u27e8\Phi iax|\Phi jbw\u27e9$ can be evaluated as

The canonical form for this element is the *m* = 0 term

The *m* = 1 term is recovered by taking the sum of the two different ways to distribute one zero-overlap orbital to the two contractions, while there is only one way to distribute two zero-overlap orbitals for the *m* = 2 term.

As a third example, consider the double excitation overlap $\u27e8\Phi x|\Phi ijabw\u27e9$, which can be evaluated as

In this case, the minus sign arises from the intrinsic −1 phase that results from the intersection of the contraction lines in The relevant canonical form of this matrix element is

### C. One-body operators

We now consider matrix elements for one-body operators of the form given in Eq. (71a). To further simplify the subsequent expressions, we can introduce intermediate matrices that account for partial contraction with the one-body operator. In particular, we introduce the partially contracted intermediate terms

Note that the individual contraction in the *F*_{0} term, indicated by the *X* matrix, may correspond to a zero-overlap orbital pair, as indicated by the notation $F\xaf0$. The ^{ye}[*XFX*]_{rs} (and similar) terms correspond to two contractions and can thus be assigned two zero overlap orbitals. The different possibilities of assigning these zero-overlap contractions is denoted as $[X\xafFX]rsyw$ or $[XFX\xaf]rsyw$ and similarly for terms involving the *Y* contraction. Crucially, these intermediate values also correspond to orbital pairs and can be precomputed once for each pair of reference determinants, leading to a one-off computational cost that scales as $ONref2\u2009n3$. As a result, the summation over the *p*, *q* indices is avoided for the subsequent evaluation of matrix elements between excited determinants, and the computational scaling of these one-body terms is the same as the overlap matrix elements.

To illustrate the application of this approach, we take the simplest one-body matrix element $\u27e8\Phi x|f^|\Phi iaw\u27e9$. For *m* zero-overlap orbitals in the reference determinants, this matrix element can be expanded as

Using the contractions defined in Eqs. (67a) and (67b), and the intermediate terms defined in Eq. (80), the canonical form for *m* = 0 is given as

If *F*_{0}, ^{ww}*X*_{ia}, and ^{ww}[*XFY*]_{ia} are precomputed and stored once for the reference determinants, the subsequent cost of evaluating $\u27e8\Phi x|f^|\Phi iaw\u27e9$ is *independent* of the number of electrons or basis functions and scales as $O1$. From this canonical form, expressions for *m* = 1 and 2 (or higher) can be identified as

Next, consider the coupling of two single excitations $\u27e8\Phi iax|f^|\Phi jbw\u27e9$, corresponding to the fully contracted terms

In this case, the canonical form is given by

Again, the form for *m* zero-overlap orbitals can be recovered as the sum of every way to distribute the zero-overlap contractions over the *X*, *Y*, or *F*_{0} terms in each product. We omit the explicit form of these *m* > 0 expressions to maintain brevity.

Finally, consider the one-body coupling of a reference determinant and a double excitation

The corresponding canonical form for *m* = 0 is

from which expressions for *m* > 0 can be obtained as described above.

### D. Two-body operators

Finally, we consider matrix elements for two-body operators, with the general form given in Eq. (71b). Again, we can define intermediate matrices in the orbital basis that can be pre-computed for a given pair of nonorthogonal determinants. First, we introduce analogs of the Coulomb and exchange matrices, respectively, defined as

Partially contracted intermediate matrices can then be defined as, e.g.,

with the constant value

The one-off computational cost of evaluating each of these intermediate matrices is dominated by the cost of evaluating Eqs. (88) and (90) for each reference determinant, and so the overall scaling is $ONref\u2009n4$.

Like the one-body intermediate matrices, when the zero-overlap orbitals are distributed among the contractions, these zeros may be independently assigned to any of the *X*, *Y*, or (*J* − *K*) terms, as each of these contain one contraction. For example, with *m* = 1, one must consider the sum of three terms given as, e.g.,

where

On the contrary, the *V*_{0} constant term corresponds to two contractions and can be assigned two zero-overlap orbitals. Noting the symmetry ^{x}$v$_{pqrs} = ^{x}$v$_{qpsr}, we denote these possibilities as

Here, the factor of two in Eq. (93a) arises because the zero-overlap orbital can be assigned to either the first or second contraction to give the same result.

To illustrate the application of this approach for excited configurations, we first consider the two-body matrix element

Combining each contraction and exploiting the intermediates in Eqs. (89) and (90) yields the canonical form (*m* = 0) for this matrix element as

Distributing the zero-overlap orbitals among each contraction then leads to explicit expressions for all *m* values as

Next, we consider the two-body coupling of two singly excited determinants as

With a total of 24 possible ways to fully contract the corresponding matrix element, we maintain brevity by directly presenting the canonical form for *m* = 0 as

From here, expressions for the cases with *m* > 0 can be recovered by distributing the zero-overlap orbitals over the contractions associated with the *V*, *X*, (*J* − *K*), or *Y* matrices. While partially contracted intermediate expressions could also be evaluated to avoid the nested summation on the last line in Eq. (98), this will generally incur an unacceptably large storage cost.

Finally, we consider the two-body coupling of a reference determinant and a double excitation

Again, there are a total of 24 possible ways to fully contract the corresponding matrix element, so we advance directly to the canonical form for *m* = 0, given as

We note that this formula has previously been identified in Ref. 45 for *m* = 0, but our approach now allows this result to be extended for *m* > 0.

### E. Illustration of scaling

Through the use of intermediate matrices in Secs. VI A–VI C, the cost of evaluating overlap and one-body matrix elements between excited configurations can be made independent of the number of electrons *N* or basis functions *n*. In both cases, the computational cost of evaluating the intermediate matrices scales as $ONref2\u2009n3$, which is determined by the cost of evaluating the screened overlap terms in Eqs. (68a) and (68b). The subsequent cost of evaluating matrix elements between excited configurations then scales as $O1$. In contrast, the cost of applying the generalized Slater–Condon rules for overlap or one-body operators is dominated by the evaluation of the occupied orbital overlap matrix, giving a scaling of $On2N$ for *every* excited coupling term.

For two-body operators, a similar application of two-body intermediates requires transformations of the two-electron integrals into an asymmetric MO representation for each pair of determinants and for every possible combination of four contractions in the four possible forms given in Eq. (67a) or (67b). These two-body intermediates therefore carry a storage overhead that scales as $O44Nref2\u2009n4$, which will generally be unacceptably large. We imagine that there may be certain applications that require only a subset of these two-electron intermediates, which could then be stored to achieve the subsequent $O1$ scaling. In general, however, the *n*^{4} scaling for two-body matrix elements cannot be overcome, although the overall cost is still reduced by avoiding biorthogonalization of the occupied orbitals in the excited configurations.

To explicitly illustrate the computational speed-up that might be expected using the extended nonorthogonal Wick’s theorem, consider the evaluation of matrix elements between two configuration interaction singles (CIS) wave functions built from different reference Slater determinants,

In particular, consider the evaluation of the overlap

a one-body operator $f^$, e.g., the transition dipole moment,

and a two-body operator $v^$ such as the two-electron repulsion

The relevant overlap, one-body, and two-body nonorthogonal matrix elements are given in Eqs. (76), (84), and (98), respectively. In each case, only one set of intermediates need to be evaluated using the biorthogonalized reference determinants, followed by a total of *N*^{2}*n*^{2} nonorthogonal matrix elements between singly excited configurations (assuming *n* ≫ *N*). Here, we only consider one pairing of the reference determinants.

The computational scaling associated with first evaluating the intermediate terms from the reference determinants and then evaluating all the nonorthogonal coupling terms between the excited configurations is summarized in Table I. For the overlap and one-body terms, the computational scaling is dominated by the evaluation of the intermediates, which scales as $On3$. In comparison, the total cost of biorthogonalizing the occupied orbitals for every pair of excited configurations scales as $ON3n4$. For two-body operators, the cost for both approaches is dominated by the contraction with the two-body integrals for each pair of excited configurations. This leads to a total scaling $ON2n6$ for both cases, but the extended nonorthogonal Wick’s theorem may still be faster as the prefactor is reduced for each element.

. | Intermediate . | Subsequent cost . | ||
---|---|---|---|---|

Coupling term . | Slat.–Con. . | Wick’s . | Slat.–Con. . | Wick’s . |

⟨^{x}Ψ_{CIS}|^{w}Ψ_{CIS}⟩ | N/A | n^{3} | N^{3}n^{4} | N^{2}n^{2} |

$\u27e8\Psi CISx|f^|\Psi CISw\u27e9$ | N/A | n^{3} | N^{3}n^{4} | N^{2}n^{2} |

$\u27e8\Psi CISx|v^|\Psi CISw\u27e9$ | N/A | n^{3} | N^{2}n^{6} | N^{2}n^{6} |

. | Intermediate . | Subsequent cost . | ||
---|---|---|---|---|

Coupling term . | Slat.–Con. . | Wick’s . | Slat.–Con. . | Wick’s . |

⟨^{x}Ψ_{CIS}|^{w}Ψ_{CIS}⟩ | N/A | n^{3} | N^{3}n^{4} | N^{2}n^{2} |

$\u27e8\Psi CISx|f^|\Psi CISw\u27e9$ | N/A | n^{3} | N^{3}n^{4} | N^{2}n^{2} |

$\u27e8\Psi CISx|v^|\Psi CISw\u27e9$ | N/A | n^{3} | N^{2}n^{6} | N^{2}n^{6} |

## VII. CONCLUDING REMARKS

Intuitive derivations and efficient implementations of nonorthogonal matrix elements are becoming increasingly important for the development of nonorthogonal configuration interaction methods and inter-state coupling terms for state-specific excited state wave functions. However, until now, the evaluation of these matrix elements has relied on the generalized Slater–Condon rules, while the second-quantized nonorthogonal Wick’s theorem fails when there are any zero-overlap orbital pairs between the reference determinants. In this work, we have extended the nonorthogonal Wick’s theorem to the case where two determinants have nonorthogonal orbitals but have a zero many-electron overlap. This new theory, which we call the extended nonorthogonal Wick’s theorem, provides the most generalized framework for evaluating matrix elements between two nonorthogonal determinants using second quantization.

Among the primary advantages of our new approach is the ease of deriving and evaluating matrix elements between excited configurations from nonorthogonal reference determinants. To illustrate this feature, we have derived a series of overlap, one-body, and two-body coupling terms between nonorthogonal excited configurations. For overlap terms and one-body operators, these excited nonorthogonal matrix elements can be expressed in terms of one-body intermediate terms that can be precomputed and stored for a given pair of reference determinants. As a result, the subsequent cost of evaluating the coupling of excited configurations scales as $O1$, which is the same as the conventional Slater–Condon rules or Wick’s theorem for orthogonal reference determinants. However, in the current form of the theory, the cost of evaluating two-body coupling terms scales as $On4$, which is the same as applying the generalized Slater–Condon rules for each pair of nonorthogonal excited configurations. From a computational perspective, the extended nonorthogonal Wick’s theorem will therefore provide the greatest acceleration in tasks that require at most one-body coupling terms, such as the evaluation of transition dipole moments for orbital-optimized excited states,^{22–25} or the coupling elements of a one-body reference Hamiltonian in nonorthogonal perturbation theories.^{14,42–44}

Looking forward, we hope that this work will encourage and accelerate future developments in nonorthogonal electronic structure theory by creating a unifying theory for deriving challenging nonorthogonal matrix elements. As part of this vision, we are working on an open-source C++ library for evaluating typical nonorthogonal matrix elements, and we intend to report on this project soon.

## ACKNOWLEDGMENTS

H.G.A.B. was supported by New College, Oxford, through the Astor Junior Research Fellowship. The author is grateful to Rebecca Lloyd for careful proof-reading.

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

### APPENDIX A: THOULESS THEOREM

Thouless’ theorem^{39} allows any Slater determinant to be represented using only single excitations from another determinant, i.e.,

Here, we follow Ref. 40 and provide a brief derivation of this theorem. First, the two sets of second-quantization operators can be related as

where

are sub-blocks of the orbital overlap matrix. Taking a suitable transformation among the orbitals of $\Phi w$ allows these to be further reduced to

The occupied orbitals can then be represented by the transformation

where now

The overall transformation in Eq. (A1) is then given by

Finally, if the occupied orbitals are represented in a biorthogonal basis satisfying Eq. (10), then the *Z* matrix elements are given simply as

### APPENDIX B: SIMILARITY TRANSFORMED OPERATORS

To evaluate common matrix elements using the single excitation operators defined in Eq. (25), we follow Ref. 40 and consider the similarity transformed operators

First, consider the expansion of $(d[Cn]\mu )\u2020$ as

Expanding $(a\mu )\u2020$ using Eq. (6) and inserting the definition of $Z\u0303Cn$ from Eq. (25) lead to

Noting the commutation relation $[b\u0303a\u2020x,b\u0303p\u2020x]=b\u0303a\u2020x\delta ip$ then allows this expression to be simplified as

### APPENDIX C: SIMILARITY TRANSFORMED CONTRACTIONS

To evaluate the contractions and with respect to the Fermi vacuum ⟨^{x}Φ|⋯|^{x}Φ⟩, we first expand the transformed operators in terms of the $b\u0303px$ and $b\u0303p\u2020x$ operators. The only non-zero contractions between these MO operators are and . We can therefore expand the contracted product as

Inserting the definition of ^{xw}*Z*_{ai} and $S\u0303akxw$ from Eqs. (14) and (11), respectively, then leads to the form

Resolving the identity allows us to obtain the additional relationship

which, when inserted into Eq. (C2), leads to

and

Finally, Eq. (37a) is recovered by modifying the separation of summation indices in Eq. (C4) to give the form