Electronic structure methods that exploit nonorthogonal Slater determinants face the challenge of efficiently computing nonorthogonal matrix elements. In a recent publication [H. G. A. Burton, J. Chem. Phys. **154**, 144109 (2021)], I introduced a generalized extension to the nonorthogonal Wick’s theorem that allows matrix elements to be derived between excited configurations from a pair of reference determinants with a singular nonorthogonal orbital overlap matrix. However, that work only provided explicit expressions for one- and two-body matrix elements between singly- or doubly-excited configurations. Here, this framework is extended to compute generalized nonorthogonal matrix elements between higher-order excitations. Pre-computing and storing intermediate values allows one- and two-body matrix elements to be evaluated with an $O(1)$ scaling relative to the system size, and the LIBGNME computational library is introduced to achieve this in practice. These advances make the evaluation of all nonorthogonal matrix elements almost as easy as their orthogonal counterparts, facilitating a new phase of development in nonorthogonal electronic structure theory.

## I. INTRODUCTION

Nonorthogonality occurs throughout modern quantum chemistry, providing chemically intuitive and compact representations of challenging electronic structures. For example, nonorthogonal configuration interaction (NOCI) uses a linear combination of nonorthogonal Slater determinants, each with bespoke orbitals representing an important electronic configuration, to capture strong static correlation^{1–10} or provide a diabatic representation of electron transfer states.^{11–13} Alternatively, nonorthogonal transition coupling elements arise for orbitally optimized excited states, which can improve the accuracy of challenging core^{14,15} or charge transfer excitations.^{16–23} Beyond providing chemical intuition, nonorthogonality naturally arises in auxiliary-field quantum Monte-Carlo,^{24,25} in coupling terms between states at different molecular geometries,^{26,27} and in spin-projected techniques.^{10,28–30}

Despite presenting many potential benefits, the development of nonorthogonal electronic structure methods is hindered by the difficulty of computing matrix elements between Slater determinants with mutually nonorthogonal orbitals. For orthogonal orbitals, the second-quantized generalized Wick’s theorem^{31} allows arbitrary matrix elements to be easily derived and evaluated using pre-computed one- and two-electron integrals. In contrast, matrix elements for nonorthogonal determinants currently require the first-quantized generalized Slater–Condon rules,^{32} with a computational scaling of at least $O(N3)$ with respect to the number of electrons *N*.^{4}

While the generalized Slater–Condon rules can be used to couple a small number of Slater determinants,^{4–7,33} their additional computational cost quickly makes larger calculations unfeasible. In particular, many recent developments involve multiple orthogonal excited configurations from different nonorthogonal reference determinants, introducing a large number of nonorthogonal matrix elements. This situation arises for the coupling terms between state-specific multi-determinant wave functions^{20–23,34} or post-NOCI techniques for capturing dynamic correlation,^{2,29,30} including perturbative approximations such as NOCI-MP2^{35–37} and NOCI-PT2.^{9} Developing efficient implementations of these methods requires a theory for computing nonorthogonal matrix elements between excited configurations that do not scale with the system size. While equations for specific applications have been independently derived many times (e.g., Refs. 5 and 35), an entirely generalized approach has not yet been developed.

A nonorthogonal extension to Wick’s theorem has been well established for many years^{38,39} and has seen limited use in quantum chemistry.^{29,30,40,41} In the nonorthogonal Wick’s theorem, matrix elements between two Slater determinants can be evaluated using contractions defined from the non-Hermitian transition density matrix.^{38} Combining this framework with Lowdin’s general formula^{42,43} allows nonorthogonal matrix elements to be readily evaluated between excited-configurations from different reference determinants.^{24,25,44} However, the derivation of the nonorthogonal Wick’s theorem is restricted to reference determinants with a non-singular orbital overlap matrix.^{30} This restriction excludes the case where the orbitals are mutually nonorthogonal but the determinants themselves have a zero overlap, which occurs when the orbital overlap matrix is singular but not diagonal.^{45}

