A well-scaling natural orbital theory

The NSO Basis.

We consider time-reversal invariant saturated systems with nondegenerate, singlet ground states. The inverse approach (6⇓–8) starts from the N-representability conditions on the 2-DM. Instead, we take a forward approach: We introduce a specific form for the trial wavefunction and derive the 2-DM explicitly. Our starting point is that of conventional FCI, except that our one-particle basis is the complete set of NSOs of the trial function Ψ, ψk(x)=ϕk(r)χk(σ), with r space and σ-spin coordinates. The NOs ϕk(r) can be real and are independent of the spin function. The complete set of N-electron orthonormal SDs Φn(x1,x2,⋯,xN), n=k1,k2,⋯,kN, formed from its NSOs supports representation of any trial wavefunction Ψ(x1,x2,⋯,xN) as the expansionΨ(x1,⋯,xN)=∑n CnΦn(x1,⋯,xN).[1]As the ground-state wave function can be chosen to be real, so can the trial functions and the normalized Cn (∑nCn2=1). The NSOs vary with the trial function or the coefficients in the search for the ground state. The combinatorial complexity of determining the ground-state energy by variation of the Cn is composed of the separate combinatorial complexities of the signs and magnitudes of the coefficients Cn. We use distinct reductive approximations for their signs and magnitudes. The signs depend on the sign convention chosen for the SDs. We use the Leibniz form for the SDs,Φn(x1,⋯,xN)=1N!∑p sgn{Pp}Ppψk1(x1)⋯ψkN(xN).[2]The sum is over the elements of the symmetric group of order N, the permutations Pp. The sign of Φn is fixed by the ordering k1<k2<⋯<kN in the product of the NSOs ψki in [2]. The SDs and their coefficients can then be specified by listing the NSOs occupied in the SDs, i.e., by the index n.

The 1-DM, the Orthogonality Constraint, the Pair-Difference Constraint.

The 1-DM of Ψ,ρ(x′,x)=N∫dx2⋯dxNΨ(x′,x2,⋯,xN)Ψ(x,x2,⋯,xN),[3]becomesρ(x′,x)=∑i≠j[∑m∌i,jCi,mCj,m]ψi(x′)ψj(x)[4]after [1] and [2] are inserted into [3]. In [4] the subindex m specifies the N−1 NSOs present in Φi,m and Φj,m, excluding ψi and ψj. As the ψ are the NSOs of Ψ, the eigenfunctions of ρ(x′,x), the bracketed quantity in [4] must vanish for i≠j. Regard the coefficients Ci,m and Cj,m as the components of vectors Ci and Cj and the bracket as their scalar product Ci⋅Cj, which must vanish. There are two realizations of this orthogonality constraint (OC). In the first, and most general form, the presence of Φi,m in Ψ does not exclude the presence of Φj,m. The individual terms in the scalar product need not vanish. The second, the pair-difference constraint (PDC), is a special case of the OC, in which the presence of Φi,m excludes Φj,m so that either Ci,m or Cj,m is zero for each m, and the sum vanishes term by term. Under the PDC, those Φ present in the expansion of Ψ must differ from one another by at least two NSOs. The OC is necessary and sufficient for N representability, whereas the PDC is only sufficient. We impose the PDC on the Φ as a simplifying variational approximation. The PDC was proved to be well-satisfied in the FCI result for H₈ using the minimal basis set STO-6G (Slater-type basis functions expressed with six contracted Gaussian functions).

Under the OC or PDC, ρ(x′,x) takes the diagonal formρ(x′,x)=∑k p1(k) ψk(x′)ψk(x).[5]Here, the p1(k)=∑nCn2 νk,n, where νk,n=1 if k∈n and 0 otherwise, are the eigenvalues of ρ(x′,x), the occupation numbers or occupation probabilities (1-OP) of its eigenfunctions ψk. They satisfy the necessary and sufficient conditions 0≤p1(k)≤1 and ∑kp1(k)=N. In general, only M>N occupation numbers p1(k) are nonnegligible, and only the corresponding active NSOs need be included in the representation of any trial function, providing a natural cutoff. The 1-DM is thus of algebraic complexity in the 1-OPs and the NSOs.

