Geminal-Based Configuration Interaction
Thomas M. Henderson Gustavo E. Scuseria



Department of Chemistry, Rice University, Houston, TX 77005-1892
Department of Physics and Astronomy, Rice University, Houston, TX 77005-1892


Department of Chemistry, Rice University, Houston, TX 77005-1892
Department of Physics and Astronomy, Rice University, Houston, TX 77005-1892



The antisymmetrized geminal power (AGP) wave function has a long history and considerable conceptual appeal, but in many situations its accuracy is wanting.
Here, we consider a form of configuration interaction (CI) based upon the AGP wave function and taking advantage of its killing operators to construct an excitation manifold.
Our geminal CI reduces to standard single-determinant–based CI in the limit in which AGP reduces to a single determinant.
It substantially improves upon AGP in the reduced BCS Hamiltonian, which serves as a prototype for the kinds of strong pairing correlations relevant in Bardeen-Cooper-Schrieffer–style superconductivity.
Moreover, our geminal CI naturally generalizes to add correlation to more general geminal-based wave functions than AGP.



Geminal-Based Configuration Interaction
    Gustavo E. Scuseria
    December 30, 2023
=======================================





§ INTRODUCTION

Suppose one has an approximate wave function |Φ⟩ which one wishes to improve.  A simple route toward doing so is to write a new wave function

    |χ⟩ = |Φ⟩ + ∑ c_μ |Φ_μ⟩

where the states |Φ_μ⟩ are linearly-independent of |Φ⟩ and the coefficients c_μ are parameters to be optimized.  By the variational principle, |χ⟩ almost certainly has lower energy than does |Φ⟩; at the very least, it cannot be higher.

The trick, of course, is to make a clever choice of states |Φ_μ⟩.  We would like |Φ_μ⟩ to be linearly independent of |Φ⟩ and we would like Hamiltonian and matrix elements between |Φ⟩ and |Φ_μ⟩ and between states |Φ_μ⟩ and |Φ_μ^'⟩ to be readily evaluated.  If |Φ⟩ has a symmetry of the exact wave function, we would like |Φ_μ⟩ to do likewise.  And of course we would like to choose states |Φ_μ⟩ which provide the most improvement to |Φ⟩ for the least amount of work.

One convenient way to choose the states |Φ_μ⟩ is to look for operators K_μ(Φ) whose action on |Φ⟩ is to annihilate it:

    K_μ(Φ) |Φ⟩ = 0,

which we shall refer to as killing operators of |Φ⟩.  Then the subset of all operators K_μ^†(Φ) which are not themselves killing operators provides a convenient way to construct the states |Φ_μ⟩, which we can create as simply

    |Φ_μ⟩ = K_μ^†(Φ)   |Φ⟩.

These states are automatically orthogonal to |Φ⟩ (and therefore linearly independent of it).  Further, matrix elements can be readily evaluated because we can adopt a normal order in which all operators K_μ^† appear adjacent to the bra and all operators K_μ appear adjacent to the ket; by reordering the various operators in a given matrix element with the aid of commutation or anticommutation rules, we can extract the non-zero portion of the matrix element as that piece which contains neither K nor K^†, and we are left with expectation values only with respect to |Φ⟩ which we presume can be easily computed.  Additionally, if we wish to enforce that the state |Φ_μ⟩, like the state |Φ⟩, is an eigenstate of some symmetry operator S with eigenvalue s, then we merely need to further limit our choice of killing operators to those which commute with S.  The (hermitian adjoints of) the killing operators thus have many desirable properties which make them well-suited to add corrections to |Φ⟩.

This observation is by no means new.  Indeed, the standard particle-hole construction in configuration interaction (CI) or coupled cluster theory proceeds entirely along these lines, and similar concerns motivate development in multi-reference theories.  In this manuscript we seek to follow the same path, using the antisymmetrized geminal power (AGP) wave function<cit.> as the reference |Φ⟩.  To do so, we first need to construct killing operators of the AGP wave function; this, in turn, requires us to say a few words about AGP.




§ AGP AND ITS KILLING OPERATORS

The AGP wave function can be written as

    |AGP⟩ = (Γ^†)^N |-⟩

where |-⟩ is the physical vacuum, N is the number of electron pairs in the system, and Γ^† creates a geminal:

    Γ^† = ∑η_p   P_p^†.

Here, P_p^† creates an electron in level p and another in a level p̅ paired with p:

    P_p^† = c^†_p   c^†_p̅.