In Ref. 46 (henceforth known as Paper I), I introduced an extension to the nonorthogonal Wick’s theorem that allows matrix elements to be computed for nonorthogonal determinants with a zero many-body overlap. This theory allows overlap, one-body, and two-body coupling terms to be derived for excited configurations from arbitrary reference determinants, and provides a second-quantized derivation of the generalized Slater–Condon rules. I then demonstrated that overlap and one-body coupling terms between nonorthogonal excited configurations can be evaluated with $O(1)$ scaling using a set of pre-computed intermediates. However, these explicit derivations were difficult to extend for coupling terms beyond double excitations and the $O(1)$ scaling did not appear to be achievable for two-electron coupling terms.

In the current work, I combine the results of Paper I with the determinantal expansion based on Löwdin’s general formula,^{42,43} giving matrix elements between excitations from arbitrary Slater determinants with any number of singular values in the reference orbital overlap matrix. Furthermore, by introducing pre-computed intermediates, I derive explicit expressions for the overlap, one-body, and two-body coupling terms between arbitrary excitations that scale as $O(1)$ with respect to the number of electrons or basis functions. Finally, I introduce the LIBGNME computational library to provide an open-source implementation of these expressions. These advances significantly reduce the computational cost compared with the previous state-of-the-art, creating a new paradigm for efficient nonorthogonal matrix elements in quantum chemistry.

In Sec. II, I review the nonorthogonal Wick’s theorem and introduce the key results from Paper I, before extending these to arbitrary excitations in Sec. III. Section IV presents numerical data illustrating the computational scaling of these generalized nonorthogonal matrix elements. The conclusions and scope for future applications are summarized in Sec. V.

## II. BACKGROUND THEORY

In this section, I summarize Paper I’s key extensions to the nonorthogonal Wick’s theorem that are required for the derivations presented in this work. Interested readers are directed to Paper I for a comprehensive explanation of these results and to Ref. 31 for details about the generalized Wick’s theorem.

### A. Biorthogonalising the molecular orbitals

This work considers matrix elements between excited configurations constructed from a general pair of mutually nonorthogonal determinants, e.g.,

where $O\u0302$ is an arbitrary operator. Here, $|\Phi x\u3009$ and $|\Phi w\u3009$ are reference Slater determinants built from different molecular orbitals ${\varphi px}$ and ${\varphi pw}$ that are expanded using a common *n*-dimensional atomic orbital basis set {*χ*_{μ}} as

The $C\u22c5p\mu \u22c5x$ denote the orbital coefficients in the molecular orbital basis of determinant $|\Phi x\u3009$ using the nonorthogonal tensor notation of Head-Gordon *et al.*^{47} The corresponding second-quantized operators $bp\u2020x$ and $bpx$ create and annihilate, respectively, an electron in the molecular orbital *p* for the orbital set corresponding to determinant $|\Phi x\u3009$.

To evaluate matrix elements using these nonorthogonal determinants, the molecular orbitals ${\varphi px}$ and ${\varphi pw}$ must be transformed into a biorthogonal set using the Löwdin pairing approach.^{48,49} First, the overlap matrix between the two sets of occupied orbitals is computed as

where *g*_{μν} = ⟨*χ*_{μ}|*χ*_{ν}⟩ is the overlap matrix for the atomic orbital basis set. A singular value decomposition is then performed to diagonalize this overlap matrix, giving the modified occupied orbital coefficients $(C\u0303w)i\u22c5\u22c5\mu $ and $C\u0303i\u22c5\u22c5\mu x$ that satisfy

When one, or more, of the diagonal overlap terms becomes zero (i.e., $S\u0303ixw=0$), the many-body overlap must vanish $\u27e8\Phi x|\Phi w\u27e9=0$. The number of these biorthogonal zero-overlap orbital pairs corresponds to the nullity of the matrix ^{xw}** S** and is denoted by the integer

