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.

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 correlation1–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 core14,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 theorem31 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 functions20–23,34 or post-NOCI techniques for capturing dynamic correlation,2,29,30 including perturbative approximations such as NOCI-MP235–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 years38,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 formula42,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.

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.

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

(1)

where Ô is an arbitrary operator. Here, |Φx and |Φw are reference Slater determinants built from different molecular orbitals {ϕpx} and {ϕpw} that are expanded using a common n-dimensional atomic orbital basis set {χμ} as

(2)

The Cpμx denote the orbital coefficients in the molecular orbital basis of determinant |Φx using the nonorthogonal tensor notation of Head-Gordon et al.47 The corresponding second-quantized operators bpx and bpx 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 {ϕpx} and {ϕ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

(3)

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̃w)iμ and C̃iμx that satisfy

(4)

When one, or more, of the diagonal overlap terms becomes zero (i.e., S̃ixw=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 xwS 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. 4 and 32. The product of the remaining non-zero singular values defines the reduced overlap

(5)

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

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

(6)

where

(7)

A matrix element such as

(8)

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,

(9)

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̃xw multiplied by the product of nonorthogonal contractions,

(10a)
(10b)

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 mk assigned to each contraction (i.e., Xqp(mk)wx and Ypq(mk)xw), which take values of 0 or 1 and satisfy kmk = 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

(11a)
(11b)

where

(12a)
(12b)
(12c)
(12d)

(Warning: for convenience throughout this work, the sign ofYpqxwhas 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 mk 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 xwM and xwP matrices, the number of biorthogonal zero-overlap orbital pairs affects every matrix element through the summation over the allowed combinations of mk 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:

  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 {mk} such that kmk = m;

  3. For every set of {mk} 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 S̃xw.

To demonstrate this process, consider the matrix element Φiax|Φjbw=Φx|bixbaxbbwbjw|Φw. Applying Step 1 gives two contributions,

(13)

Each term corresponds to a product of two fundamental contractions with a phase of +1. Taking a sum over the different m1, m2 values satisfying m1 + m2 = m, and multiplying by the reduced overlap S̃xw, then gives

(14)

which can also be represented as the determinant,

(15)

This matrix element reduces to zero if m > 2.

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

(16)

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̃xw. 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

(17)

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.

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

(18)

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

(19)

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.

Representing one-body operators in the molecular orbital basis for determinant |Φx as

(20)

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

(21)

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

(22)

where Lx 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

(23)

where, for convenience, the element along this row corresponding to operator index p is placed first, and the dummy variables mk 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,

(24)

which takes different values depending on whether the contraction Xqp(mi)xx is assigned to a zero-overlap orbital pair, as determined by mi. 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

(25a)
(25b)
(25c)
(25d)

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

(26)

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 (mi, mj) 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 mi 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.

Finally, two-body operators can be represented as

(27)

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

(28)

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

(29)

where each constituent equation is derived in detail below.

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

(30)

where the intermediate terms have been defined,

(31a)
(31b)

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,

(32)

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,

(33a)
(33b)
(33c)
(33d)

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 32n2. 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

(34)

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

(35)

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

(36a)
(36b)
(36c)
(36d)

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 Φijabx|v̂|Φklcdw by the checkerboard pattern,

Every effective one-body operator defined for the L2 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

(37a)
(37b)
(37c)
(37d)
(37e)
(37f)
(37g)
(37h)
(37i)
(37j)

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 (mi, mj, mk, ml). Taking into account the two-electron symmetry relation, the maximum storage requirement is 160n4. In practice, the storage may be reduced if not all combinations of (mi, mj, mk, ml) 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

(38)

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(L1)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.

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

FIG. 1.

Comparison of the average timing for computing one-body matrix elements in H2O using the extended nonorthogonal Wick’s theorem (blue) and the generalized Slater–Condon rules (red). A total of 1600 single–single and 78 400 double–double coupling terms are computed for each basis set and only high-spin excitations are considered. Data are plotted on a log–log scale.

FIG. 1.

Comparison of the average timing for computing one-body matrix elements in H2O using the extended nonorthogonal Wick’s theorem (blue) and the generalized Slater–Condon rules (red). A total of 1600 single–single and 78 400 double–double coupling terms are computed for each basis set and only high-spin excitations are considered. Data are plotted on a log–log scale.

Close modal

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 n4 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 Φiax|v̂|Φjbw 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.

FIG. 2.

Comparison of the average timing for computing two-body matrix elements in H2O using the extended nonorthogonal Wick’s theorem (blue) and the generalized Slater–Condon rules (red). A total of 1600 single–single and 78 400 double–double coupling terms are computed for each basis set and only high-spin excitations are considered. Data are plotted on a log–log scale.

FIG. 2.

Comparison of the average timing for computing two-body matrix elements in H2O using the extended nonorthogonal Wick’s theorem (blue) and the generalized Slater–Condon rules (red). A total of 1600 single–single and 78 400 double–double coupling terms are computed for each basis set and only high-spin excitations are considered. Data are plotted on a log–log scale.

Close modal

This work has extended the nonorthogonal Wick’s theorem38,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 established38,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 160n4 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 (mi, mj, mk, ml) combinations are required, or if excitations are only considered within an active orbital space. Alternatively, we can expect techniques such as the Cholesky decomposition54 or resolution-of-the-identity55 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 functions20–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.

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.

The author has no conflicts to disclose.

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

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

1.
H.
Koch
and
E.
Dalgaard
,
Chem. Phys. Lett.
212
,
193
(
1993
).
2.
S.
Ten-no
,
Theor. Chim. Acta
98
,
182
(
1998
).
3.
P. Y.
Ayala
and
H. B.
Schlegel
,
J. Chem. Phys.
108
,
7560
(
1998
).
4.
A. J. W.
Thom
and
M.
Head-Gordon
,
J. Chem. Phys.
131
,
124113
(
2009
).
5.
E. J.
Sundstrom
and
M.
Head-Gordon
,
J. Chem. Phys.
140
,
114103
(
2014
).
6.
N. J.
Mayhall
,
P. R.
Horn
,
E. J.
Sundstrom
, and
M.
Head-Gordon
,
Phys. Chem. Chem. Phys.
16
,
22694
(
2014
).
7.
H. G. A.
Burton
and
A. J. W.
Thom
,
J. Chem. Theory Comput.
15
,
4851
(
2019
).
8.
B. C.
Huynh
and
A. J. W.
Thom
,
J. Chem. Theory Comput.
16
,
904
(
2020
).
9.
H. G. A.
Burton
and
A. J. W.
Thom
,
J. Chem. Theory Comput.
16
,
5586
(
2020
).
10.
N.
Lee
and
A. J. W.
Thom
,
J. Chem. Theory Comput.
18
,
710
(
2022
).
11.
K. T.
Jensen
,
R. L.
Benson
,
S.
Cardamone
, and
A. J. W.
Thom
,
J. Chem. Theory Comput.
14
,
4629
(
2018
).
12.
M.
Wibowo
,
R.
Broer
, and
R. W. A.
Havenith
,
Comput. Theor. Chem.
1116
,
190
(
2017
).
13.
T. P.
Straatsma
,
R.
Broer
,
A.
Sánchez-Mansilla
,
C.
Sousa
, and
C.
de Graaf
,
J. Chem. Theory Comput.
18
,
3549
(
2022
).
14.
K. J.
Oosterbaan
,
A. F.
White
, and
M.
Head-Gordon
,
J. Chem. Phys.
149
,
044116
(
2018
).
15.
K. J.
Oosterbaan
,
A. F.
White
, and
M.
Head-Gordon
,
J. Chem. Theory Comput.
15
,
2966
(
2019
).
16.
A. T. B.
Gilbert
,
N. A.
Besley
, and
P. M. W.
Gill
,
J. Phys. Chem. A
112
,
13164
(
2008
).
17.
D.
Hait
and
M.
Head-Gordon
,
J. Chem. Theory Comput.
16
,
1699
(
2020
).
18.
K.
Carter-Fenk
and
J. M.
Herbert
,
J. Chem. Theory Comput.
16
,
5067
(
2020
).
19.
G.
Levi
,
A. V.
Ivanov
, and
H.
Jónsson
,
J. Chem. Theory Comput.
16
,
6968
(
2020
).
20.
J. A. R.
Shea
and
E.
Neuscamman
,
J. Chem. Phys.
149
,
081101
(
2018
).
21.
T. S.
Hardikar
and
E.
Neuscamman
,
J. Chem. Phys.
153
,
164108
(
2020
).
22.
L. N.
Tran
,
J. A. R.
Shea
, and
E.
Neuscamman
,
J. Chem. Theory Comput.
15
,
4790
(
2019
).
23.
L. N.
Tran
and
E.
Neuscamman
,
J. Phys. Chem. A
124
,
8273
(
2020
).
24.
A.
Mahajan
and
S.
Sharma
,
J. Chem. Phys.
153
,
194108
(
2020
).
25.
A.
Mahajan
and
S.
Sharma
,
J. Chem. Theory Comput.
17
,
4786
(
2021
).
26.
H.-T.
Chen
,
J.
Chen
,
D. V.
Cofer-Shabica
,
Z.
Zhou
,
V.
Athavale
,
G.
Medders
,
M. F. S. J.
Menger
,
J. E.
Subotnik
, and
Z.
Jin
,
J. Chem. Theory Comput.
18
,
3296
(
2022
).
27.
S.
Lee
,
Y.
Horbatenko
,
M.
Filatov
, and
C. H.
Choi
,
J. Phys. Chem. Lett.
12
,
4722
(
2021
).
28.
G. E.
Scuseria
,
C. A.
Jiménez-Hoyos
,
T. M.
Henderson
,
K.
Samanta
, and
J. K.
Ellis
,
J. Chem. Phys.
135
,
124108
(
2011
).
29.
T.
Tsuchimochi
and
S.
Ten-no
,
J. Chem. Phys.
144
,
011101
(
2016
).
30.
T.
Tsuchimochi
and
S.
Ten-no
,
J. Chem. Theory Comput.
12
,
1741
(
2016
).
31.
I.
Shavitt
and
R.
Bartlett
,
Many-Body Methods in Chemistry and Physics
(
Cambridge University Press
,
2009
).
32.
I.
Mayer
,
Simple Theorems, Proofs, and Derivations in Quantum Chemistry
(
Springer
,
2003
).
33.
H. G. A.
Burton
and
A. J. W.
Thom
,
J. Chem. Theory Comput.
12
,
167
(
2016
).
34.
F.
Kossoski
,
Y.
Damour
, and
P.-F.
Loos
,
J. Phys. Chem. Lett.
13
,
4342
(
2022
).
35.
S. R.
Yost
,
T.
Kowalczyk
, and
T.
Van Voorhis
,
J. Chem. Phys.
139
,
174104
(
2013
).
36.
S. R.
Yost
and
M.
Head-Gordon
,
J. Chem. Phys.
145
,
054105
(
2016
).
37.
S. R.
Yost
and
M.
Head-Gordon
,
J. Chem. Theory Comput.
14
,
4791
(
2019
).
38.
J.-P.
Blaizot
and
G.
Ripka
,
Quantum Theory of Finite Systems
(
MIT Press
,
1986
).
39.
J.
Hendeković
,
M.
Pavlović
, and
F.
Sokolić
,
Chem. Phys. Lett.
77
,
382
(
1981
).
40.
C. A.
Jiménez-Hoyos
,
R.
Rodríguez-Guzmán
, and
G. E.
Scuseria
,
Phys. Rev. A
86
,
052102
(
2012
).
41.
J.
Nite
and
C. A.
Jiménez-Hoyos
, “
Efficient multi-configurational wavefunction method with dynamical correlation using non-orthogonal configuration interaction singles and doubles (NOCISD),” Theoretical and Computational Chemistry
, chemRxiv:11369646.v1.
42.
P.-O.
Löwdin
,
Phys. Rev.
97
,
1474
(
1955
).
43.
P.-O.
Löwdin
,
Phys. Rev.
97
,
1490
(
1955
).
44.
J.
Rodriguez-Laguna
,
L. M.
Robledo
, and
J.
Dukelsky
,
Phys. Rev. A
101
,
012105
(
2020
).
45.
P.
Ring
and
P.
Schuck
,
The Nuclear Many-Body Problem
(
Springer-Verlag
,
1980
).
46.
H. G. A.
Burton
,
J. Chem. Phys.
154
,
144109
(
2021
).
47.
M.
Head-Gordon
,
P. E.
Maslen
, and
C. A.
White
,
J. Chem. Phys.
108
,
616
(
1998
).
48.
A. T.
Amos
and
G. G.
Hall
,
Proc. R. Soc. London, Ser. A
263
,
483
(
1961
).
49.
G. G.
Hall
,
Proc. R. Soc. London, Ser. A
205
,
541
(
1951
).
50.
S. C.
Leasure
and
G. G.
Balint-Kurti
,
Phys. Rev. A
31
,
2107
(
1985
).
51.
A.
Szabo
and
N. S.
Ostlund
,
Modern Quantum Chemistry
(
Dover Publications, Inc.
,
1989
).
52.
H. G. A.
Burton
, LIBGNME: C++ library for generalized nonorthognal matrix elements, https://github.com/hgaburton/libgnme.
53.
T. H.
Dunning
, Jr.
,
J. Chem. Phys.
90
,
1007
(
1989
).
54.
H.
Koch
,
A.
Sánchez de Merás
, and
T. B.
Pedersen
,
J. Chem. Phys.
118
,
9481
(
2003
).
55.
R. A.
Kendall
and
H. A.
Früchtl
,
Theor. Chim. Acta
97
,
158
(
1997
).