The 2-DM, the Sign Conjecture, the ξ-Approximation.

The 2-DM of Ψ,π(x1′x2′;x1x2)=N(N−1)∫dx3⋯dxN×Ψ(x1′,x2′,x3,⋯,xN)Ψ(x1,x2,x3,⋯,xN),becomesπ(x1′x2′;x1x2)=∑i<i′,j<j′,mi,i′,j,j′∉mCii′mCjj′m(ψi(x1′)ψi′(x2′)−ψi′(x1′)ψi(x2′))(ψj(x1)ψj′(x2)−ψj′(x1)ψj(x2)).[6]π separates into a part πd diagonal in the indices, i.e., with ii′=jj′, and an off-diagonal part, πod, with ii′≠jj′:πd(x1′x2′;x1x2)=12∑i≠j p11(ij)(ψi(x1′)ψj(x2′)−ψj(x1′)ψi(x2′))(ψi(x1)ψj(x2)−ψj(x1)ψi(x2)),[7]πod(x1′x2′;x1x2)=∑i<i′≠j<j′,mi,i′,j,j′∉mCii′mCjj′m(ψi(x1′)ψi′(x2′)−ψi′(x1′)ψi(x2′))(ψj(x1)ψj′(x2)−ψj′(x1)ψj(x2)).[8]Electron correlation is expressed through πod. The analogous off-diagonal part of ρ(x′,x) is suppressed by the OC, an advantage of the NSO basis. Note that the PDC has eliminated 3-index terms from πod in [8].

The p11(ij)= ∑nCn2 νi,nνj,n in πd are joint two-state occupation probabilities (2-OPs). Mazziotti has reported (10, 15) necessary and sufficient conditions on the 2-OPs that arise from the positivity conditions on the q-OPs, i.e., the p11⋯1(i1,i2,⋯,iq)=∑nCn2 νi1,nνi2,n⋯νiq,n, at any order 2≤q≤N. These conditions derive from the positivity conditions (15) on the diagonal elements of the 2-DM (16). Limiting ourselves to the (2,2) and (2,3) conditions, the following conditions for the 2-OPs hold:sup(p1(i)+p1(j)−1,0)≤p11(ij) ≤ p1(<),[9]sup(p1(i)+p1(j)+p1(k)−1,0)≤p11(ij)+p11(ik)+p11(jk).[10]p1(<) is the lesser of p1(i) and p1(j). In addition, the sum rule∑j(≠i) p11(ij)=(N−1)p1(i),[11]must be satisfied. Conditions 9–11 were first established by Weinhold and Bright Wilson (17). They are necessary but not sufficient conditions for N representability (7, 10, 16, 18). Establishing a complete set of conditions is quantum Merlin–Arthur (QMA)-hard in N because the number of (2,q) positivity conditions increases combinatorially with increasing q≤N. Fortunately, numerical calculations on atoms and molecules indicate that sufficiently accurate lower-bound ground-state energies often result by imposing (2,q)-positivity conditions with q≤3 (9, 19). This suggests that even in the most difficult situations, fermionic problems in atoms and molecules should require only a finite and small set of positivity conditions. Here we shall limit ourselves to conditions 9–11, as we found in our numerical calculations that they are sufficient to produce accurate lower-bounds. If higher-order conditions were found to be necessary, it would not be hard for us to add a few more.