*m*. The value of

*m*determines which instance of the generalized Slater–Condon rules is required for a matrix element between the reference states $|\Phi x\u3009$ and $|\Phi w\u3009$, as described in Refs. 4 and 32. The product of the remaining non-zero singular values defines the reduced overlap

which also appears in the generalized Slater–Condon rules.^{4}

### B. Extended nonorthogonal Wick’s theorem

Without loss of generality, I will represent operators in the molecular orbital basis for determinant $|\Phi x\u3009$. For example, a one-body operator $f\u0302$ is given by

where

A matrix element such as

can then be evaluated using an extension to Wick’s theorem (see Ref. 31). In particular, each term in Eq. (8) is computed as the sum of all possible ways to fully contract the second-quantized operators, for example,

The phase of each term is given by (−1)^{h} where *h* is the number of intersecting contraction lines.

For the extended nonorthogonal Wick’s theorem outlined in Paper I, a general matrix element in this expansion (e.g., ) is equal to the reduced overlap $S\u0303xw$ multiplied by the product of nonorthogonal contractions,

where $Xqpwx$ and $Yqpxw$ are intermediate values defined for a given pair of determinants [Eq. (11)]. If there are zero-overlap orbital pairs in the biorthogonal basis (*m* > 0), then a sum must be taken over every possible way to assign the *m* zeros to the contractions in each product. This distribution is denoted by indices *m*_{k} assigned to each contraction (i.e., $Xqp(mk)wx$ and $Ypq(mk)xw$), which take values of 0 or 1 and satisfy *∑*_{k}*m*_{k} = *m*. Accounting for these zero-overlap orbital pairs is the central result of Paper I compared with the standard nonorthogonal Wick’s theorem.^{38,45} The individual contractions, represented in the original (un-transformed) orbital basis, are then

where

(*Warning**: for convenience throughout this work, the sign of* $Ypqxw$ *has been reversed relative to Paper I.*)

Notably, the contractions in Eq. (10) are defined with respect to the original orbital coefficients such that the definition of excited configurations is not affected by the biorthogonal transformation. The number of zero-overlap orbital pairs corresponds to the underlying biorthogonal basis and there is no relationship between the *p*, *q* and *m*_{k} indices. Furthermore, the $Xqp(mk)wx$ and $Ypq(mk)xw$ intermediates can be evaluated and stored once for each pair of nonorthogonal reference determinants.

While the biorthogonal orbital coefficients only enter in the definition of the ^{xw}** M** and

^{xw}

**matrices, the number of biorthogonal zero-overlap orbital pairs affects every matrix element through the summation over the allowed combinations of**

*P**m*

_{k}values. This distribution is essential for recovering the different instances of the generalized Slater–Condon rules.

^{46}If

*m*is larger than the total number of contractions, then the corresponding matrix element is strictly zero.

In summary, a nonorthogonal matrix element between excited configurations can be evaluated through the following procedure:

Assemble all fully contracted combinations of the second-quantized operator and excitations, and compute the corresponding phase factors;

For each term, sum every possible way to distribute

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

*m*_{k}} in each term, construct the relevant contribution as a product of fundamental contractions defined in Eqs. (10) and (11);Multiply the combined expression by the reduced overlap $S\u0303xw$.

To demonstrate this process, consider the matrix element $\u27e8\Phi iax|\Phi jbw\u27e9=\u27e8\Phi x|bi\u2020xbaxbb\u2020wbjw|\Phi w\u27e9$. Applying Step 1 gives two contributions,

Each term corresponds to a product of two fundamental contractions with a phase of +1. Taking a sum over the different *m*_{1}, *m*_{2} values satisfying *m*_{1} + *m*_{2} = *m*, and multiplying by the reduced overlap $S\u0303xw$, then gives

