Amplitude determinant coupled cluster with pairwise doubles

In variational Monte Carlo, the energy expression is

which may be approximated by sampling a set of occupations $\xi$ from the wave function’s probability distribution ${|\left<{\vec{n}}|{\Psi}\right>|}^{2}/\left<{\Psi}|{\Psi}\right>$. The energy is then estimated as an average of local energies.

with $n_{s}$ being the number of Monte Carlo samples in $\xi$.

For our purposes, it is convenient to write the $ab$ $initio$ Hamiltonian as

where $\left(pq|{rs}\right)$ are the usual two-electron coulomb integrals in $\left(11|{22}\right)$ order, ${a}_{p}^{\dagger}$ (${a}_{\overline{p}}^{\dagger}$) and ${a}_{p}$(${a}_{\overline{p}}$) are creation and destruction operators of an $\alpha$($\beta$) electron in the $p$th spatial orbital, and $t_{pq}$ are modified one-electron integrals,

where $h_{pq}$ are the standard one-electron integrals.

From the energy and Hamiltonian expressions, we see that if an ansatz can evaluate $\left<{\vec{n}}|{\Psi}\right>$ efficiently, this is sufficient for use with VMC in Fock space. For the pCCD ansatz with RHF as reference, its cluster amplitudes can be write in matrix form

in which index $i,j,k$ represent occupied orbitals in the reference and $a,b,c$ represent virtual orbitals in the reference.

The amplitude expression for pCCD ansatz can be written as

in which $\hat{T}$ is the cluster operator

Suppose the occupation number vector is a quadruply excited configuration, ${\vec{n}}_{ij}^{ab}$, relative to the reference determinant. Then only the second order cluster operator term will not vanish. Therefore,

in which “$\mathrm{Perm}$” represents the permanent of a matrix and the factor 2 comes from the number of ways of permuting $i$, $j$, $a$, $b$ indices, which cancels the ${1}/{2}$ from the exponentials Taylor series. This result generalizes, with the amplitude of pCCD for a given occupation number vector given by a permanent of an ${n}_{ex}$/2 by ${n}_{ex}$/2 part of the cluster amplitude matrix, with ${n}_{ex}$ the excitation level in the occupation number vectorAyers:2013:ap1rog . For $N$th level excitations, there are $N$! ways to permute the indices, therefore the exponential’s 1$/N$! is canceled and no constants appear before the permanent.

Thus, the difficulty in evaluating a pCCD occupation number coefficient is equivalent to the evaluation of a permanent. However, there are no known polynomial cost methods for permanent evaluation, and so pCCD, like CC in general, is incompatible with VMC.

We propose a different, but related ansatz, in which coefficients are defined to be determinants rather than permanents. Using the $n_{ex}$/2 by $n_{ex}$/2 matrix

we define the wave function coefficient to be

These determinants can be evaluated through LU decomposition with $\mathcal{O}\left({n}_{ex}^{3}\right)$ cost, making this pADCCD ansatz compatible with VMC.

We should note that our choice of determinant instead of permanent forms a different approximation than pCCD, and since pCCD is not exact, using a determinant is not necessarily better or worse. Moreover, as we will prove later, size consistency still holds, and so one of CC’s most important properties is maintained. And since we make the lower triangle negative, the coefficients obtained from pADCCD match those of pCCD through quadruple excitations. Unlike pCCD, pADCCD cannot be written in the compact $\mathrm{exp}\left(\hat{T}\right)$ form, but with size consistency nonetheless maintained, this seems a minor inconvenience.

We will derive the energy expression for pADCCD in this subsection. Before considering our particular wave function, notice that the subset of Hamiltonian terms that contain only number and hole operators (i.e. ${a}_{i}^{\dagger}{a}_{i}$ and ${a}_{b}^{\dagger}{a}_{b}$) will add a wave-function-independent contribution to the local energy,

The remaining one-electron terms in the Hamiltonian of the type ${a}_{i}^{\dagger}{a}_{a}$ will not contribute to the local energy as it breaks seniority symmetry and pADCCD is strictly seniority zero. Likewise, the only two-electron terms in the Hamiltonian that will give a non-vanishing contribution to the local energy are the terms of the type ${a}_{i}^{\dagger}{a}_{a}{a}_{\overline{i}}^{\dagger}{a}_{\overline{a}}$. Therefore, the expression for the local energy is,

in which $\left|{\vec{n}}_{i\overline{i}}^{a\overline{a}}\right>$ is obtained by exciting two electrons ($\alpha$ and $\beta$) from the $ith$ occupied orbital to the $ath$ virtual orbital in $\left|\vec{n}\right>$.

In this subsection, we will prove that our pADCCD ansatz is strictly size consistent. The definition of size-consistency states that the energy calculated with two non-interacting systems A and B together as a “super system” should be equal to the sum of the energies of systems A and B calculated separately. To prove this, consider two non-interacting systems A and B. Then the coefficient matrix becomes block diagonal

