Generalized nonorthogonal matrix elements for arbitrary excitations

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 nonorthogonal extension to Wick's theorem that allows matrix elements to be derived between excited configurations from different reference determinants. 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 arbitrary excitations. Pre-computing and storing intermediate values allows one- and two-body matrix elements to be evaluated with an $\mathcal{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 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 structure.][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(N 3 ) 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][5][6][7]33 their additional computational cost quickly makes larger calculations unfeasible. In articular, 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 a) Electronic mail: hgaburton@gmail.comfunctions [20][21][22][23]34 or post-NOCI techniques for capturing dynamic correlation, 2,29,30 including perturbative approximations such as NOCI-MP2 [35][36][37] and NOCI-PT2.9 Developing efficient implementations of these methods requires a theory for computing nonorthogonal matrix elements between excited configurations that does 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.In this work, I derive a theory that achieves this aim for arbitrary one-and two-body nonorthogonal matrix elements.
The first significant step towards an entirely general theory for nonorthogonal matrix elements was the nonorthogonal extension to Wick's theorem. 38However, that approach is restricted to reference determinants with a non-zero manybody overlap, excluding the case where the orbitals are mutually nonorthogonal but the determinants themselves have a zero overlap, 39 and has seen limited use in quantum chemistry. 29,30,40Similarly, Mahajan and Sharma introduced an approach for efficiently evaluating Hamiltonian matrix elements between an arbitrary multi-determinant wave function and a single reference determinant built from different orbitals. 24,25While promising, this theory also fails for orthogonal reference determinants with mutually nonorthogonal orbitals.
In Ref. 41 (henceforth known as Paper 1), I introduced a generalized nonorthogonal extension to Wick's theorem that allows matrix elements to be computed for any pair of nonorthogonal determinants, even if their many-body overlap is zero.This theory allows arbitrary overlap, one-body, and twobody coupling terms to be derived for excited configurations from different reference determinants, and enables a secondquantized 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 extend the results of Paper 1 by introducing a simplified framework for computing nonorthogo-arXiv:2208.10208v1[physics.chem-ph]22 Aug 2022 nal matrix elements between arbitrary excited configurations from different reference determinants.This framework combines the extended nonorthogonal Wick's theorem derived in Paper 1 with a determinantal expansion of Wick's theorem based on Löwdin's general formula, 42,43 giving universal expressions for arbitrary excitations.Furthermore, by introducing pre-computed intermediates, I derive explicit expressions for the overlap, one-body, and two-body coupling terms between arbitrary excitations which scale as O(1) with respect to the number of electrons or basis functions.These advances significantly reduce the computational cost compared to the previous state-of-the-art, creating a new paradigm for efficient nonorthogonal matrix elements in quantum chemistry.
In Section II, I introduce the key results from Paper 1, before extending these to arbitrary excitations in Section III.Section IV presents numerical data illustrating the computational scaling of these generalized nonorthogonal matrix elements.The primary conclusions and scope for future applications are summarized in Section V.