which can also be represented as the determinant,

This matrix element reduces to zero if *m* > 2.

## III. EXTENSION FOR ARBITRARY EXCITATIONS

### A. Unification with Löwdin’s general formula

Paper I demonstrated the explicit derivation of nonorthogonal matrix elements for the overlap, one-body, and two-body operators between excited configurations from reference determinants with zero-overlap orbital pairs.^{46} However, those derivations are difficult to generalize for higher-order excitations due to the increasingly large number of fully contracted terms in the expansion of a matrix element. A generalized approach can be achieved using Löwdin’s general formula,^{42,43} which expresses the matrix element for an arbitrary product of creation and annihilation operators as a determinant and automatically includes every fully-contracted term with the correct phase factor. Following previous nonorthogonal derivations,^{25,29,30,44,50} Löwdin’s general formula can be applied to the nonorthogonal contractions defined in Eq. (10) to give

Here, terms in the upper triangle of the determinant correspond to contractions where the annihilation operator is to the left of the creation operator and are multiplied by a factor of −1 to give the correct sign once the determinant is expanded.

We now demonstrate how this representation of the nonorthogonal Wick’s theorem is extended to reference determinants with zero-overlap orbital pairs. First, in contrast to the conventional Wick’s theorem, the entire determinant in Eq. (16) is multiplied by the reduced overlap $S\u0303xw$. Second, the distribution of *m* zero-overlap biorthogonal orbital pairs over products of individual contractions is achieved by exploiting the fact that every term in the expansion of a determinant includes one element from each column of the corresponding matrix. Therefore, a summation can be taken over every possible way to distribute the *m* zero-overlap orbital pairs among the columns in Eq. (16). Inserting the nonorthogonal contractions defined in Eqs. (10) and (11) then gives

where *L* is the total number of creation and annihilation pairs in the operator string. In Eq. (17), the lower triangle of the matrix (including the diagonal) contains “*X*-type” matrix elements while the upper triangle contains “*Y*-type” matrix elements. Reversing the sign of the “*Y*-type” terms defined in Eq. (11b) relative to Paper I ensures that these elements enter the determinant in Eq. (17) with a positive sign contribution, simplifying subsequent derivations. This modified representation of the extended nonorthogonal Wick’s theorem now provides matrix elements between excited configurations from nonorthogonal reference determinants with an arbitrary number of zero-overlap orbital pairs, as demonstrated in Secs. III B–III D.

### B. Overlap elements

The overlap between a general pair of excited configurations is given by

Using the approach presented in Sec. III A, these overlap matrix elements can be expanded as

where *L* is the combined number of excitations in the bra and ket states. Since fundamental contractions *X* and *Y* can be precomputed and stored once for a given pair of reference determinants, the computational cost of evaluating Eq. (19) is only controlled by the total number of excitations *L* and the number of zero-overlap orbitals in the biorthogonal basis *m*. Specifically, the cost of evaluating each determinant in Eq. (19) scales as $O(L3)$, while the number of ways to distribute the *m* zeros over the *L* contractions is $Lm$, giving an overall scaling of $O(L3Lm)$. In comparison, the conventional first-quantized approach evaluates the determinant of the occupied orbital overlap matrix between a pair of excited configurations, giving an $O(N3)$ scaling.^{32} Since the number of excitations is smaller than the total number of electrons by definition, the evaluation of the extended nonorthogonal Wick’s theorem using Eq. (19) is generally much more efficient.

### C. One-body coupling

Representing one-body operators in the molecular orbital basis for determinant $|\Phi x\u3009$ as

allows a general coupling term between two excitations from nonorthogonal reference determinants to be expressed as

Applying Löwdin’s general formula for the extended nonorthogonal Wick’s theorem give the expansion,

where *L*_{x} is the number of excitations in the bra state and *L* is the combined number of excitations in the bra and ket states. Following Ref. 42, a Laplace expansion can then be performed along the row corresponding to the one-body operator index *q* to give

