Geminal-Based Configuration Interaction

TL;DR

This paper introduces a geminal-based configuration interaction method that enhances the accuracy of AGP wave functions, especially for systems with strong pairing correlations, by constructing an excitation manifold that generalizes standard CI.

Contribution

It develops a geminal CI approach that improves upon AGP and can incorporate more general geminal wave functions, advancing the modeling of strongly correlated systems.

Findings

Geminal CI outperforms AGP in the reduced BCS Hamiltonian.

The method reduces to standard CI when AGP simplifies to a single determinant.

It naturally extends to more general geminal wave functions.

Abstract

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.

Equations

$$|\chi ⟩=|\Phi ⟩+\sum c_{\mu }|{\Phi }_{\mu }⟩$$

$$K_{\mu }(\Phi )|\Phi ⟩=0,$$

$$|{\Phi }_{\mu }⟩=K_{\mu }^{†}(\Phi )|\Phi ⟩.$$

$$|\mathrm{AGP}⟩={\left({\Gamma }^{†}\right)}^N|-⟩$$

$${\Gamma }^{†}=\sum {\eta }_p{P}_p^{†}.$$

$$P_p^{†}=c_p^{†}{c}_{\bar{p}}^{†}.$$

$$P_p$$

$$N_p$$

$$K_{pq}=a{P}_p^{†}{P}_q+b{P}_q^{†}{P}_p+c{N}_p+d{N}_q+e{N}_p{N}_q$$

$$K_{pq}|\mathrm{AGP}⟩={\left({\Gamma }^{†}\right)}^N{K}_{pq}|-⟩$$

$$[P_p,P_q^{†}]$$

$$[N_p,P_q^{†}]$$

$$K_{pq}$$

$$+\frac{1}{2}{\eta }_p{\eta }_q\left(N_p{N}_q-N_p-N_q\right).$$

$$|\mathrm{GCI}⟩=|\mathrm{AGP}⟩+\sum _{\begin{array}{c}p>q\\ {\eta }_p\ne {\eta }_q\end{array}}{c}_{pq}{K}_{pq}^{†}|\mathrm{AGP}⟩.$$

$$H=\sum _p{ϵ}_p{N}_p-G\sum _{pq}{P}_p^{†}{P}_q.$$

$$|\Psi ⟩$$

$${\Gamma }_{\mu }^{†}$$

$$\%E_{c,\mathrm{HF}}=\frac{E-E_{\mathrm{HF}}}{E_{\mathrm{exact}}-E_{\mathrm{HF}}}$$

$$D_{pq}={\eta }_p{E}_q^p-{\eta }_q{E}_p^q$$

$$E_q^p=c_p^{†}{c}_q+c_{\bar{p}}^{†}{c}_{\bar{q}}$$

$$\displaystyle|\mathrm{GCI}\rangle=\big{(}1$$

$$\displaystyle+\frac{1}{2}\,\sum c_{pqrs}\,\{D_{pq}^{\dagger},D_{rs}^{\dagger}\}+\ldots\big{)}|\mathrm{AGP}\rangle.$$

$$D_{ab}$$

$$D_{ia}$$

$$D_{ij}$$

$$|\Psi ⟩=\prod {\left({\Gamma }_{\mu }^{†}\right)}^{N_{\mu }}|-⟩$$

$${\Gamma }_{\mu }^{†}=\sum _{p\in \mu }{\eta }_p{P}_p^{†}.$$

Full text

I Introduction

Suppose one has an approximate wave function $|\Phi\rangle$ which one wishes to improve. A simple route toward doing so is to write a new wave function

[TABLE]

where the states $|\Phi_{\mu}\rangle$ are linearly-independent of $|\Phi\rangle$ and the coefficients $c_{\mu}$ are parameters to be optimized. By the variational principle, $|\chi\rangle$ almost certainly has lower energy than does $|\Phi\rangle$; at the very least, it cannot be higher.

The trick, of course, is to make a clever choice of states $|\Phi_{\mu}\rangle$. We would like $|\Phi_{\mu}\rangle$ to be linearly independent of $|\Phi\rangle$ and we would like Hamiltonian and matrix elements between $|\Phi\rangle$ and $|\Phi_{\mu}\rangle$ and between states $|\Phi_{\mu}\rangle$ and $|\Phi_{\mu^{\prime}}\rangle$ to be readily evaluated. If $|\Phi\rangle$ has a symmetry of the exact wave function, we would like $|\Phi_{\mu}\rangle$ to do likewise. And of course we would like to choose states $|\Phi_{\mu}\rangle$ which provide the most improvement to $|\Phi\rangle$ for the least amount of work.

One convenient way to choose the states $|\Phi_{\mu}\rangle$ is to look for operators $K_{\mu}(\Phi)$ whose action on $|\Phi\rangle$ is to annihilate it:

[TABLE]

which we shall refer to as killing operators of $|\Phi\rangle$. Then the subset of all operators $K_{\mu}^{\dagger}(\Phi)$ which are not themselves killing operators provides a convenient way to construct the states $|\Phi_{\mu}\rangle$, which we can create as simply

[TABLE]