since there are no “cross-excitations” between the two systems. Then due to the properties of determinants, we have

and it follows that,

Thus the energy expression becomes

and so pADCCD is size consistent.

Before discussing the variational minimization of the energy, we first lay the groundwork by developing efficient evaluations of terms that we call derivative ratios, which will make the optimization relatively simple to describe. Defining derivative notation $\left|{\Psi}^{x}\right>\equiv\partial\left|\Psi\right>/\partial{\mu}_{x}$, where ${\mu}_{x}$ is the xth wave function parameter, we first consider the “bare” derivative ratio

For amplitudes, these ratios are:

in which $\Theta\left(\vec{n}\right)$ is the inverse of $\tilde{T}\left(\vec{n}\right)$.

In addition to the bare derivative ratios, we also define the energy derivative ratios,

which are related to derivatives of the local energy by

Based on the expression of local energy, its derivative is:

and the evaluation of local energy derivatives are similar to that of bare derivative ratios.

We variationally optimize the energy of the pADCCD wave function using the linear method (LM)Nightingale:2001:linear_method ; UmrTouFilSorHen-PRL-07 . The LM works by repeatedly solving the Schrödinger equation in a special subspace of the full Hilbert space defined by the span of the wave function and its first derivatives. More precisely, we construct the generalized eigenvalue problem

where $\left|{\Psi}^{x}\right>$ and $\left|{\Psi}^{y}\right>$ are shorthand for derivatives of $\left|{\Psi}\right>$ with respect to the xth and yth wave function parameters ${\mu}_{x}$ and ${\mu}_{y}$, respectively, and $\left|{\Psi}^{0}\right>\equiv\left|\Psi\right>$. After solving this eigenvalue problem for $c$, one updates the parameters by

The Hamiltonian and overlap matrices are built by Monte Carlo sampling

Inspecting the above Monte Carlo approximations to the Hamiltonian and overlap matrices in the first derivative subspace makes clear that only derivative ratios $\mathcal{D}_{n}\left({\mu}_{x}\right)$ and $\mathcal{G}_{n}\left({\mu}_{x}\right)$ from Eqs. 17 and 19 are needed. Therefore, the LM can be applied efficiently to optimize the amplitudes by evaluating these ratios.

Like DOCI and pCCD, pADCCD is not invariant to the choice of orbital basis, since all open-shell determinants are neglected. Thus optimal orbitals need to be found to fully minimize the energy. In this section we introduce an orbital optimization method for pADCCD. First consider the one-body anti-Hermitian operator

which, when exponentiated, creates unitary orbital rotations (here $\sigma$ indexes spin).

Given this rotation operator, we can generalize the energy to be

We expand this energy to second order in $\vec{\kappa}$,

where we work at $\vec{\kappa}=0$ by transforming the basis in which we express the Hamiltonian (i.e. by transforming the one- and two-electron integrals).

We define the energy gradient and hessian,

both of which are functions of the one- and two-electron reduced density matrices

We evaluate both $\bm{\gamma}$ and $\bm{\Gamma}$ the same way we evaluate the energy.

It should be noted that the one-electron reduced density matrix $\gamma$ is diagonal in the basis in which we define the pairing. The two-electron reduced density matrix $\Gamma$ is also very sparse thanks to pADCCD being seniority zero. Detailed expressions for the orbital gradient and hessian can be found in the Appendix.

Using $\vec{\omega}$ and $\bm{A}$, we can minimize the energy with respect to orbital rotations by Newton-Raphson,

After each step, we update the MO-coefficients by

where $\bm{C}^{\left(0\right)}$ is the old MO-coefficients matrix, $\bm{C}$ is the new MO-coefficients matrix, and $\bm{U}$ is the MO-rotation matrix. We then recompute our MO integrals to work in the basis in which $\vec{\kappa}$ is again zero.

In each iteration, we first perform a linear method update for the cluster amplitudes. Then we evaluate the one- and two-electron reduced density matrices with the updated wave function. Finally, $\vec{\omega}$ and $\bm{A}$ are built and NR is performed to update the orbital basis.

Having presented our optimization method, we now analyze its cost. In each iteration, we need to loop over occupied and virtual orbitals to compute local energies, RDMs, and derivative ratios. As the evaluation of the relevant determinants scale as ${n}_{ex}^{3}$, in which $n_{ex}$ is the pair excitation level, the overall cost scales as ${n}_{s}{n}_{o}{n}_{u}{n}_{ex}^{3}$, in which $n_{s}$, $n_{o}$ and $n_{u}$ are the number of samples, number of occupied and unoccupied orbitals in reference determinant, respectively. Although this $\mathcal{O}\left(N^{6}\right)$ scaling looks much steeper than the $\mathcal{O}\left(N^{3}\right)$ scaling of pCCD, one should note the highly excited configurations are rarely sampled, and so in many molecules the ${n}_{ex}^{3}$ term in the scaling may behave more like a constant. In such a regime, pADCCD’s scaling may appear closer to $\mathcal{O}\left(N^{3}\right)$.