where, for convenience, the element along this row corresponding to operator index *p* is placed first, and the dummy variables *m*_{k} have been rearranged. On the second row of Eq. (23), and those below it, the columns in the minor submatrices have been swapped so that the column corresponding to the one-body operator (index *p*) occupies the position originally held by the corresponding cofactor (e.g., index *i*, *j*, *d*, or *c*), ensuring that these rows all have the same overall sign contribution.

While a naïve implementation of Eq. (23) scales as $O(n2)$ due to the summations over *p* and *q*, this scaling can be removed by introducing a series of intermediate quantities. The first line requires the intermediate,

which takes different values depending on whether the contraction $Xqp(mi)xx$ is assigned to a zero-overlap orbital pair, as determined by *m*_{i}. For the remaining lines in Eq. (23), the rules for multiplying a determinant by a scalar can be exploited to move the summation over *p* and *q* inside the corresponding minor submatrices. Introducing the intermediate quantities

and substituting these into the corresponding column for index *p* gives the overall matrix element expression,

This expression comprises the term $F0(m1)xx$ multiplied by the total overlap of the excited configurations, minus terms where each column of the submatrix is replaced by the corresponding intermediate, e.g., $Fai(m1,m2)xw$.

Computationally, the $F$ intermediates in Eq. (25) can be pre-computed with scaling $O(16n3)$ and stored with scaling $O(16n2)$, where the factor of 16 accounts for the four possible values of (*m*_{i}, *m*_{j}) and the four combinations *xx*, *xw*, *wx*, and *ww*. Once these intermediates have been evaluated, the total cost of computing a one-body matrix element only depends on the total number of excitations in the bra and ket states *L*, which controls the size of the determinants in Eq. (26) and the number of columns that must be replaced with the intermediates, e.g., $Fai(m1,m2)xw$. The cost of evaluating each determinant in Eq. (26) scales as $O(L3)$ and there are *L* + 1 determinants that must be computed. The summation over the *m*_{i} values depends on the total number of way to distribute *m* zero-overlap orbitals between the *L* + 1 contractions, given by $L+1m$. Consequently, the overall scaling for each one-body matrix element is $O(L3(L+1)L+1m)$, which is constant with respect to the number of electrons or the size of the basis set. In comparison, applying the generalized Slater–Condon rules for each pair of determinants requires the biorthogonalization of the *full* set of occupied orbitals, with an iterative $O(N3)$ scaling, followed by an $O(n2)$ contraction of the co-density matrices with the one-electron integrals.^{4,32} Evidently, the new approach described here is significantly more efficient when a large number of coupling elements are required.

### D. Two-body coupling

Finally, two-body operators can be represented as

where the two-electron integrals $vprqsx=(pr|qs)x$ are expressed in the MO basis for reference $|\Phi x\u3009$ and are represented using Mulliken notation.^{51} A general coupling term between two excitations from nonorthogonal reference determinants is then given by

Deriving expressions for these matrix elements using the Laplace expansion approach is much more involved than the one-body operators and is much harder for the reader to follow. Instead, an alternative approach can be used that constructs effective zero- or one-body operators by partially contracting subunits containing two or four of the indices *p*, *r*, *q*, *s*, and the one-body framework established in Sec. III C can then be applied. The final expression is given by

where each constituent equation is derived in detail below.

First, contractions containing the subunits and lead to the “zeroth-order” contribution,

where the intermediate terms have been defined,

The contribution from Eq. (30) is analogous to the one-body $F0(m1)xx$ contribution in the first line in Eq. (26).

Next, the partially contracted subunits and can be used to define a new effective one-body operator,

where the intermediate term in Eq. (31a) has been employed and the factor of two arises from the equivalence of terms such as and under the permutation of the dummy indices *p*, *q*, *r*, *s*. The contribution from these partially contracted subunits can then be evaluated using the rules established for one-body operators in Sec. III C by introducing the additional intermediate terms,