II. BACKGROUND THEORY
In this Section, I summarize the key aspects of the extended nonorthogonal Wick's theorem derived in Paper 1 that are required for the derivations presented in this work.Interested readers are directed to Paper 1 for a comprehensive explanation of these results and to Ref. 31 for a detailed description of 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 Ô is an arbitrary operator.Here, | x Φ and | w Φ are reference Slater determinants built from different molecular orbitals { x φ p } and { w φ p } that are expanded using a common ndimensional atomic orbital basis set {χ µ } as The x C µ• •p denote the orbital coefficients in the molecular orbital basis of determinant | x Φ using the nonorthogonal tensor notation of Head-Gordon et al. 44 The corresponding secondquantized operators x b † p and x b p create and annihilate, respectively, an electron in the molecular orbital p for the orbital set corresponding to determinant | x Φ .
To evaluate matrix elements using these nonorthogonal determinants, the molecular orbitals { x φ p } and { w φ p } must be transformed into a biorthogonal set using the Löwdin pairing approach. 45,46First, 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 ( w C) When one, or more, of the diagonal overlap terms becomes zero (i.e., xw S i = 0), the many-body overlap must vanish x Φ| w Φ = 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 | x Φ and | w Φ , as described in Refs.32 and 4. 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 | x Φ .For example, a one-body operator f 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 secondquantized 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 1, a general matrix element in this expansion (e.g., ) is equal to the reduced overlap xw S multiplied by the product of nonorthogonal contractions where wx X qp and xw Y qp 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., wx X (m k ) qp and xw Y (m k ) pq ), which take values of 0 or 1 and satisfy k m k = m.The individual contractions, represented in the original (untransformed) orbital basis, are then defined as where (Warning: for convenience throughout this work, the sign of xw Y pq has been reversed relative to Paper 1.) 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 wx X (m k ) qp and xw Y (m k ) pq 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 P matrices, the number of biorthogonal zero-overlap orbital pairs affects every matrix element through the summation over the allowed combinations of m k values.This distribution is essential for recovering the different instances of the generalized Slater-Condon rules. 41f m is larger than the total number of contractions, then the corresponding matrix element is strictly zero.
In summary, an arbitrary nonorthogonal matrix element between excited configurations can be evaluated through the following procedure: 1. Assemble all fully contracted combinations of the second-quantized operator and excitations, and compute the corresponding phase factors; 2. For each term, sum every possible way to distribute m zeros among the contractions {m k } such that k m k = m; 3. 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); 4. Multiply the combined expression by the reduced overlap xw S .
To demonstrate this process, consider the matrix element 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 xw S , then gives which can also be represented as the determinant This matrix element reduces to zero if m > 2.

A. Unification with Löwdin's general formula
Paper 1 demonstrated the explicit derivation of matrix elements for the overlap, one-body, and two-body operators between excited configurations up to double excitations using the extended nonorthogonal Wick's theorem. 41However, 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 fully 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.Since Löwdin's general formula is only based on the different diagrammatic combinations of contractions and their corresponding phases, it can be readily extended 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.
In contrast to the conventional form of Löwdin's general formula, the entire determinant is multiplied by the reduced overlap xw S .Furthermore, the distribution of m zero-overlap biorthogonal orbital pairs over products of individual contractions can be 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 1 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 using Löwdin's general formula now allows matrix elements between arbitrary nonorthogonal excited configurations to be evaluated, as demonstrated in Sections III B-III D.