The πd of [7] contains only 2-OPs and products of two distinct NSOs; it has at most algebraic complexity ∼M3 deriving from condition 10. Thus, when only conditions 9–11 are imposed, the combinatorial complexity of the ground-state problem resides entirely in the πod of [8]. We extract the sign s(ii′m) of the coefficient Cii′m in [8] and, relating its magnitude to the joint N-OP p11⋯1(ii′m)≡Cii′m2, we rewrite[8] asπod(x1′x2′;x1x2)=∑i<i′≠j<j′,mi,i′,j,j′∉ms(ii′m)s(jj′m)p11⋯11/2(ii′m)p11⋯11/2(jj′m)(ψi(x1′)ψi′(x2′)−ψi′(x1′)ψi(x2′))(ψj(x1)ψj′(x2)−ψj′(x1)ψj(x2)).[12]We suppose that a variational approximation exists in whichs(ii′m)s(jj′m)=s(ii′)s(jj′),∀m.[13]This sign conjecture reduces the sign complexity to algebraic, scaling as M2. πod simplifies toπod(x1′x2′;x1x2)=∑i<i′≠j<j′s(ii′)s(jj′)[∑mi,i′,j,j′∉mp11⋯11/2(ii′m)p11⋯11/2(jj′m)](ψi(x1′)ψi′(x2′)−ψi′(x1′)ψi(x2′))(ψj(x1)ψj′(x2)−ψj′(x1)ψj(x2)).[14]The quantities p11⋯11/2(ii′m) and p11⋯11/2(jj′m) are mth components of vectors p11⋯11/2(ii′) and p11⋯11/2(jj′). The bracketed quantity in [14] is their scalar product. Express it as∑mi,i′,j,j′∉mp11⋯11/2(ii′m)p11⋯11/2(jj′m)=p11001/2(ii′jj′)p00111/2(ii′jj′)ξ(ii′jj′),[15]where p1100(ii′jj′) is the square magnitude of the vector p11⋯11/2(ii′) and p0011(ii′jj′) that of p11⋯11/2(jj′). p1100(ii′jj′) is the probability that ψi and ψi′ are occupied while ψj and ψj′ are not:p1100(ii′jj′)=∑mi,i′,j,j′∉mp11⋯1(ii′m)=∑n Cn2 νi,nνi′,n(1−νj,n)(1−νj′,n),and the reverse is true for p0011(ii′jj′).

The Schwarz inequality 0≤ξ(ii′jj′)≤1 imposes bounds on ξ(ii′jj′), the cosine of the hyperangle between the vectors. The upper bound ξ=1 is exact for N=2. Substituting [15] into [14] yieldsπod(x1′x2′;x1x2)=∑i<i′≠j<j′s(ii′)s(jj′)[p1100(ii′jj′) p0011(ii′jj′)]1/2ξ(ii′jj′)(ψi(x1′)ψi′(x2′)−ψi′(x1′)ψi(x2′))(ψj(x1)ψj′(x2)−ψj′(x1)ψj(x2)),[16]in which only ξ(ii′jj′) retains combinatorial complexity:ξ(ii′jj′)=∑′p11⋯11/2(ii′m)p11⋯11/2(jj′m)(∑′p11⋯1(ii′m) ∑′p11⋯1(jj′m))1/2,where the primed sums are over all m with i,i′,j,j′∉m.

Inserting 4-OPs likep1111(ii′kl)=∑nCn2 νii′,n νkl,n;k<l≠i,i′,j,j′in place of the N-OPs in ξ reduces the complexity of πod to algebraic. The resulting approximation,ξ(ii′jj′)≈∑k<l″p11111/2(ii′kl)p11111/2(jj′kl)(∑k<l″p1111(ii′kl) ∑k<l″p1111(jj′kl))1/2,[17]is not variational, but obeys the 0,1 bounds of the Schwarz inequality. It is exact for N=4, and scales as M4. In [17] the doubly primed sums are over the indices k<l, which must differ from i,i′,j,j′. Bounds on the p1111 that are the generalizations of [9–11] for 3-OPs and 4-OPs can be formulated.

The OP-NOFT Energy Functional.