Each of these intermediates can be computed with a computational scaling of $O(8n3)$, and the overall storage cost for a pair of reference determinants is 32*n*^{2}. The one-body terms that are analogous to the first line in Eq. (26), which represent the contractions and have already been included in Eq. (30). Therefore, the unique contributions of these effective one-body operators to the full two-body matrix element in Eq. (28) is

The summation within the parentheses includes the replacement of each column of the overlap determinant with the corresponding intermediate terms from Eq. (33).

The remaining terms include cases where all the indices *p*, *q*, *r*, *s* are contracted with operators that occur in the bra or ket excitation strings. An effective one-body operator can be constructed by considering partial contractions with the form, e.g.,

where the indices *p*, *r* and *p*, *s* are contracted with bra or ket excitation operators, giving effective operators with the form, e.g.,

where the permutation of the dummy indices *p*, *q*, *r*, *s* has been exploited to incorporate the “exchange-like” terms. The operators in Eq. (36) contain a phase factor *ϕ*_{ab} that takes the value +1 if *a* and *b* correspond to the same excitation or are separated by an odd number of excitations, and −1 if *a* and *b* correspond to different excitations separated by an even number of excitations, as illustrated for the example $\u27e8\Phi ijabx|v\u0302|\Phi klcdw\u27e9$ by the checkerboard pattern,

Every effective one-body operator defined for the *L*^{2} combinations of a creation and annihilation operator pair (excluding *p*, *q*, *r*, *s*) must be considered. The contribution to the full two-body matrix element can then be evaluated using the one-body approach described in Sec. III C by introducing a final set of intermediates

The standard two-electron integral symmetry relation $Jab,cd(mi,mj,mk,ml)xw,yz=Jcd,ab(mk,ml,mi,mj)yz,xw$ allows intermediates for the remaining combinations of indices to be obtained. Each intermediate can be computed with a computational scaling of $O(16n5)$, where the factor of 16 comes from the possible combinations of (*m*_{i}, *m*_{j}, *m*_{k}, *m*_{l}). Taking into account the two-electron symmetry relation, the maximum storage requirement is 160*n*^{4}. In practice, the storage may be reduced if not all combinations of (*m*_{i}, *m*_{j}, *m*_{k}, *m*_{l}) are required (i.e., *m* < 4) or if excitations are only considered within an active orbital space.

When applying the one-body framework (Sec. III C) to the effective operators defined in Eq. (36), the terms that arise from contracting the *q* and *s* indices are discarded as these are already taken into account by Eq. (34). These terms correspond to the first line in Eq. (26). Therefore, the remaining unique contributions to the two-body matrix element are

Individual lines in this expression correspond to the effective one-body operators constructed from a different pair of the creation and annihilation operators selected from excitation strings [i.e., (*ai*), (*aj*), …, (*kc*)]. The sum within a line corresponds to the replacement of each column by the corresponding intermediate for the remaining excitation indices, where the excitation indices used to build the effective one-body operators have been removed. Every contribution also includes a summation over the possible ways to distribute the *m* zero-overlap orbital pairs over the contractions. The factor of 1/2 accounts for the double counting of contributions, for example, the term $Jai,bj(m1,m2,m3,m4)xx,xx$ appears for the effective one-body operator with indices (*ai*) and (*bj*) due to the symmetry $Jai,bj(m1,m2,m3,m4)xx,xx=Jbj,ai(m1,m2,m3,m4)xx,xx$.