The AGP wave function has a symmetry known as seniority in which levels p and p̅ are either both occupied or both empty in every determinant in the expansion of the AGP wave function.  General Hamiltonians, of course, do not have seniority as a symmetry, so typically we will need killing operators which do not necessarily conserve seniority.  We shall have more to say on this point later (and indeed, such operators for AGP appeared in the literature many years ago<cit.>), but for now we will limit our attention to seniority-conserving killing operators, which we may construct by making use of the operators P_p^† as well as 


    P_p    = c_p̅  c_p,
    
    
    N_p    = c^†_p   c_p + c^†_p̅  c_p̅


which also preserve seniority.

The only one-body operators we have at hand which preserve both seniority and particle number are the number operators N_p; as they are Hermitian, if a linear combination of the N_p's is a killer, so is its adjoint (we assume that the parameters η_p and the coefficients in any linear combination of operators intended for use as a killing operator are all real, which guarantees that the wave function also respects complex conjugation symmetry).  Thus, the relevant seniority- and number-conserving killing operators of the AGP wave function are two-body or higher.  We can write these killing operators as

    K_pq = a   P_p^†  P_q + b   P_q^†  P_p + c   N_p + d   N_q + e   N_p   N_q

where we assume p  q (for p = q, the operator is Hermitian and thus of no use to us).

A simple way to guarantee that K_pq is a killing operator is to adjust the coefficients a, b, c, d, and e such that K_pq commutes with Γ^†; if this is the case, then

    K_pq |AGP⟩ = (Γ^†)^N   K_pq |-⟩

and as K_pq annihilates the physical vacuum, it would then also annihilate the AGP state.

Using the commutation relations


    [P_p,P_q^†]    = δ_pq (1 - N_p),
    
    
    [N_p,P_q^†]    = 2  δ_pq  P_q^†,


one can verify that up to an overall multiplicative factor, the unique operator K_pq which has these properties is

    K_pq   = η_p^2   P_p^†  P_q + η_q^2   P_q^†  P_p
    
       + 1/2 η_p  η_q  (N_p   N_q - N_p - N_q).

Note that K_pq = K_qp.  For η_p ≠η_q, we have K_pq^†≠ K_pq, so we can hope to write a sort of geminal configuration interaction (GCI) state as

    |GCI⟩ = |AGP⟩ + ∑_p > q
    η_p ≠η_q c_pq  K_pq^†  |AGP⟩.




This is the wave function we intend to investigate in this work.  We will limit our attention to AGP wave functions which respect spin symmetry, so that the orbitals p and p̅ are the ↑-spin and ↓-spin spinorbitals for the same spatial orbital.  Further, because we have chosen to consider only the seniority-conserving kiling operators, we should limit ourselves to Hamiltonians for which seniority is a symmetry; accordingly, we will specialize to the pairing or reduced BCS Hamiltonian, which can be written entirely in terms of the seniority zero generators P, P^†, and N and for which the exact energy can be obtained via a form of Bethe ansatz which simply requires solving the nonlinear Richardson equations.<cit.>  Note that even for general Hamiltonians, only the portion of the Hamiltonian which itself conserves seniority contributes to AGP expectation values involving H, K, and K^†.




§ THE PAIRING HAMILTONIAN

The pairing Hamiltonian can be written as

    H = ∑_p ϵ_p   N_p - G  ∑_pq P_p^†  P_q.

While other choices are possible, we take ϵ_p = p and G>0.  Note that for positive G the Hamiltonian has an attractive two-body interacton rather than the more familiar repulsive two-body interaction.  The mean-field symmetry which spontaneously breaks is, accordingly, number symmetry rather than spin symmetry, giving rise to the well-known Bardeen-Cooper-Schrieffer (BCS) wave function.<cit.>  Note also that the interaction is infinite in range.

This Hamiltonian was introduced as a phenomenological model to describe pairing correlations in which electron pairs interact with the hole pairs they have left behind, but has also been used to model ultrasmall superconducting grains.<cit.>



Remarkably, the Hamiltonian can be solved exactly.  Its eigenstates can be written as


    |Ψ⟩   = ∏Γ_μ^† |-⟩,
    
    Γ_μ^†   = ∑_p 1/2  ϵ_p - R_μ  P_p^†,


where the energy is given by the sum of the rapidities R_μ, which satisfy a set of coupled nonlinear equations.  Note that the exact wave function is simply an antisymmetric product of non-orthogonal geminals.