The trial energy E[Ψ]=〈Ψ|H^|Ψ〉, the expectation value of the Hamiltonian H^, is an explicit functional of the 1- and 2-DM:E[Ψ]=E[ρ,π]=tr{ρh^}+tr{πw^}.Here h^ is the single-particle kinetic-energy operator plus the external potential, and w^ is the 2-electron Coulomb interaction. E[ρ,π] splits into two parts, Ed diagonal and Eod off-diagonal in the SD:E=Ed+EodEd=tr{ρh^}+tr{πdw^}Eod=tr{πodw^}.[18]The Hartree–Fock wave function minimizes Ed; πod introduces electron correlation into Eod. The explicit forms of Ed and Eod follow from [5], [7], and [16]:Ed=∑i p1(i)hii+∑i<j p11(ij)[Jij−Kij],[19]where hii=〈ψi|h^|ψi〉, and Jij=〈ψiψi|w^|ψjψj〉 and Kij=〈ψiψj|w^|ψjψi〉 are the Coulomb and exchange integrals, respectively. The second term on the right-hand side of [19] originates from πd, the form of which is represented exactly in our theory. The second term on the right-hand side of [19] contains only positive contributions and is essential; the integral relation connecting π and ρ depends only on πd and guarantees that the E is self-interaction-free.Eod=∑i<i′≠j<j′s(ii′)s(jj′)p11001/2(ii′jj′)p00111/2(ii′jj′)ξ(ii′jj′)[Kii′jj′−Kii′j′j],[20]where Kii′,jj′=〈ψiψi′|w^|ψjψj′〉. Eqs. 18–20 define the OP-NOFT energy functional within the PDC. Including the complexity of efficient evaluation of the matrix elements, it scales as M5 if the N-representability conditions for the 3- and 4-OPs can be limited to those deriving from the (3,q) and (4,q) positivity conditions with q≤4.

Proof of the Sign Conjecture.

A variational sign approximation must be a statement about the sign s(ii′m) or s(jj′m) of each coefficient appearing in [8]. To prove [13], we must find at least one statement in which the m dependences of s(ii′m) and s(jj′m) cancel. We have found two and present one here and one in SI Appendix, section S1. The former is valid for the general case of matrix elements [Kii′,jj′−Kii′,j′j] of arbitrary sign, the latter only for positive ones.

Assigning each index l in Cn a sign s(l) and taking s(n) as their product to form s(n)=∏l∈ns(l) is a variational approximation. Consequently s(ii′m)=s(i)s(i′)s(m), ands(ii′m) s(jj′m)=s(i)s(i′)s(j)s(j′),[21]so that [13] is proved, with s(ii′)=s(i)s(i′).

This approximation treats the form and phase, 0 or π, of each NSO as independent variables. The choice of signs for each index is not specified in [21]. Most energy minimization schemes start with random initial NSOs; similarly the choice of signs in [21] should be random, half positive and half negative. The number of the initial NSOs should be greater than the anticipated value of M to allow for unequal numbers of positive and negative signs of the active NSOs.

This complete factorization of s(n) and thence of s(ii′) is a restrictive approximation. A variational approximation yielding unfactorized s(ii′) could be more accurate. In SI Appendix, section S1 we have introduced a different variational approximation and rule for the signs which leads to [13] without factorization for positive matrix elements. Random assignment of signs in [13] and the sign rule of SI Appendix, section S1 yield identical results where tested, the significance of which is discussed there.

Theory.

The SDs in [1] can be classified by their seniority, the number A of singly occupied one-particle states they contain. For N even and for a global spin singlet (S=0) state, the N-particle Hilbert space divides into sectors of increasing even seniority starting with A=0, where all SDs contain only doubly occupied states. For molecular systems CI expansions converge rapidly with seniority, and DOCI A=0 calculations describe static correlation rather well, as demonstrated in ref. 20.