Once all the required intermediates have been computed and stored, the computational cost of evaluating a two-body matrix element using this approach is determined by cost of computing Eqs. (30), (34), and (38). For each term, there are $L+2m$ ways to distribute *m* zero-overlap orbital pairs among the *L* + 2 contractions. Equation (30) involves the computation of only one *L* × *L* determinant, giving an $O(L3L+2m)$ scaling. Equation (34) requires the computation of *L* determinants where each column in the overlap determinant is replaced by intermediates of the form Eq. (33), giving a scaling of $O(L4L+2m)$. Finally, Eq. (38) involves the computation of *L* − 1 determinants with dimensions (*L* − 1) × (*L* − 1) for each pair of creation and annihilation operators in the excitation strings, giving a scaling of $O(L2(L\u22121)4L+2m)$. The asymptotic scaling for individual two-body matrix elements using these intermediates is, therefore, $O(L6L+2m)$. Although this computational cost increases rapidly with the number of excitations in a two-body matrix element, it has $O(1)$ scaling with respect to the number of basis functions or electrons. In contrast, applying the generalized Slater–Condon rules for two-body operators scales as $O(n4)$ for each matrix element and rapidly makes the computation of a large number of matrix elements unfeasible.

## IV. ILLUSTRATION OF COMPUTATIONAL SCALING

The extended nonorthogonal Wick’s theorem for arbitrary excitations has been implemented in a developmental open-source C++ package LIBGNME, available for download at Ref. 52. The primary advantage of the extended nonorthogonal Wick’s theorem over the generalized Slater–Condon rules is the $O(1)$ scaling with respect to the number of basis functions or electrons that can be achieved once all the required intermediates have been pre-computed. This scaling means that the computation of nonorthogonal matrix elements for excited configurations becomes almost as straightforward as the conventional orthogonal Slater–Condon rules or Wick’s theorem.

The acceleration relative to the generalized Slater–Condon rules can be demonstrated by comparing the average time required to compute one- and two-body matrix elements between singly excited or doubly excited configurations with increasingly large correlation-consistent basis sets.^{53} Illustrative calculations using LIBGNME have been performed for two nonorthogonal reference determinants corresponding to the spin–flip pair of spin-symmetry-broken unrestricted Hartree–Fock solutions for the ground state of H_{2}O at a bond angle of 104.5 and *R*(O–H) = 1.35 Å.

Average timings for each one-body matrix element are compared between the extended nonorthogonal Wick’s theorem and the generalized Slater–Condon rules in Fig. 1. This average is computed using all the high-spin single or double excited coupling terms within a (10, 13) active space constructed from the lowest-energy molecular orbitals. This choice corresponds to 40 single and 280 double excitations, giving a total of 1600 single–single and 78 400 double–double coupling terms for each basis set. The statistical distribution is assessed using 48 replicas of each calculation. As expected, these data demonstrate the $O(1)$ scaling of the extended nonorthogonal Wick’s theorem with respect to the basis set size, while the generalized Slater–Condon rules scale asymptotically as $O(n2)$. For small basis sets, the scaling of the generalized Slater–Condon rules becomes constant as the computational cost is dominated by biorthogonalizing the occupied orbitals in each pair of excited configurations. In the largest basis set considered (cc-pV5Z), the extended nonorthogonal Wick’s theorem provides a computational acceleration of nearly 3 orders of magnitude relative to the generalized Slater–Condon rules.

Analogous average timings for two-body matrix elements between singly excited and doubly excited configurations are presented in Fig. 2. Here, the generalized Slater–Condon rules show a near-ideal $O(n4)$ computational scaling, while the extended nonorthogonal Wick’s theorem has a constant cost for all basis sets, as predicted. The larger *n*^{4} scaling of two-body compared with one-body matrix elements means that the extended nonorthogonal Wick’s theorem offers an even greater advantage over the generalized Slater–Condon rules, as demonstrated by the six orders of magnitude acceleration achieved for the $\u27e8\Phi iax|v\u0302|\Phi jbw\u27e9$ matrix elements in the cc-pVQZ basis set. These numerical results highlight the significant computational advance provided by the generalized nonorthogonal matrix elements compared with the previous state-of-the-art.