These states are automatically orthogonal to $|\Phi\rangle$ (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_{\mu}^{\dagger}$ appear adjacent to the bra and all operators $K_{\mu}$ 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^{\dagger}$, and we are left with expectation values only with respect to $|\Phi\rangle$ which we presume can be easily computed. Additionally, if we wish to enforce that the state $|\Phi_{\mu}\rangle$, like the state $|\Phi\rangle$, 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 $|\Phi\rangle$.

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 functionColeman (1965) as the reference $|\Phi\rangle$. 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.

II AGP And Its Killing Operators

The AGP wave function can be written as

[TABLE]

where $|-\rangle$ is the physical vacuum, $N$ is the number of electron pairs in the system, and $\Gamma^{\dagger}$ creates a geminal:

[TABLE]

Here, $P_{p}^{\dagger}$ creates an electron in level $p$ and another in a level $\bar{p}$ paired with $p$:

[TABLE]

The AGP wave function has a symmetry known as seniority in which levels $p$ and $\bar{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 agoWeiner and Goscinski (1983)), 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}^{\dagger}$ as well as

[TABLE]

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 $\eta_{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

[TABLE]

where we assume $p\neq 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 $\Gamma^{\dagger}$; if this is the case, then

[TABLE]

and as $K_{pq}$ annihilates the physical vacuum, it would then also annihilate the AGP state.

Using the commutation relations

[TABLE]

one can verify that up to an overall multiplicative factor, the unique operator $K_{pq}$ which has these properties is

[TABLE]

Note that $K_{pq}=K_{qp}$. For $\eta_{p}\neq\eta_{q}$, we have $K_{pq}^{\dagger}\neq K_{pq}$, so we can hope to write a sort of geminal configuration interaction (GCI) state as

[TABLE]

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 $\bar{p}$ are the $\uparrow$-spin and $\downarrow$-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^{\dagger}$, 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.Richardson (1963); Dukelsky and Pittel (2004) 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^{\dagger}$.

III The Pairing Hamiltonian

The pairing Hamiltonian can be written as

[TABLE]

While other choices are possible, we take $\epsilon_{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.Bardeen, Cooper, and Scrieffer (1957) 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.Braun and von Delft (1998); Dukelsky and Sierra (1999); Sierra et al. (2000)

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

[TABLE]

where the energy is given by the sum of the rapidities $R_{\mu}$, 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 theoryDukelsky et al. (2003); Henderson et al. (2014) and even extended coupled cluster.Henderson, Bulik, and Scuseria (2015) This can readily be seen by Fig. 1, which shows the fraction

[TABLE]

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\to\infty$ limit.111The manner of plotting in Fig. 1 suggests errors in AGP for small $G$, but in fact AGP is exact for $G\to 0$; Fig. 1 merely reflects the fact that AGP delivers the wrong $\mathcal{O}(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.

IV 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. 2, the left panel of which can be directly compared to Fig. 1. 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. 2, which shows the fraction $\%E_{c,\mathrm{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}^{\dagger}$ 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 $\langle K_{pq}\,K_{rs}^{\dagger}\rangle$ 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. 3 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. 3 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. 4 shows the errors in the occupation numbers $n_{p}=\langle N_{p}\rangle$ for the 16-site, half-filled pairing Hamiltonian at $G\sim 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.Hirsch et al. (2002)

V 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

[TABLE]

for $p\neq q$ in terms of operators

[TABLE]

which become the standard spin-integrated unitary group generators when levels $p$ and $\bar{p}$ are the $\uparrow$- and $\downarrow$-spinorbitals corresponding to the same spatial orbital. Killers of AGP were first discussed in quantum chemistry,Weiner and Goscinski (1983) and recently applied by Dukelsky and collaborators to superfluid nuclei,Dukelsky et al. albeit not in the same $D_{pq}$ form as above. We may think of the operators $D$, $D^{\dagger}$, 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^{\dagger}\,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

[TABLE]

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}^{\dagger}$ with $D^{\dagger}_{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 $\eta_{p}$, then the AGP wave function reduces to a single determinant. Using the conventional notation in which $a$ and $b$ index unoccupied orbitals with $\eta=0$ and $i$ and $j$ index occupied orbitals with $\eta\neq 0$, then

[TABLE]

The adjoint operators $D_{ij}^{\dagger}$ annihilate the reference determinant because $i\neq j$, while the adjoint operators $D_{ia}^{\dagger}$ 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.Ripoche et al. (2017) (and also to a spin-projected analog due to Tshuchimochi and Ten-noTsuchimochi and Ten-no (2016a, b)). A coupled-cluster generalization of our geminal CI would likewise approximate symmetry-projected coupled cluster theory.Duguet (2015); Duguet and Signoracci (2017); Tsuchimochi and Ten-no (2017); Qiu et al. (2017, 2019) 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 5 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

[TABLE]

where the geminals created by the different $\Gamma_{\mu}^{\dagger}$ are strongly orthogonal (i.e. are expanded in mutually orthogonal single-particle orbital spaces):

[TABLE]

This geminal wave function is a kind of cluster mean-field state.Jiménez-Hoyos and Scuseria (2015) It includes AGP on the one hand (when there is a single geminal so that $\prod\left(\Gamma_{\mu}^{\dagger}\right)^{N_{\mu}}=\left(\Gamma^{\dagger}\right)^{N}$) and the antisymmetrized product of strongly orthogonal geminalsSurján (1999); Rassolov (2002); Surján et al. (2012); Limacher et al. (2013); Pernal (2014) (APSG) on the other (when the various $N_{\mu}$ 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 methodKállay and Surján (1999) and Goddard’s generalized valence bond CI.Carter and Goddard III (1988); Faglioni and Goddard III (1999) While we have limited our attention in this manuscript to the AGP case, the theory is considerably more general.

Acknowledgements.

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.