The PDC is equivalent to a restriction on seniorities in that seniorities differing only by 4 are allowed. Recent CI calculations for systems with even numbers of electrons showed that the seniority 2 sector largely decouples from the seniority 0 sector, supporting the accuracy of the PDC (20). These considerations also apply to systems with an odd number of electrons, in which seniorities would be odd but still differ only by 4.

We now formulate OP-NOFT explicitly in the A=0 sector to illustrate further how an OP-NOFT functional is constructed and to prepare for numerical implementation; it becomes OP-NOFT-0, in which the PDC is automatically satisfied. Tracing out the spins, [5] becomesρ(r′,r)=2∑k p1(k)ϕk(r′)ϕk(r).[22]k now labels M(>N/2) active doubly occupied NO states, and the following conditions hold:0≤p1(k)≤1 and 2∑k p1(k)=N.[23]In [22] and [23] p1(k) is the occupation number of either of the paired NSOs having the NO ϕk.

In the 2-DM, double occupancy results in a major simplification of the structure of πod. We make the orbital and spin components of the NSO indices explicit. They take the form is, with i now the orbital index and s=± the spin index. The only index pairs that can enter πod in [16] are i+,i− and j+,j−. The only sets of two index pairs that can enter the right-hand side of [17] are i+,i−, k+,k− and j+,j−, k+,k−. The occupation numbers νi+ and νi− are equal, with values 0 or 1, so that all 4-NSO OPs in [17] and [16] are identical to the corresponding spin-independent 2-NO OPs, e.g., p1111(i+,i−,k+,k−)=p11(ik). The signs in [16] depend on a single orbital index, s(i+,i−)=s(i), and the ξ depend on two-orbital indices, ξ(i+,i−,j+,j−)=ξ(ij). With these simplifications, the 2-DM of [7] and [16] becomesπ(r1′r2′;r1r2)=πd(r1′r2′;r1r2)+πod(r1′r2′;r1r2),[24]after tracing out the spins, withπd(r1′r2′;r1r2)=2∑ij p11(ij)(2ϕi(r1′)ϕj(r2′)ϕi(r1)ϕj(r2)−ϕi(r1′)ϕj(r2′)ϕj(r1)ϕi(r2)),[25]πod(r1′r2′;r1r2)=2∑i≠j s(i)s(j)[p10(ij)p01(ij)]1/2ξ(ij)ϕi(r1′)ϕi(r2′)ϕj(r1)ϕj(r2).[26]The sum in [25] includes the term i=j, for which p11(ii)=p1(i), and ξ(ij) in [26] is nowξ(ij)≈∑k(≠i,j) p111/2(ik)p111/2(jk)[∑k(≠i,j) p11(ik) ∑k(≠i,j) p11(jk)]1/2.[27]The one- and two-orbital OPs of OP-NOFT-0 lie within the same bounds as in the general case, [9] and [10], and their sum rules become, respectively, [23] and2∑j(≠i) p11(ij)=(N−2)p1(i).[28]The π of [24] satisfies two important sum rules:∫dr2 π(rr2;r′r2)=(N−1)ρ(r,r′)∫dr1 dr2 π(r1r2;r1r2)w(r12)≥0.The OP-NOFT-0 form for π, [24–27], is exact when N=2 with ξ=1. It is equivalent to DOCI for N=4 if the signs are correct. We show numerically that this is the case for the sign choice of SI Appendix, section S1 for the Paldus H₄ test, as reported in SI Appendix, section S5. For all of the H₄ configurations studied, the OP-NOFT-0 correlation energy coincides with that of DOCI to numerical precision. That the signs are correct for H₂ and H₄ confirms the validity of the sign rule in those cases and suggests a broader utility. When N>4, the ξ-approximation of [27] and the limitation to the (2,2) and (2,3) positivity conditions break the equivalence to DOCI.