B. Overlap Elements
The overlap between a general pair of excited configurations is given by Using Lowdin's general representation for the extended nonorthogonal Wick's theorem presented in Section 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 the 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(L 3 ), while the number of ways to distribute the m zeros over the L contractions is L m , giving an overall scaling of O(L 3 L m ).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(N 3 ) scaling. 32Since 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 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(n 2 ) 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 xx X (m i ) qp 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 wx F ww F and substituting these into the corresponding column for index p gives the overall matrix element expression This expression comprises the term xx F (m 1 ) 0 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., xw F (m 1 ,m 2 ) ai .Computationally, the F intermediates in Eq. ( 25) can be pre-computed with scaling O(16n 3 ) and stored with scaling O(16n 2 ), 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.xw F (m 1 ,m 2 ) ai .The cost of evaluating each determinant in Eq. ( 26) scales as O(L 3 ) 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+1 m .Consequently, the overall scaling for each one-body matrix element is O(L 3 (L + 1) L+1 m ), 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 biorthogonalisation of the full set of occupied orbitals, with an iterative O(N 3 ) scaling, followed by an O(n 2 ) contraction of the co-density matrices with the one-electron integrals. 4,32Evidently, 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 x v prqs = x (pr|qs) are expressed in the MO basis for reference | x Φ and are represented using Mulliken notation. 47A 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 Section III C can then be applied.The final expression is given by 30) + Eq. ( 34) + Eq. ( 38), where each constituent equation is derived in detail below.
First, contractions containing the subunits where the intermediate terms have been defined x ṽ(m i ) pr = qs The contribution from Eq. ( 30) is analogous to the one-body xx F (m 1 ) 0 contribution in the first line in Eq. ( 26).
Next, the partially-contracted subunits where the intermediate term in Eq. (31a) has been employed and the factor of two arises from the equivalence of terms such x b s 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 Section III C by introducing the additional intermediate terms Each of these intermediates can be computed with a computational scaling of O(8n 3 ) and the overall storage cost for a pair of reference determinants is 32n 2 .The one-body terms that are analogous to the first line in Eq. ( 26), which represent the contractions x b s , 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 x Φ ab i j |v| w Φ cd kl by the chequerboard 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 Section III C by introducing a final set of intermediates The standard two-electron integral symmetry relation xw,yz J allows intermediates for the remaining combinations of indices to be obtained.Each intermediate can be computed with a computational scaling of O(16n 5 ), where the factor of 16 comes from the possible combinations of (m i , m j , m k , m l ) and the maximum storage requirement is 256n 4 .In practice, this storage requirement is 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 (Section 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 ai,l j 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), (a j), . . ., (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 xx,xx J (m 1 ,m 2 ,m 3 ,m 4 ) ai,b j appears for the effective one-body operator with indices (ai) and (b j) due to the symmetry xx,xx J (m 1 ,m 2 ,m 3 ,m 4 ) ai,b j . 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 the the cost of computing Eqs.(30), (34) and (38).For each term, there are L+2 m ways to distribute m zero-overlap orbital pairs among the L + 2 contractions.Eq. ( 30) involves the computation of only one L × L determinant, giving an O(L 3 L+2 m ) scaling.Eq. ( 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(L 4 L+2 m ).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(L 2 (L − 1) 4 L+2 m ).The asymptotic scaling for individual two-body matrix elements using these intermediates is therefore O(L 6 L+2 m ).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(n 4 ) 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 opensource C++ package LibGNME, available for download at Ref. 48.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-or doubly-excited configurations with increasingly large correlation-consistent basis sets. 49Illustrative 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 singly-or doubly-excited coupling terms within a (13,10) active space comprising the lowestenergy molecular orbitals.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(n 2 ).For small basis sets, the scaling of the generalized Slater-Condon rules becomes constant as the computational cost is dominated by biorthogonalising 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-and doubly-excited configurations are presented in Fig. 2. Here, the generalized Slater-Condon rules show a near-ideal O(n 4 ) 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 to 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 x Φ a i | b| w Φ b j matrix elements in the cc-pVQZ basis set.These numerical results highlight the significant computational advance provided by the generalized nonorthogonal matrix elements compared to the previous state-of-the-art.

V. DISCUSSION AND OUTLOOK
This work has extended the generalized nonorthogonal matrix elements introduced in Ref. 41 to derive coupling elements for arbitrary excited configurations.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(n 2 ) or O(n 4 ) for one-or two-body matrix elements, respectively.
While 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 gener-alized to cases where zero-overlap orbital pairs occur in the biorthogonalisation of the reference determinants, providing a unification with the various instances of the generalized Slater-Condon rules. 32Consequently, 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(n 5 ) scaling as a standard two-electron integral transformation, the cost of storing all these intermediates for a given pair of determinants scales as 256 n 4 .The storage of two-electron integrals is already challenging for larger basis sets, so increasing this by a factor of 256 represents a significant computational overhead.Fortunately, this cost can be reduced if there are no zerooverlap 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, techniques such as the Cholesky decomposition 50 or resolution-of-the-identity 51 may provide further computational savings.
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 framework introduced in Ref. 41, and extended in this work, now offers a practical 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 NOCIS. 14,15Furthermore, the generalisation to arbitrary excitations will enable the evaluation of coupling terms between excited state-specific multi-configurational wave functions [20][21][22][23]34 that are prohibitively expensive using the generalized Slater-Condon rules. Ultmately, the generality of this framework, and the availability of the open-source LibGNME implementation, 48 will cataylse a new phase of development in nonorthogonal electronic structure theory.

- 1 FIG. 1 : 4 FIG. 2 :
FIG. 1: Comparison of the average timing for computingone-body matrix elements in H 2 O using the extended nonorthogonal Wick's theorem (blue) and the generalized Slater-Condon rules (red).Data are plotted on a log-log scale.