Generalized nonorthogonal matrix elements: Unifying Wick’s theorem and the Slater–Condon rules

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.¹⁷ In each method, the nonorthogonality of different determinants can capture strong static correlation effects by breaking and restoring symmetries of the Hamiltonian,¹⁵ 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.³¹ 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³⁵ 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.³⁶ 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,³⁹ 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,¹⁴ 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³⁹ 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.

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

Here, we employ the nonorthogonal tensor notation of Head-Gordon et al.⁴⁶ 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 $Φw$ is defined as

where $−$ is the physical vacuum and the molecular orbital creation operators $bi†x$ satisfy the standard fermionic anticommutation rules.⁴⁷ Using the expansion [Eq. (1)], the MO creation and annihilation operators can be represented in terms of the covariant atomic spin-orbital creation $aμ†$ and annihilation a_μ operators as

The covariant atomic spin-orbital operators have only one non-zero anticommutator,⁴⁷ 

where g_μν = ⟨χ_μ|χ_ν⟩ defines the corresponding covariant metric tensor (overlap matrix).⁴⁶ 

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μ)†$ and a^μ defined as

where g^μν is the contravariant metric tensor corresponding to the inverse covariant overlap matrix,⁴⁶ 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

The conventional form of Thouless theorem allows two nonorthogonal determinants to be related by an exponential operator of single excitations as³⁹ 

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 ^xwZ_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.

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̃ixw=0$⁠. Taking the case with m zero overlaps between orbitals k₁, …, 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̃$ 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̃,b̃kx]=0$ to shift the $exp(Z̃)$ 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̃ikxw=0$ for all i. The commutativity of these single excitation operators with the $b̃kx$ annihilation operators leads to the simplified relationship

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

Expanding the product of $(exp(z^k)−1)$ 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̃Cn$ 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₁ and k₂ leads to

where $Z̃(k1,k2)=Z̃+z^k1+z^k2$ and $Z̃(k1)=Z̃+z^k1$⁠.

Efficiently deriving matrix elements using the conventional Wick’s theorem requires the introduction of contractions, defined for two creation or annihilation operators ^xb_p and ^xb_q as

where {^xb_p^xb_q} represents a normal-ordered operator string with respect to the reference Fermi vacuum ⟨^xΦ|⋯|^xΦ⟩.³⁶ 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.,

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 $⟨Φx|(aμ)†(aν)†⋯aσaτ|Φw⟩$⁠. 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μ)†$⁠, 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

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]μ$ and $(d[Cn]μ)†$ are expressed purely in terms of the $b̃px$ creation and $b̃p†x$ 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 $⟨Φx|(aμ)†(aν)†⋯aσaτ|Φw⟩$ requires the derivation of contractions between the a^μ and $(aμ)†$ atomic spin-orbital operators rather than the transformed $d[Cn]μ$ and $(d[Cn]μ)†$ 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 $∑Cn1$ 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 $∑Cn̸∋k1$ 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−1n$⁠. 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̃xw$ 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$Φ⟩.

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.³⁶ 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̃xw$ appears as an overall prefactor. The order of summation over the k₁, k₂ indices and $Cn$ can then be swapped to give

The third term in square brackets $∑Cn̸∋k1,k21$ is simply the number of ways to pick n orbitals from the m − 2 zero-overlap orbitals that remain when orbitals k₁ and k₂ are removed, given by $m−2n$⁠. 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

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₁ and m₂ values under the constraint m₁ + m₂ = 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:

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

2. For each term, sum every possible way to distribute m zeros among the contractions such that ∑_im_i = m.

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

4. Multiply the final combined expression with the reduced overlap $S̃xw$⁠.

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³⁶ 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³² and to derive matrix elements between excited configurations with respect to the nonorthogonal reference determinants.

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.

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³² in the form presented in Ref. 5 as

Here, we note that ^xwM^νμ = ^xwW^νμ when there are no zero-overlap orbitals, while $Pνμxw=Pkνμxw$ for one-zero overlap in orbital k. The matrices ^xwM and ^xwP_k correspond, respectively, to the weighted and unweighted co-density matrices discussed in Ref. 5.

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 $⟨Φx|(aμ)†(aν)†aτaσ|Φw⟩$ 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³² to be recovered as

Note that for the m = 1 case with one zero-overlap in orbital k, we have exploited the identities $PkτμxwPkσνxx=PkσμxwPkτνxx$ and $PkτμxwPkσνxw=PkσμxwPkτνxw$ to introduce the simplification ^xwM → ^xwW.

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., $⟨Φij…ab…x|⋯|Φkl…cd…w⟩$⁠. 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²⟩ coupling terms in NOCI expansions.⁶ 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²² or state-specific complete active space SCF.²⁴ 

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.

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