The expectation value E=〈Ψ|H^|Ψ〉 becomesE=2∑i p1(i)〈ϕi|h^|ϕi〉+∑ij p11(ij)(2Jij−Kij)+∑i≠j s(i)s(j)p101/2(ij) p011/2(ij) ξ(ij) Kij,[29]where Jij and Kij are, respectively, positive Hartree and exchange integrals defined in terms of the NOs, Jij=〈ϕiϕi|w^|ϕjϕj〉 and Kij=〈ϕiϕj|w^|ϕjϕi〉. With [27] for ξ(ij), E in [29] is a functional of the NOs and the 1- and 2-state OPs. Kollmar introduced a similar J-K functional but simplified the 2-DM (21). The signs are chosen by a sign rule and are not variables. p10 is related to p11 and p1 by p1(i)=p11(ij)+p10(ij) and is eliminated from the functional. Each sum in the denominator of ξ(ij) in [27] is simplified by the sum rule of [28] to, e.g.,∑k(≠i,j)p11(ik)=12(N−2)p1(i)−p11(ij).As stated above, we assume that the (2,2) and (2,3) positivity conditions are sufficient in practice. Under this circumstance, the infimum of E with respect to the NOs and the OPs, subject to the constraints 23 and 9, 10, and 28, yields a variational approximation to the ground-state energy, apart from the ξ-approximation 27, for N>4.

Eq. 29 is a generalization of the NOFT formulations of 1-DM functional theories, which require only 1-state OPs. The extra complexity from 2-state OPs and implicit 4-state OPs is more than compensated by the substantial gain in accuracy it makes possible. The computational cost of calculating E from [29] scales like Hartree–Fock energy-functional minimization with a greater prefactor [M3 vs. (N/2)3] due to fractional occupation of NOs.

Numerical Results for Simple Molecular Systems.