At the end of each iteration, we need to reset $\hat{\kappa}$ to 0 via a basis rotation. The one- and two-electron integrals needed to represent Hamiltonian in the new basis can be evaluated at an $\mathcal{O}\left({N}^{5}\right)$ cost. However, as the basis rotation is required only once per LM iteration, rather than once per sample, its cost is typically negligible compared to that of the sampling effort involved in optimizing cluster amplitudes.

pADCCD results were obtained using our own software for VMC in Hilbert space, with one- and two-electron integrals for the Hamiltonian taken from PySCFPyscf_brief . The full configuration interaction (FCI) results were obtained from MolproMOLPRO_brief and CISD results from Psi4Psi4:2012 . pCCD and DOCI results were kindly shared by Peter A. Limacherperson_comm . We froze N and O 1s orbitals at the RHF level. Sample size is taken to be 3.6$\times$10⁶. All statistical uncertainties were converged to less than 0.01eV in all cases.

Before showing our examples, we first check the size consistency of pADCCD by calculating the energy of up to 5 well separated H₂ molecules in a STO-3G basis. As shown in Table 1, the error per molecule does not grow with the increase of system size like in CISD. Instead, it remains zero within statistical uncertainty. This result demonstrates pADCCD’s size consistency, as we proved in section II.4.

We begin our results with a simple example, the LiH molecule in the cc-pVDZ basiscc-pvdz . The system has only two valence electrons, and both DOCI and pCCD deliver almost exact results compared to FCI. Due to our method’s similarity with pCCD, we also expect nearly exact results. Figure 1 shows our results at 17 bond lengths between 0.9$\mathring{A}$ and 4.4$\mathring{A}$. The non-parallelarity error (NPE) defined as the difference between the largest and smallest error with respect to FCI along the potential surface is about 2 mE_h, confirming our expectations.

We next turn our attention to the dissociation of hydrogen flouride in a 6-31G basis6-31g . Figure 2 shows the absolute energy of pADCCD, along with DOCI, pCCDperson_comm and FCI results. The results show that pADCCD, DOCI and pCCD are energetically almost indistinguishable. A close analysis shows that the NPE of pADCCD and pCCD with respect to FCI is about 17 and 18 mE_h, respectively. Unlike LiH, where seniority zero based wave functions are near exact, we see a large energy gap between all seniority zero based ansatzes and FCI. This is a remainder that while seniority zero wave functions are often effective for strong correlations, they do not capture all the details of weak correlation.

In Figure 2 we also plot the results using only RHF canonical orbitals and no further orbital optimization. As one can see, the results are quite poor when using RHF orbitals, especially when one stretches the bonds and the optimal orbitals become more and more localized.

Our next example is the symmetric double dissociation of H₂O, as shown in Figure 3. Again, pADCCD provides nearly identical energies compared to DOCI and pCCD, and the NPE with respect to FCI is about 19 mE_h for pADCCD and 16 mE_h for pCCD. The coincidence of the pADCCD with pCCD and DOCI is thus true not just for one pair of strongly correlated electrons, but for two pairs as well. However, we can still see from the plot that like pCCD and DOCI, a significant amount of dynamic correlation is clearly missing in all these methods.

In order to show the importance of orbital optimization, we plot optimized and RHF canonical orbitals for H₂O at bond length 2.2$\mathring{A}$ in Figure 4. It is very clear that at this stretched geometry, the optimized orbitals are much more localized than canonical orbitals. This localization is also seen in pCCDScuseria:2011:seniority . Indeed, the qualitative difference between optimized orbitals and canonical orbitals emphasizes the necessity of orbital optimization.

Our final example is the dissociation of N₂. As Figure 5 reveals, the difference between pCCD and pADCCD is more noticeable in N₂, which is to be extected now that hextuples (the first excitation level at which the ansatz forms are different) are needed for a qualitatively correct description of the dissociation. Indeed, we find that disabling hextuple and higher excitations raises the pADCCD energy by 0.5eV. Although their NPEs now differ noticeably, 56 mE_h for pCCD versus 65 mE_h for pADCCD, they are of the same order of magnitude and neither are close to quantitative. Achieving a more quantitative accuracy clearly requires a more flexible cluster expansion, but, as is well known , this route leads to qualitatively incorrect variational violations when pursued in a traditional CC approach. As variational violations are not possible in a VMC-based approach, it will be interesting in future to investigate more flexible expansions within the amplitude determinant framework.