## V. DISCUSSION AND OUTLOOK

This work has extended the nonorthogonal Wick’s theorem^{38,45} and the results of Ref. 46 to derive coupling elements between arbitrary excited configurations from a pair of reference Slater determinants with a singular nonorthogonal orbital overlap matrix. By pre-computing and storing various intermediate terms, subsequent one- and two-body matrix elements can be evaluated with a computational cost that scales $O(1)$ with respect to the number of basis functions or electrons. These developments provide a significant improvement over the commonly used generalized Slater–Condon rules, which asymptotically scale as $O(n2)$ or $O(n4)$ for one- or two-body matrix elements, respectively.

While the nonorthogonal Wick’s theorem is well established^{38,45} and explicit expressions for certain nonorthogonal matrix elements have previously been reported, this work presents an entirely general and unified framework for developing nonorthogonal techniques. The current approach shares many similarities with the nonorthogonal matrix elements derived by Mahajan and Sharma,^{24,25} including the use of the determinantal expansion of Wick’s theorem. However, in contrast to their work, the derivations presented here are entirely generalized to cases where singular values occur in the biorthogonalization of the reference orbitals, providing a unification with the various instances of the generalized Slater–Condon rules.^{32} Consequently, this framework establishes the most general and flexible version of Wick’s theorem for computing matrix elements between arbitrary Slater determinants.

In practice, the bottleneck for these generalized nonorthogonal matrix elements is the computation and storage of the intermediate terms required for two-body operators. While the evaluation of the intermediates in Eq. (37) has the same $O(n5)$ scaling as a standard two-electron integral transformation, the cost of storing all these intermediates for a given pair of determinants scales as 160*n*^{4} when the two-electron symmetry is taken into account. The storage of two-electron integrals is already challenging for larger basis sets, so increasing this by a factor of 160 represents a significant computational overhead. This cost can be reduced if there are no zero-overlap orbital pairs in the reference determinants (i.e., *m* = 0), such that only a subset of (*m*_{i}, *m*_{j}, *m*_{k}, *m*_{l}) combinations are required, or if excitations are only considered within an active orbital space. Alternatively, we can expect techniques such as the Cholesky decomposition^{54} or resolution-of-the-identity^{55} to provide further reductions in the storage and computational cost of computing the two-electron intermediates. These areas for improvement are beyond the scope of the current work but will form the focus for future computational research and development.

Until now, the development of advanced nonorthogonal techniques in electronic structure theory has been hindered by the computational cost of computing arbitrary nonorthogonal matrix elements. The generalized extension to the nonorthogonal Wick’s theorem,^{38,45} introduced in Ref. 46 and further developed in this work, now offers a universal route to overcome this challenge. Moving forwards, the ability to rapidly compute nonorthogonal matrix elements between excited configurations will allow on-the-fly implementations of NOCI extensions such as NOCI-PT2,^{9} NOCI-MP2,^{36,37} and nonorthogonal configuration interaction singles (NOCIS).^{14,15} Furthermore, the generalization to arbitrary excitations will enable the evaluation of coupling terms between excited state-specific multi-configurational wave functions^{20–23,34} that are prohibitively expensive using the generalized Slater–Condon rules. Ultimately, the generality of this framework, and the availability of the open-source LIBGNME implementation,^{52} will catalyze a new phase of development in nonorthogonal electronic structure theory.

## ACKNOWLEDGMENTS

H.G.A.B. was supported by New College, Oxford through the Astor Junior Research Fellowship. The author thanks Nicholas Lee and Alex Thom for helpful feedback during the preparation of this manuscript, and Rebecca Lloyd for careful proofreading.

## AUTHOR DECLARATIONS

### Conflict of Interest

The author has no conflicts to disclose.

### Author Contributions

**Hugh G. A. Burton**: Conceptualization (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Project administration (lead); Software (lead); Writing – original draft (lead); Writing – review & editing (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.