To test OP-NOFT-0, we studied several diatomic molecules, the Paldus H₄ test, and linear chains of H atoms with open boundary conditions. We included all electrons (core and valence) and expanded the NOs in the Gaussian 6–31G** (Pople's notation for a double-zeta split-valence Gaussian basis set), STO-6G, and cc-pVTZ (correlation consistent valence triple-zeta Gaussian basis set with polarization functions) bases. The constrained minimization was performed by damped Car–Parrinello dynamics (22), as detailed in SI Appendix, section S2.

We started the minimization from NOs and OPs obeying the constraints but otherwise random. The signs were taken from the sign rule of SI Appendix, section S1, Table S2. They were kept fixed during optimization.

At convergence, the active subset of NOs had p1≥10−3. The remaining NOs contributed negligibly to the energy. The same active NOs and signs were found for several test cases starting instead from a sufficiently large set of random NOs, half with positive and half with negative signs.^‡ It is significant for the rule of [21] for arbitrary matrix-element signs that for the systems tested, the Brillouin–Wigner perturbation-theory-based rule of SI Appendix, Table S1, and the alternative of random initial assignment of signs to pairs yield the same results for positive matrix elements. The procedure of [21] also yields half positive and half negative signs for the pairs when signs are assigned randomly to the individual NOs with no reference to the matrix-element signs.

We report the dissociation energy curves of the dimers H₂, LiH, and HF in SI Appendix, section S3, Figs. S1–S3. We performed restricted Hartree–Fock, DFT [PBE (24) and/or PBE0 (25)], complete active space self-consistent field (CASSCF), and coupled cluster with single, double, and perturbative triple electron-hole excitations [CCSD(T)] calculations with the same basis. For H₂, our functional depends only on 1-state OPs and the signs s(i); it reduces to the exact expression of Löwdin and Shull (26). When the s(i) are chosen from SI Appendix, Table S2, the OP-NOFT-0 dissociation energy curve coincides with that of CASSCF at all interatomic separations, implying that this sign rule is exact for H₂. Even in a system as simple as H₂, spin-restricted Hartree–Fock and DFT fail badly at dissociation because these single-reference theories cannot recover the Heitler–London form of the wavefunction.

OP-NOFT-0 becomes identical to DOCI for N=4 when the sign choice is correct. The conditions 9, 10, and 28 simplify in this case as discussed in SI Appendix, section S4. Expression 27 for ξ is exact, but the A=0 restriction is not. We performed the H₄ Paldus test using a minimal 1s basis set for OP-NOFT-0, DOCI, and FCI. The DOCI/OP-NOFT-0 equivalence and the accuracy of the OP-NOFT-0 signs are confirmed by the results reported in SI Appendix, section S5. Whereas DOCI only captures about 25–90% of the configuration-dependent correlation energy in the Paldus test, the OP-NOFT-0 dissociation curve of LiH almost coincides with CASSCF in SI Appendix, Fig. S2, which indicates that its 1s electrons are nearly inert so that higher seniorities contribute negligibly to its ground-state energy.

HF, a 10-electron system, provides a complete test of the relation of OP-NOFT-0 to DOCI. The energies obtained with the basis set 6–31G** are shown in SI Appendix, Fig. S3. OP-NOFT-0 and DOCI are above both CASSCF and CCSD(T), although OP-NOFT-0 lies below DOCI because, although size-consistent, the ξ-approximation is nonvariational here. The OP-NOFT-0 signs are correct. The results are discussed further in SI Appendix, section S3.

Linear H chains are relatively simple systems whose energy surfaces present a serious challenge for single-reference methods. Fig. 1 shows symmetric dissociation energy curves of H₈. OP-NOFT-0 provides a consistent description of the energy close to and everywhere above the CASSCF reference. The breakdown of CCSD(T) at large separations is caused by its single-reference character. The deviation of OP-NOFT-0 from CASSCF should be attributed mainly to the restriction to the A=0 sector, a conclusion supported by the seniority-restricted CI calculations of ref. 20. Close comparison with those calculations is not entirely straightforward, as ref. 20 used the slightly smaller 6–31G basis and a fixed, symmetric, or broken symmetry, molecular orbital basis, whereas we used self-consistent NOs.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Symmetric dissociation curve of a linear H₈ chain. The squares indicate one-half of the energy of a H₁₆ chain (black square: Hartree–Fock energy; green squares: OP-NOFT-0 energy). All Hartree–Fock results are from spin-restricted calculations.

The OP-NOFT-0 1-DM displays the entanglement due to correlation through variation of the occupation numbers and the von Neumann entanglement entropy with interatomic separation shown in SI Appendix, section S6, Fig. S5. The increase of entanglement entropy with separation signals a dramatic increase of correlation corresponding to multireference character. The OP-NOFT-0 2-DM gives access to electron pair correlations.

To test the dependence of the accuracy of OP-NOFT-0 on electron number, we studied the symmetric dissociation of H₁₆.^§ Results for the energy of H₁₆ divided by 2 are shown as squares in Fig. 1. OP-NOFT-0 works equally well for this longer chain. The total energy at dissociation is twice that of H₈, and the slightly increased binding energy per atom at equilibrium arises from an increase in the correlation energy, as expected from more effective screening in the larger system.

The N₂ molecule is a severe test because of its triple bond. OP-NOFT-0, DOCI, and FCI results obtained with the minimal basis set STO-6G are compared in Fig. 2. They are similar to those for HF, with OP-NOFT-0 and DOCI both above FCI but OP-NOFT-0 below DOCI. The accuracy of our DOCI results was improved by use of the optimized OP-NOFT-0 NOs for the DOCI basis.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Dissociation curve of the N₂ molecule.

Note that in all of the systems studied, the positivity condition (2,2) was found to be sufficient at near-equilibrium up to intermediate separations dominated by dynamic correlation because the (2,3) condition was automatically satisfied there. Moreover, only in the case of H₈, H₁₆, and N₂ at large separations did inclusion of the (2,3) positivity condition turn out to be essential to prevent runaway from the ground-state solution.