Although the exact solution is readily available, this Hamiltonian is very difficult for conventional electronic structure methods such as single-reference coupled cluster theory<cit.> and even extended coupled cluster.<cit.>  This can readily be seen by Fig. <ref>, which shows the fraction

    %E_c,HF = E - E_HF/E_exact - E_HF

of the exact correlation energy with respect to Hartree-Fock recovered as a function of G/G_c. where G_c refers to the value of G at which the mean-field undergoes spontaneous symmetry breaking toward a number-broken mean-field solution.  But while sophisticated coupled cluster approaches fail for this problem, AGP provides good albeit imperfect results for almost all G and is an eigenstate in the G →∞ limit.[The manner of plotting in Fig. <ref> suggests errors in AGP for small G, but in fact AGP is exact for G → 0; Fig. <ref> merely reflects the fact that AGP delivers the wrong 𝒪(G^2) contribution to the total energy]  Note that in the repulsive case the situation is very different, and traditional coupled cluster performs exceptionally well.  The attractive pairing Hamiltonian thus forms an interesting test case for our geminal CI.  Our hope is that if we can correctly describe this exactly solvable model, we can also correctly describe related but more physical problems for which no exact solution exists but which pose similar computational challenges.




§ RESULTS

To test the accuracy of our geminal CI ansatz, we return to the half-filled 16-site pairing Hamiltonian, which is shown in Fig. <ref>, the left panel of which can be directly compared to Fig. <ref>.  Although the geminal CI is not exact, it has successfully remedied almost all of the deficiencies of AGP.  This is made much clearer by examination of the right panel of Fig. <ref>, which shows the fraction %E_c,AGP of correlation with respect to AGP recovered by the geminal CI wave function and shows that GCI removes, in this case, more than 99% of the energy error made by AGP.

The right panel also reveals a curious (and unpleasant) feature of the geminal CI.  While for small G the AGP reference becomes more like a single determinant (see below) for which we expect the number of double excitations to be N   (M-N) where N is the number of pairs and M is the number of levels, the number of killing operators we construct is M   (M-1)/2.  Several must be linearly-dependent for small G, and these linearly-dependent states must be removed from the diagonalization.  This, we expected.  What was unexpected is that we always encounter one linearly-dependent mode even for larger G, and this appears to be true independent of the number of pairs N and number of levels M.  Put differently, there always appears to be at least one linear combination of operators K_pq^† which is a killing operator of AGP and which must therefore be excluded.  We do not at present have a better prescription for finding this mode than by diagonalizing the metric ⟨ K_pq  K_rs^†⟩ and eliminating the mode with zero eigenvalue.  Note that this zero eigenvalue is not simply a small eigenvalue, but rather one which vanishes to the double precision in which we have implemented our geminal CI.



To get a sense of how our method performs as we adjust system size, Fig. <ref> shows the fraction of correlation energy recovered in several different pairing Hamiltonians.  Clearly our results deteriorate somewhat as system size increases, but the accuracy remains uniformly excellent.  Results away from half-filling (not shown) are similar.  Neither AGP nor truncated CI are fully extensive methods, and we ought not expect geminal CI to be extensive either.  The pairing Hamiltonian, unfortunately, is not a good place to test the error in extensivity, because the exact energy is not linear in system size in the first place since the interaction has infinite range.  Nonetheless, as the right panel of Fig. <ref> makes clear, a large portion of the size-dependence in the geminal CI is due to the underlying AGP rather than to the geminal CI per se.



Although obtaining accurate energies is perhaps the chief priority, we would also like similar accuracy for the calculation of properties (or, in other words, we would like to get both the energy and the wave function right).  To assess this, Fig. <ref> shows the errors in the occupation numbers n_p = ⟨ N_p ⟩ for the 16-site, half-filled pairing Hamiltonian at G ∼ 2   G_c.  This value is large enough that we expect to see occupations differ appreciably from 0 and 2 (the mean-field result) without being too large.  Clearly the geminal CI significantly reduces the error in the occupation numbers, and thus in typical single-particle properties.  We should note that the plot has a symmetry around n_p = 1 which is just a consequence of a symmetry of the half-filled pairing Hamiltonian.<cit.>




§ DISCUSSION

We have seen that for the pairing Hamiltonian, our geminal CI is remarkably accurate.  Here, we wish to discuss the extension to more realistic Hamiltonians which do not have seniority as a symmetry.

For such Hamiltonians, the relevant killing operators are one-body in nature, and can be written as

    D_pq = η_p   E^p_q - η_q   E^q_p

for p ≠ q in terms of operators

    E^p_q = c_p^†  c_q + c_p̅^†  c_q̅

which become the standard spin-integrated unitary group generators when levels p and p̅ are the ↑- and ↓-spinorbitals corresponding to the same spatial orbital.  Killers of AGP were first discussed in quantum chemistry,<cit.> and recently applied by Dukelsky and collaborators to superfluid nuclei,<cit.> albeit not in the same D_pq form as above.  We may think of the operators D, D^†, and N as a basis of one-body operators and so in principle can use them to write a normal order of the form D^†  N   D.

Our killing operators K_pq are essentially just 1/2   D_pq^2; while not strictly identical, K_pq and 1/2   D_pq^2 act identically on seniority-zero states.  For problems such as the molecular Hamiltonian, we can generalize our geminal CI to the form

    |GCI⟩ = (1    + ∑ c_pq  D_pq^†
       + 1/2 ∑ c_pqrs {D_pq^†,D_rs^†} + …) |AGP⟩.

The expectation value of the Hamiltonian with respect to this wave function can readily be taken and minimized with respect to the variational parameters, and we shall report on its applications in due time.  The anticommutator in the quadratic term is present to ensure that it corresponds to a genuinely two-body operator (because the commutator of D_pq^† with D^†_rs is a non-zero one-body operator).  

Let us say a few words about the single-determinant limit.  If we have only N non-zero parameters η_p, then the AGP wave function reduces to a single determinant.  Using the conventional notation in which a and b index unoccupied orbitals with η = 0 and i and j index occupied orbitals with η≠ 0, then


    D_ab   = 0,
    
    
    D_ia   = η_i   E^i_a,
    
    
    D_ij   = η_i   E^i_j - η_j   E^j_i.


The adjoint operators D_ij^† annihilate the reference determinant because i ≠ j, while the adjoint operators D_ia^† are, up to an overall factor, the usual particle-hole excitation operators.  Thus, in the single-determinant limit of AGP, our geminal CI reduces to standard single-reference configuration interaction when using the D_pq killers (and to the CI analog of pair coupled cluster doubles when using the K_pq killers).



Our geminal CI bears close resemblance to some methods which merge correlation and symmetry projection, essentially because AGP is a symmetry-projected mean-field method (AGP is equivalent to number-projected BCS).  Although the two methods are distinct, geminal CI is particularly closely related to the projected configuration interaction of Ripoche et al.<cit.> (and also to a spin-projected analog due to Tshuchimochi and Ten-no<cit.>).  A coupled-cluster generalization of our geminal CI would likewise approximate symmetry-projected coupled cluster theory.<cit.>  These methods write a correlated wave function which uses the killing operators of BCS and then act a number projector on the correlated state.  In contrast, our method works directly with those operators which annihilate the projected mean-field state.  Figure <ref> shows that for a small pairing model, geminal CI provides accuracy similar to symmetry-projected coupled cluster with singles and doubles without needing to worry about the highly technical details of how one might carry out the symmetry projection of a correlated state.

Finally, a word on the name “geminal CI”.  We have chosen this name because while we have considered only the AGP case, the ideas discussed in this manuscript generalize to other geminal wave functions.  Consider, for example, a geminal wave function

    |Ψ⟩ = ∏(Γ_μ^†)^N_μ |-⟩

where the geminals created by the different Γ_μ^† are strongly orthogonal (i.e. are expanded in mutually orthogonal single-particle orbital spaces):

    Γ_μ^† = ∑_p ∈μη_p   P_p^†.

This geminal wave function is a kind of cluster mean-field state.<cit.>  It includes AGP on the one hand (when there is a single geminal so that ∏(Γ_μ^†)^N_μ = (Γ^†)^N) and the antisymmetrized product of strongly orthogonal geminals<cit.> (APSG) on the other (when the various N_μ are all equal to 1).  Such wave functions are annihilated by our killing operators K_pq and D_pq and thus form a perfectly valid reference state for our geminal CI.  An APSG-based version of our geminal CI would be related to both Kállay and Surján's APSG-based CI method<cit.> and Goddard's generalized valence bond CI.<cit.>  While we have limited our attention in this manuscript to the AGP case, the theory is considerably more general.


This work was supported by the U.S. National Science Foundation under Grant CHE-1762320. G.E.S. is a Welch Foundation Chair (C-0036).  We thank Jorge Dukelsky for providing his AGP code and Armin Khamoshi for providing code to evaluate AGP density matrices.




