Black-Box Description of Electron Correlation with the Spin-Extended

This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes.

Article pubs.acs.org/JCTC

Black-Box Description of Electron Correlation with the Spin-Extended Configuration Interaction Model: Implementation and Assessment Takashi Tsuchimochi* and Seiichiro Ten-no* Graduate School of System Informatics, Kobe University, Kobe 657-8501, Japan ABSTRACT: In our recent Communication (J. Chem. Phys. 2016, 144, 011101), we proposed Wick’s theorem for nonorthogonal determinants and applied it to spin-extended configuration interaction with singles and doubles (ECISD) based on spin-projected unrestricted Hartree−Fock (SUHF), given that SUHF is a special case of nonorthogonal CI. It was shown that ECISD is accurate for bond-dissociation processes of several representative molecules. In the present work, we give a detailed derivation of ECISD and report an efficient implementation with two physically motivated preconditioning schemes in the generalized Davidson diagonalization for ECISD, whose Hamiltonian and overlap are not diagonal dominant due to SUHF’s nonorthogonal character. In the first approach, we exploit the properties of corresponding pair orbitals and spin-projection operator and rotate the spin-projected CI basis so that the Hamiltonian is approximately diagonal. The second scheme is based on the normal ordered Hamiltonian, which suggests neglecting the expensive two-particle terms in the preconditioning. To enable frozen-core approximations in ECISD, core orbitals were introduced in SUHF. We also show the validity of our method with various numerical examples for static correlations, apart from the left−right correlation in bond-dissociation processes: the ground state energies of the Be isoelectronic series, excitation energies of representative small molecules, and spectroscopic constants of the strongly correlated BN singlet state. Several aspects of ECISD were studied.

1. INTRODUCTION In electronic structure theory, electron correlation is often divided into two conceptual categories: dynamic and static correlations. The dynamic correlation effect accounts for electron collisions that are neglected in the mean-field theory of Hartree−Fock (HF). When HF offers a qualitatively correct description as a zeroth order wave function, it weakly couples with other determinants in the configuration space, and various single-reference approaches such as perturbation theory and coupled cluster (CC) theory can deliver an accurate picture of dynamic correlation to allow quantitative assessments of chemical systems.1,2 In electronically degenerate systems, on the other hand, HF is no longer a good approximation, and it is well-known that most single-reference methods break down. In such cases, multiple determinants play an essential role at zeroth order, giving rise to non-negligible static (nondynamic, strong) correlation due to strong couplings among determinants through a Hamiltonian. In order to achieve “chemical accuracy” in such strongly correlated systems, one is required to be able to describe both dynamic and static correlations correctly. Since it is often difficult to retrieve static correlation after treating dynamic correlation, multireference (MR) approaches, such as complete active space self-consistent field (CASSCF),3−5 are typically employed to obtain static correlation first, followed by perturbation theory,6−8 configuration interaction (CI),9,10 or CC11−14 © XXXX American Chemical Society

to add the residual dynamic correlation. Such MR treatment has been shown to provide very accurate results, but at the same time, it is computationally very demanding compared to its singlereference variants. To reduce the computational effort in capturing static correlation, a variety of approximations have been proposed by many authors, including density matrix renormalization group,15−17 constrained-pairing mean-field theory,18,19 twoelectron reduced density matrix method,20 seniority-based CC approaches,21−23 and projected HF (PHF) theory.24,25 Nonorthogonal CI (NOCI) is another promising route, where one employs a superposition of multiple nonorthogonal Slater determinants to obtain static correlation as well as some dynamic correlation.26−30 The number of determinants required in NOCI is often considerably small compared to that in regular MR approaches such as CASSCF.31,32 This is because nonorthogonal determinants effectively include higher excitation effects due to the exponential ansatz of the Thouless form.33 In this regard, PHF is rephrased as a special type of NOCI, where a set of nonorthogonal determinants are prepared by a well-defined coordinate associated with the unitary rotation operator in a projection operator. Hence, while there are several proposals to Received: February 5, 2016

A

DOI: 10.1021/acs.jctc.6b00137 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation include the residual correlation effects on top of PHF,34−39 it would be most natural to follow the developments in NOCI. However, it is somewhat surprising to us that although NOCI has been known for a long time it has seldom been used as a zeroth order wave function for dynamic correlation. This is perhaps due to the fact that, by increasing the number of nonorthogonal determinants, NOCI will eventually coincide with the exact full CI (FCI) limit. Therefore, it is expected that, with enough determinants, not only static correlation but also a good amount of dynamic correlation can be gained with NOCI by itself, without resorting to a different scheme beyond it. It should be pointed out that, with this expectation, other authors have proposed a NOCI-like approach with a superposition of different PHF states to include dynamic correlation with successful applications.34−36,40 However, since the computational cost for evaluating matrix elements between nonorthogonal determinants at least scales as 6 (N3e ), with Ne being the number of electrons, it will certainly become impossible to carry out a large NOCI calculation with millions of determinants or more. From this point of view, it is therefore deemed more advantageous to develop hierarchical post-NOCI methods, such as CI and CC on top of a NOCI reference in the basis of particle-hole excited configurations from each nonorthogonal Slater determinant. As mentioned, however, we noticed very few works along this line. One of us proposed such a scheme based on two spin-restricted nonorthogonal determinants with very accurate results.41 Yost et al. recently developed a “perturb-then-diagonalize” scheme that uses the first-order wave functions of non-aufbau SCF states as a NOCI basis.42 Another reason as to why the developments in this direction have been hindered may be that such post-NOCI entails efficient evaluation of matrix elements of the Hamiltonian and overlap operators between “excited” nonorthogonal determinants. This may be achieved by invoking the work of Tsuchimochi and Van Voorhis, who utilized perturbative orbital rotations of the principal ground-state determinant.43,44 However, the physical meaning of each obtained term is not immediately clear and, moreover, the algebra would become intractable for double and higher excitations. A more favorable approach to attack this problem is to introduce Wick’s theorem.45 In our previous work (ref 46), henceforth called Paper I, we showed that one can second-quantize Löwdin’s seminal work47,48 to formulate the Wick theorem for a set of nonorthogonal determinants, thereby enabling the systematic development of post-NOCI approaches in a similar fashion as that for single-reference schemes. The theorem was applied to particle-hole configuration interaction based on spin-projected unrestricted HF (SUHF),25 historically called spin-extended HF (EHF),48,49 where one deliberately breaks and restores the Ŝ2 symmetry in a HF determinant through the spin-projection operator. When truncated at double excitations, spin-extended configuration interaction singles and doubles (ECISD) provides a well-defined wave function that is accurate for strongly correlated systems while being an eigenfunction of spin operator Ŝ2.46 In the present article, we continue to develop and further investigate ECISD. Our focus is on a detailed derivation, computationally practical implementation, and numerical and analytical assessments of the method. We extend the previous work and introduce an efficient scheme for Davidson diagonalization50 for ECISD. In Paper I, it was reported that the convergence of ECISD is somewhat slow, typically taking more than 2 to 3 times the cycles that it takes for standard CISD, as the ECISD Hamiltonian is dense. To resolve this issue, we revisit

corresponding pair orbitals,51−54 which ensure a special sparsity in the spin-projected matrix elements and thus help to reduce the number of cycles in the Davdison diagonalization with modified preconditioning. Another way to achieve a faster convergence is to employ a different Hamiltonian partitioning based on the normal-ordered Hamiltonian. We will compare the applicability of each preconditioning scheme. To further ease the computational effort, we discuss core orbitals in SUHF and adopt the frozen-core approximation. ECISD is then applied to various MR systems, as benchmark calculations, apart from bond dissociations, in which the long-range left−right correlation plays a role, namely, singlet−triplet splitting energies, energy diagrams of the ground and excited states of BeH2, correlation energies in Be isoelectronic series, and spectroscopic constants of the singlet BN molecule. We will also assess several properties of the method and discuss the redundancy of single excitations in ECISD, which actually helps the convergence of Davidson diagonalization. We also address the importance of solving the ECISD equation compared to projecting a converged unrestricted CISD (UCISD) wave function. Finally, approximate size-intensivity in excitation energies and the relation between ECISD and multireference configuration interaction (MRCI) are mentioned.

2. THEORY 2.1. Spin-Extended Configuration Interaction. To allow the current work to be as self-contained as possible, we will briefly review our ECISD formulation. We start by writing the spinprojection operator as P̂ =

2S + 1 8π 2

∫Ω dΩ w(Ω) R̂(Ω)

(1)

where S is the spin quantum number, Ω = (α, β, γ) is a set of Euler angles, w(Ω) is Wigner’s D-matrix as weights, and R̂ (Ω) is a rotation operator defined as ̂ ̂ ̂ R(̂ Ω) = e−iαSz e−iβSy e−iγSz

(2)

The spin-projection operator in this form has many advantages over the widely used Löwdin operator48,51 because, while the latter is an intractable many-body operator, the former is represented by an integration of the orbital rotation operator R̂ (Ω) with single-particle operators as its exponent.55 In SUHF, |Φ⟩ is a UHF-like determinant whose orbitals are optimized by minimizing the SUHF energy functional ESUHF =

† ̂ |̂ Φ⟩ ̂ |̂ Φ⟩ ⟨Φ|P ̂ HP ⟨Φ|HP = † ⟨Φ|P|̂ Φ⟩ ⟨Φ|P ̂ P|̂ Φ⟩

(3)

where we have used the facts that P̂ is Hermitian, idempotent, and commutable with Ĥ . Therefore, SUHF falls into the “variation-after-projection” (VAP) category as opposed to “projection-after-variation” (PAV), in which a genuine UHF wave function is spin-projected, which is a widely employed approach in quantum chemistry. There have been many attempts to extend SUHF (as a special case of the more general PHF framework) to include the residual dynamic correlation ΔEc = Eexact − ESUHF.34−40 The most straightforward approach is to introduce a simple CI wave function ansatz ⎛ |Ψ⟩ = P ⎜⎜̂ c0|Φ⟩ + ⎝ B

⎞ cijab|Φijab⟩ + ···⎟⎟ ⎠ ig = ⎜ vo oo −1 vv, † −1 ⎟ (R g ) ⎠ ⎝ R g (R g )

⟨Φia|H̅ |ΨECISD⟩

+ ⟨H̅ ⟩g (>ac>ki + >ai>kc)ckc + (>kc-ai + >ai-kc (16)

ak

+ >ki-ac − >ac-ki + =̅ ic )ckc + ⟨H̅ ⟩g

where o and v stand for occupied and virtual spaces of |Φ⟩, respectively. One special feature of (>g)pq is that it reduces to the delta function when R̂ g = 1̂, in which case one goes back to the orthogonal results. Next, we write the normal-ordered Hamiltonian and properly transform it. Our Hamiltonian is written in a basis as Ĥ =

∑ hpqap†aq + pq

1 4

⎞ ⎛1 ×⎜ >ai>kc>ld + >ac>ki>ld⎟cklcd + (>ai>ld-kc ⎝2 ⎠ 1 + >kc>ld-ai − >ac>ld-ki + >ac>ki-ld 2 1 1 kl ak al + >ki>ld-ac + >ai=̅ cd+>ld=̅ ic + >ki=̅ cd 4 2 ⎫ 1 kl cd ⎪ − >ac=̅ id)ckl ⎬ ⎪ 2 ⎭ (23b)

∑ vrs̅ pqap†aq†asar pqrs

(17)

where hpq and = ⟨pq||rs⟩ are the one-particle Hamiltonian matrix and antisymmetrized two-electron integrals in the Dirac notation, respectively. One could alternatively write Ĥ in b basis, but our results are invariant with this choice. The normal-ordered Hamiltonian Ĥ N,g is then v ̅ pq rs

̂ + VN,g ̂ Ĥ N,g = Ĥ − ⟨Ĥ ⟩g = FN,g

⟨Φijab|H̅ |ΨECISD⟩ = P(ij) P(ab) ∑ wgng g

⎧ ⎪⎛ 1 1 ab⎞ × ⎨⎜ ⟨H̅ ⟩g >ai>bj + >bj-ai + =̅ ij ⎟c0 ⎪⎝ 2 ⎠ 4 ⎩

(18)

where ̂ = FN,g ̂ = VN,g

⎛ ⎞ 1 + ⟨H̅ ⟩g ⎜>bj>ac>ki + >kc>ai>bj⎟ckc ⎝ ⎠ 2

∑ (Fg )pq {ap†aq} (19a)

pq

1 4

+ (>ki>bj-ac − >ac>bj-ki + >ac>ki-bj

∑ vrs̅ pq{ap†aq†asar }

1 1 bk ak >ai>bj-kc + >ac=̅ ij + >bj=̅ ic 2 2 1 1 ab ab + >ki=̅ cj + >kc=̅ ij )ckc + ⟨H̅ ⟩g 2 4 ⎛ 1 ×⎜>ac>bj>ki>ld + >ac>bd>ki>lj ⎝ 4 ⎞ 1 1 + >ai>bj>kc>ld⎟cklcd + ( >ai>kc>ld-bj ⎠ 4 2 + >ai>kc-bj +

(19b)

pqrs

with the definition of normal ordering {Â } ⟨{Â }⟩g ≡ 0

(20)

In eq 19a, we introduced the transition Fock matrix (Fg )pq = hpq +

∑ ⟨pr||qs⟩⟨ar†as⟩g rs

(21)

− >ac>bj>ld-ki + >ki>bj>ld-ac

The excitation operators are also expanded as a combination of normal products and contractions, e.g., for R̂ g|Φijab⟩ = b†a b†b bjbiR̂ g|Φ⟩

+ +

(22)

where P(ij) is the usual permutation operator. The rank in each normal product corresponds to the excitation rank. Thus, one may find the Hamiltonian elements for each normal product using exactly the same rule as that in regular Wick’s theorem. Using the generalized version56 of the nonorthogonal Wick theorem (see eq 62), it becomes readily clear that only the connected contractions have to be evaluated. Using the Wick theorem and the antisymmetry of c, we arrive at

− + +

where we omitted subscript g from >g for brevity. The Einstein convention is assumed for repeated indices except for grid point g, whose summation is explicitly written to emphasize the dependence of all elements except c. Note that in the above equations the summations over k, l and c, d are no longer restricted. Eliminating single excitations from the equations will result in ECID. The Tamm−Dancoff approximation of

⎧ ⟨Φ|H̅ |ΨECISD⟩ = ∑ wgng ⎨⟨H̅ ⟩g c0 + ⟨H̅ ⟩g >kcckc + -kcckc ⎩ g +

⎛ 1 1 kl ⎞ ⎫ ⟨H̅ ⟩g >kc>ldcklcd + ⎜-kc>ld + =̅ cd⎟cklcd ⎬ ⎝ ⎠ ⎭ 2 4

1 >ai>bj>ld-kc + >ai>bd>lj-kc 2 1 1 >bd>ki>lj-ac − >ac>bd>lj-ki 2 2 1 1 1 bk ab ab >ac>ld=̅ ij + >ki>ld=̅ cj + >kc>ld=̅ ij 2 2 8 1 1 kl al bk >ac>bj=̅ id + >ki>bj=̅ cd+>ai>ld=̅ jc 2 2 1 1 kl ak kl >ai>bj=̅ cd + >bd>lj=̅ ic + >ac>bd=̅ ij 8 8 ⎫ ⎪ 1 ab >ki>lj=̅ cd )cklcd ⎬ ⎪ 8 ⎭ (23c)

+ >ac>ld>ki-bj +

⎛ ba†bb†bjbi = {ba†bb†bjbi} + P(ij) P(ab)⎜{ba†bi}(>g)jb ⎝ ⎞ 1 + (>g)ia (>g)jb ⎟ ⎠ 2

⎧ ⎪ = ∑ wgng ⎨(-ai + ⟨H̅ ⟩g >ai)c0 ⎪ g ⎩

(23a) D


Article

Journal of Chemical Theory and Computation

2.3. Davidson Diagonalization. In the regular Davidson diagonalization for an orthogonal CI basis,50 one prepares a set of k approximate eigenvectors c(i) and considers the linear combination ∑ki dic(i) that minimizes the Rayleigh quotient, i.e., the expectation value of the Hamiltonian. The generalization for a nonorthogonal basis is simple: one takes metric S into account, and the Rayleigh quotient becomes

time-dependent PHF is obtained by neglecting the doubles space.43 In deriving eqs 23a−23c, we used the following transformed integrals to simplify the equations (while subscript g was also removed) (-g)rs =

∑ (3 g)rp (Fg )pq (9 g)qs (24a)

pq

(=̅ g)tu vw =

∑ (3 g)tp (3 g)uq vrs̅ pq(9 g)rv (9 g)sw

(k)

ε

(24b)

pqrs

with ⎛ ⟨a† b ⟩ ⟨a† b ⟩ ⎞ ⎛ > oo 0ov ⎞ p∈v q∈o g ⎟ g ⎜ p∈o q∈o g ⎟ ⎜ 3g ≔ ⎜ = ⎜⟨a a† ⟩ ⟨a a† ⟩ ⎟⎟ ⎜−> gvo 1vv ⎟ ⎠ ⎝ p∈v q∈v g ⎠ ⎝ p∈v q∈o g

⎛⟨a† a ⟩ ⟨a b† ⟩ ⎞ ⎛ oo 1 0ov ⎞ p∈o q∈v g ⎟ ⎜ p∈o q∈o g ⎜⎜ vo ⎟ 9g ≔ ⎜ = ⎜⟨a† a ⟩ ⟨a b† ⟩ ⎟⎟ ⎝ > g > gvv ⎟⎠ p ∈ o q ∈ v g p ∈ v q ∈ v g ⎝ ⎠

We should note that

(=̅ g)rspq

=

−(=̅ g)srpq

=

−(=̅ g)qp rs

∑ wgng

(25a)

g

c kc k

kc

(25b)

(j) H͠ ij = c(i) †Hc ̅

(28a)

Sij̃ = c(i) †Sc(j)

(28b)

͠ = ε(k)Sd ̃ Hd

but

(29)

The residual vector k

r(k) = (H̅ − ε(k)S) ∑ di c(i)

(30)

i

accounts for the gradient of eq 27 and gives the direction to lower ε(k). In Davidson diagonalization, the residual vector is preconditioned based on first-order perturbation ξ(k) = (H̅ 0 − ε(k)S0)−1r(k)

(31)

where H̅ 0 and S0 are some approximations to H̅ and S that are easily invertible, e.g., diagonals. Although the nonorthogonality is taken care of with eq 28b, we find it more stable to Gram− Schmidt orthogonalize ξ(k) against previous c(i) under metric S k

c(k + 1) =

∏ (1 − c(i)c(i) †S)ξ(k) (32)

i

followed by the normalization of c so that S̃ is actually the identity matrix. The above procedure is iterated until the Euclidean norm of r(k)†Sr(k) becomes smaller than a preset threshold. In order to maximize the performance of the Davidson diagonalization, one should use robust preconditioning in eq 31. However, since >g is no longer the identity matrix in ECISD, neither H̅ nor S is diagonal dominant. As a result, using diagonals of H̅ and S for preconditioning suffers from a relatively slow convergence in ECISD. In fact, poor convergence behavior with diagonal preconditioning for a projected CI matrix has been numerically reported in ref 43 for time-dependent SUHF. Below, we will address this issue and propose two schemes to increase the efficiency of the method. 2.3.1. CI Basis Rotation. The most natural basis in standard post-HF methods is the canonical orbital basis where the Fock matrix is diagonal. On the other hand, the corresponding pair orbital basis (CO)51,52 has been widely employed for spinprojection methods as well as for nonorthogonal methods due to their useful properties. We will not review the properties of CO in detail here, but a brief summary can be found in Appendix B. Essentially, each occupied α spin orbital ϕi uniquely pairs with its corresponding virtual α orbital ϕi′, the corresponding occupied β orbital ϕi,̅ and the corresponding virtual β orbital ϕi′̅ . Here, a bar and prime denote a β orbital and the corresponding occupied (k+1)

cd ld kl

(26)

⟨Φai |ΨECISD⟩

(27)

Equation 27 can be minimized by solving a small nonorthogonal CI in this reduced space

{c + > c + 12 > > c } 0

k

∑i , j di*djSij̃

where

(=̅ g)rspq ≠ (=̅ g)rspq . The reason for these transformations is as follows. When Wick’s theorem is applied to normal-ordered Hamiltonian elements such as ⟨Φai |Ĥ N,gR̂ g|Φck⟩ = ng⟨a†i aaĤ N,gb†c bk⟩g, it immediately becomes clear that creation operators a†p in F̂N,g and V̂ N,g are always contracted either with aa of the de-excitation operator or bk of the excitation operator. These contractions can be simply expressed as 3 g . In contrast, annihilation operators aq are necessarily contracted with a†i or b†c , represented by 9 g . Therefore, it is convenient to take the summations of indices as in eqs 24a and 24b to transform the Fock and two-electron integrals. This Hamiltonian transformation is similar to but different from the one used in the coupled-cluster theory with exp(T̂ 1),41,57 which does not contain redundant de-excitations, contrary to our R̂ g. The overlap components Sc are readily available from eqs 23a−23c because they are the terms that do not contract (de)excitation operators with H̅ . In other words, one can easily find Sc by collecting the terms that contain ⟨H̅ ⟩g in eqs 23a−23c and then removing ⟨H̅ ⟩g element from the expressions, e.g. ⟨Φ|ΨECISD⟩ =

=

k ∑i , j di*djH͠ ij

⟨Φab ij |ΨECISD⟩

and are similarly obtained. The above ECISD equations bare a computational cost of 6(N 6) with a prefactor of Ngrid due to the last three terms in eq 23c, which are the most intense part; they give rise to computational costs of 6(o3v 3), 6(o 4 v 2), and 6(o2v 4), respectively. Notice that these terms still remain in regular CISD, which is exactly obtained by setting >pq = δpq . Finally, although this section focused on ECISD, these equations are equally applicable to any nonorthogonal-based CI as long as Roo g is nonsingular (which is always true for spinprojection) by realizing P̂ = ∑gwgR̂ g is a generator that, upon its action onto |Φ⟩, produces a superposition of nonorthogonal determinants. However, if P̂ does not commute with H̅ , then one may have to take into account the bra-projection with ⟨Φ|P̂ † instead of ⟨Φ|. Also, it is possible to extend the nonorthogonal Wick theorem for cases where Roo g is singular. The generalization to such cases will be published elsewhere. E


Article

Journal of Chemical Theory and Computation (virtual) counterpart of a virtual (occupied) orbital. Remarkably, in the CO basis of |Φ⟩, >g is sparse and has a well-defined band structure, namely, (>g)pq is zero unless p = q or p and q constitute a corresponding pair (q = p, p,̅ p′, or p′). ̅ Hence, for example, the projection overlap between singly excited determinants ⟨Φiaσσ |P|̂ Φbjττ ⟩ =

∑ wgng(>aσ ,bτ>jτ , iσ + >aσ , iσ>jτ , bτ) g

(33)

where σ and τ are spin indices, is zero for most determinants. The only nonzero elements are ⟨Φia|P|̂ Φia⟩ = ⟨Φai ̅̅ |P|̂ Φai ̅̅ ⟩ ≃ ⟨Φia|P|̂ Φai ̅̅ ⟩

(34)

and ⟨Φii ′|P|̂ Φ jj ′⟩ = ⟨Φ ii ̅̅ ′|P|̂ Φ jj ̅̅ ′⟩ = −⟨Φii ′|P|̂ Φ jj ̅̅ ′⟩ = −⟨Φ ii ̅̅ ′|P|̂ Φ jj ′⟩ (35)

≃0

for i ≠ j, where these equalities can be easily verified by the facts that (>g)pp = (>g) p̅ p̅ and (>g)pp ′ = −(>g)p ′ p. Other overlap couplings are rigorously zero. The elements in eq 35 are a magnitude or two smaller than those of eq 34 and therefore are considered negligible. For this reason, configurations |Φai ⟩ and |Φai ̅ ⟩̅ are strongly interacting through the spin-projection operator. Similarly, each doubly excited configuration interacts mostly with its spin-flip configurations. Figure 1 depicts the structure of S in the CO basis for the linear H4 chain in the 6-31G basis as an example (the matrix dimension is 199). It is noteworthy that, although S contains significant off-diagonal contributions, it is well-structured and somewhat sparse. In the same figure, we also show H̅ in the CO basis. Importantly, both H̅ and S share essentially the same structure. Therefore, one way to define a better preconditioning for the Davidson diagonalization is to rotate the CI basis so that S, and thus H̅ , is approximately diagonal dominant. For the singles space, this task is relatively easy. What eqs 34 and 35 claim is that S in this space is approximately a band matrix, which is supported by the example shown in Figure 1 (singles are 2nd−25th elements). Therefore, one can find a diagonaldominant CI basis for singles by introducing singlet- and tripletlike combinations |1Ξia⟩ =

1 (|Φia⟩ + |Φai ̅̅ ⟩) 2

|2 Ξia⟩ =

1 (|Φia⟩ − |Φai ̅̅ ⟩) 2

Figure 1. Structures of S (top) and H̅ (bottom) of ECISD in the CO basis for the linear H4 chain in the 6-31G basis. The black elements are largely nonzero, the white ones are zero, and those that are in between are gray. The matrix dimension is 199 with 24 single excitations and 174 double excitations, and the CI basis is arranged in the following order: βα ββ |Φ⟩, |Φαα⟩, |Φββ⟩, |Φαα αα⟩, |Φβα⟩, |Φββ⟩.

(36a)

(36b)

which give small off-diagonals, ⟨1Ξai |P̂ |2Ξai ⟩ ≃ 0. For the doubles space, there are three different types of rotation, of which the first two are obtained in a manner similar to that for singles space: (1) α and β excitations to the paired virtual corresponding orbitals (a and a). ̅ 1 ̅ (|Φijaa̅ ̅ ⟩ + |Φaa |1Ξijaa⟩ = i ̅ j ⟩) (37a) 2 |2 Ξijaa⟩ =

1 ̅ (|Φijaa̅ ̅ ⟩ − |Φaa i ̅ j ⟩) 2

|1Ξiiab⟩ =

1 (|Φiabi ̅ ̅ ⟩ + |Φiai̅ b̅ ⟩) 2

(38a)

|2 Ξiiab⟩ =

1 (|Φiabi ̅ ̅ ⟩ − |Φiai̅ b̅ ⟩) 2

(38b)

(3) All possible double excitations from/to different orbitals. For the last category, there are six strongly interacting configurations for each (ijab) spin-subset in the doubles space, ab̅ ab̅ ab̅̅ ab̅ ab̅ i.e., |Φab ij ⟩, |Φi j ̅ ⟩, |Φij̅ ⟩, |Φij̅ ⟩, |Φi j ̅ ⟩, and |Φ i ̅ j ̅⟩. Unfortunately, there are no a priori known combinations for this subspace. Nevertheless, one can form and diagonalize the 6 × 6 overlap matrices for all (ijab) subsets

(37b)

S(ijab)U(ijab) = U(ijab)s(ijab)

(2) α and β excitations from the paired occupied corresponding orbitals (i and i).̅

(39)

where s(ijab) is diagonal, and perform independent rotations of H̅ and S with U(ijab). F


Article

Journal of Chemical Theory and Computation With such a rotated Hamiltonian and overlap, it is expected that the diagonal preconditioning is directly applicable with eq 31. Note that ξ(r) should be appropriately rotated after (before) the preconditioning procedure, as the Hamiltonian contraction is carried out in the original CI basis. The computational effort required for this basis rotation is basically negligible, although this scheme exploits the unique structure of ECISD matrix elements and thus is not applicable to other nonorthogonal approaches such as methods based on spinprojected GHF. 2.3.2. Hamiltonian Partitioning. An alternative, and perhaps more general, preconditioning is to use a different partitioning of the Hamiltonian matrix suitable for perturbation theory H̅ = H̅ 0 + H̅ 1

(40)

From this point of view, we require that H̅ 0 ≫ H̅ 1, which is obviously not the case if one chooses only the diagonal elements as H̅ 0 (see Figure 1). Instead, inspired by the normal-ordered Hamiltonian (eq 18), we partition H̅ as follows (41) F̅ = H̅ − V where F̅ contains all of the terms except those with explicit twopq body terms, i.e., =̅ rs . This amounts to the partitioning between ⟨Ĥ ⟩g + F̂N,g and V̂ N,g in the ECISD Hamiltonian. Also note that this conceptually corresponds to using the orbital energies as the preconditioning for regular CISD, and it is thus assumed that V is small enough compared to F̅. Indeed, Figure 2 shows the structure of F̅, which is very similar to that of H̅ , while V is small. Our preconditioning matrix can be then defined as

A (0k) = F̅ − ε(k)S

and one can determine ξ equations A (0k)ξ(k) = r(k)

(42) (k)

by solving the following linear (43)

A(k) 0

is not diagonal; therefore, eq 43 has to be solved iteratively by using a linear equation solver. The idea is that its contraction with a guess vector is computationally less intense as we avoid 6(N 6) operations inherent in the last three terms of eq 23c. However, solving the above equations still requires an iterative 6(N 5) cost. Hence, even if the convergence for macroiterations (i.e., the Davidson diagonalization) is improved, the microiterations for the linear solver can influence the total computational performance. 2.4. Frozen Core Approximation. In chemical systems, core electrons are often inert and thus do not significantly correlate with valence electrons. In most correlated methods, core orbitals are therefore kept frozen and not involved in calculations to reduce computational costs. These orbitals are typically those with the lowest orbital energies in the HF canonical representation. On the other hand, ECISD is based on SUHF, where there are no well-defined orbital energies. Yet, SUHF has occupied and virtual spaces, and the wave function of SUHF (and thus of ECISD) is invariant with respect to the orbital rotation among each space. Therefore, the simplest way to define core orbitals is to use the semicanonical orbitals of the unprojected state |Φ⟩, i.e., diagonalize the oo block of the UHF Fock matrix constructed from |Φ⟩ and choose the orbitals with the lowest eigenvalues as the core orbitals. While this scheme appears to be intuitive and reasonable, the UHF Fock matrix is symmetry-broken, as are such defined core orbitals, i.e., ϕcore ≠ ϕcore iα iβ . Removing symmetry-broken core

Figure 2. Same as Figure 1 but for F̅ (top) and V (bottom).

orbitals from the ECISD space should not influence the spinprojection, as is clear from eq 4, but we found that the scheme is prone to numerical instabilities. The issue is that an imbalanced treatment between α and β core orbitals makes metric S not exactly singular, that is, the lowest eigenvalues of S are exactly zero without the frozen core approximation, but as soon as different core orbitals are used for different spins, they gain some contamination and turn to a very small but non-negligible value around 10−6 to 10−7. Accordingly, the algorithm for Davidson’s diagonalization suffers from these small eigenvalues and sometime even diverges. More importantly, the desirable structure of the corresponding pair orbitals is lost, and the aforementioned proposal for CI basis rotation (see Section 2.3.1) becomes inapplicable. A different, more promising approach that we employ in the current work is to adopt the natural orbitals (NO) of SUHF, where the SUHF spin-integrated one-particle density matrix D (also called charge density matrix) is diagonal. Then, core orbitals may be identified as the natural orbitals with the largest occupations np and are spin-free. It should be noted that the natural orbitals of a SUHF wave function are equivalent to those of unprojected state, whereas the occupations are different except G


Article

Journal of Chemical Theory and Computation for fully occupied and vacant orbitals, np = 2 and np = 0.58 NOs are directly connected to COs of |Φ⟩ by ϕpNO = ϕpNO = ̅

1 n′p

and diagonalize it in the intrinsic core space. From the viewpoint of MCSCF, the intrinsic core space of SUHF approximately corresponds to the doubly occupied inactive space; hence, the diagonalization of f pq in this block gives effective orbital energies. If the core occupations are less than 2, then performing a spin-constrained optimization in SUHF is suggested. Spincontamination in the underlying UHF wave function |Φ⟩ can be decomposed into the contribution of each natural orbital, which can be eliminated by modifying the Fock matrix during the SCF calculation.59,60 This procedure results in either np = 2 or 0 for the chosen orbitals at convergence. The same trick applies to SUHF.25 For the evaluation of D, we can exploit the spin-free character of the unitary group generator Ê pq = a†paq + a†p̅aq̅

(ϕpCO + ϕpCO) ̅

(44)

where n′p are the occupation numbers of |Φ⟩. If a SUHF occupation np is (very close to) 2, then np′ is as well and therefore ϕNO ≡ ϕCO ≡ ϕCO p p p̅ . The structure of the corresponding pair orbitals is nonetheless retained because the corresponding virtual CO counterparts, ϕCO p′ = ϕp̅′ , belong to the “intrinsic virtual space”, where all orbitals have zero occupations. The natural orbital space is schematically depicted in Figure 3. The problem with

† p

Dpq =

⟨Φ|P ̂ Eq̂ P|̂ Φ⟩ †

⟨Φ|P ̂ P|̂ Φ⟩

p

⟨Φ|Eq̂ P|̂ Φ⟩ = ⟨Φ|P|̂ Φ⟩

(46)

where the numerator can be evaluated with the sum of the α and β transition density matrices in the same orbital basis (such as atomic orbital basis) p

⟨Φ|Eq̂ P|̂ Φ⟩ =

∑ wgng(⟨ap†aq⟩g g

+ ⟨a p†̅ aq ̅ ⟩g )

(47)

Obviously, one can define frozen virtual orbitals in the same manner, if needed, based on the intrinsic virtual block of the generalized Fock matrix f.

3. ILLUSTRATIVE CALCULATIONS 3.1. Convergence Behavior. We first investigate the convergence character of different preconditionings in the Davidson diagonalization. For this purpose, we consider the following test systems: (a) a linear chain of H6 symmetrically separated by 1 Å and the HF molecule separated by (b) 1 and (c) 2 Å. All calculations were done using a 6-31G basis set, and the CI dimensions are about 1000−1300. We used a simple unit vector c0 = 1 as the initial guess for c. In all calculations, the iteration was never restarted to allow for an unbiased comparison. Figure 4 displays the residual (Euclidean) norm ∥r(k)∥ (under metric S) as a function of the number of cycles k. We also show the CISD results with the Davidson diagonalization, which bares a remarkably robust convergence behavior for (b), where no triplet instability was found. However, when static correlation is abundant, as in (a) and (c), symmetry-broken UCISD shows a relatively slower convergence.

Figure 3. Blockwise structure of the natural orbital space. Numbers indicate the occupation number of each orbital.

this scheme, however, is that there can be several (i.e., more than expected) degenerate natural orbitals with ni = 2. These orbitals constitute the “intrinsic core space”, and the true core orbitals are left undefined. This can happen whenever the SUHF correlation is less significant and most orbitals preserve the spin-symmetry. To get around this difficulty, we construct an effective Fock matrix with D f pq = hpq +

⎛

∑ Drs⎝⟨pr|qs⟩ − ⎜

rs

⎞ 1 ⟨pr |sq⟩⎟ ⎠ 2

(45)

Figure 4. Convergence behavior with each preconditioning for (a) H6 (RH−H) = 1 Å, (b) HF (RHF = 1 Å), and (c) HF (RHF = 2 Å). Each plot indicates (△) diagonal preconditioning in the original CI basis, (□) diagonal preconditioning in the rotated CI basis where S is approximately diagonal, (*) diagonal preconditioning in the rotated CI basis where S exactly is diagonal, (⊙) use of F̅ as in eq 42, and (+) UCISD. H


Article


Figure 5. Energy differences E − EFCI for several correction schemes in dissociation curves: (left) HF, (middle) symmetric dissociation of H2O, and (right) N2.

As has been indicated,43,46 the standard diagonal preconditioning in ECISD exhibits very poor convergence character for all cases. It is not always monotonically convergent; the residual norm oscillates, as seen in (b). This behavior is ascribed to the fact that neither H̅ nor S is diagonal dominant; hence, significant off-diagonal contributions are missing at each iteration. Consequently, it typically takes many more cycles than does UCISD. When the CI basis is rotated so that S is approximately diagonal dominant, as described in Section 2.3.1, we obtain a solid improvement in the convergence behavior (denoted “Approximate rotation” in Figure 4). The number of cycles to converge is reduced by a factor of 2−3 compared to that with the original scheme, and it usually converges within 15−20 cycles. Moreover, this scheme is able to converge with fewer cycles than UCISD in (b) and (c), whose off-diagonal Hamiltonian elements seemingly become larger. On the other hand, one would expect that the ideal CI basis is the one where S is exactly diagonal. Although this scheme is practically unfeasible because it requires the diagonalization of prohibitively large S, we also tested such “exact rotation” just for comparison. Interestingly, we found that it improves the results only marginally compared to those with the approximate rotation. Another approach that we proposed is based on the normalordered Fock; it utilizes F̅ in place of H̅ 0 for the preconditioning and solves eq 43 using a linear equation solver. For all cases tested, Fock-based preconditioning achieves the fastest convergence with ECISD, but it does so at the cost of having a computational effort that is more demanding compared to that of the approximate rotation scheme by an order of magnitude (the cost of the latter is only 6(N 4) for each iteration). Therefore, we conclude that, overall, the approximate rotation scheme is the best compromise for ECISD. However, for other nonorthogonal methods, the approximate rotation scheme may not apply, as it is unique to the ECISD Hamiltonian. In such cases, it is expected that Fock-based preconditioning would be more generally useful. 3.2. Comparison of Several Size-Consistent Corrections. In general, as a truncated CI, ECISD is neither sizeconsistent nor size-extensive, as opposed to CC, because simultaneous double excitations (quadruples, sextuples, and so on) are completely neglected, although it can mitigate the size-consistency error resulting from SUHF.46 One proposal to approximately fix the size-inconsistency is to introduce the Davidson correction61,62

ΔE DC = (1 − C02)ΔEc

where C0 is the ratio that the SUHF wave function spans within the Hilbert space of ECISD. Defining a projection operator †

P|̂ Φ⟩⟨Φ|P ̂ Ô = ⟨Φ|P|̂ Φ⟩

(49)

one can extract the information for C0 C0 =

⟨Ψ|Ô |Ψ⟩ = ⟨Ψ|Ψ⟩

|⟨Φ|Ψ⟩| ⟨Φ|P|̂ Φ⟩⟨Ψ|Ψ⟩

(50)

There are also several variants of a posteriori size-consistent corrections besides eq 48, including those due to Brueckner (traditionally called renormalized Davidson correction, ΔERD),63,64 Davidson and Silver65 (ΔEDS), Pople66 (ΔEPC), and Meissner67 (ΔEMC); for the complete expressions of these corrections, see ref 68. The facile computation of C0 with eq 50 allows us to equally adopt these size-consistent corrections to ECISD. Here, we study the performances of these scheme for the dissociation curves of HF, H2O (symmetric dissociation), and N2. Figure 5 illustrates the energy differences from the FCI curves in mHartree with a 6-31G basis. The calculated signless C0 values are also shown as dashed curves. For H2O, ∠HOH is fixed to 109.57°. We have used two core orbitals in the calculations of N2, as described in Section 2.4. As is clearly seen in Figure 5, the intermediate region of the dissociation process is the most difficult to describe, where C0 is found to be minimal. The smaller C0 becomes, the less reliable the SUHF wave function is as a starting point, and vice versa. Generally speaking, the energy correction (in an absolute value) increases in the following order: ΔEMC < ΔEPC < ΔEDC < ΔERD < ΔEDS. In particular, ΔEDS significantly overestimates the quadruple effect when C0 < 0.96. On the other hand, ΔEMC underestimates the effect around equilibrium compared to that at the dissociation limit, but it is not severely affected by C0. ΔEDC and ΔEPC are somewhat stable (although this can be due to some kind of error cancellation). The nonparallelity errors (NPE) from the FCI energy curves are listed in Table 1. It should be stressed that all of the correction schemes improve the ECISD results and outperform CASPT2. In most cases, we find that ΔEDC gives the best results; thus, we will use it for the rest of this work (termed ECISD+Q), unless otherwise noted. 3.3. Singlet−Triplet Splitting Energies. One of the important quantities of organic compounds is singlet−triplet splitting, i.e., the energy difference between singlet and triplet states ΔEST = ES − E T

(48) I

(51) DOI: 10.1021/acs.jctc.6b00137 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article


is introduced. The T-amplitudes in UCCSD for ms = 0 might have been optimized preferably toward the lower triplet state energy. This is interesting because it is widely accepted that UCCSD eliminates the spin-contaminant to first order.72 He and Cremer stated that if a UHF singlet state is contaminated only with a triplet state then UCCSD energy is of a singlet state,73,74 but our results clearly indicate this is not the case. Rather, it appears that a triplet state of ms = 0 is contaminated by a singlet state in UHF; thus, UCCSD gives almost a triplet state energy. This makes much more sense because in the exact FCI limit of UCCn the ms = 0 unrestricted solution must converge to the triplet state since it is the lowest spin state. These observations are also supported by the fact that the perturbative triples correction on UCCSD worsens the results, decreasing ΔEST toward zero. Note that UCCSD energy may be slightly contaminated with higher spin states, but their energies are much higher than those of both singlet and triplet states. In Table 3, we also list the percentage of correlation energy covered by each method (the FCI correlation energy is set to 100%). For both restricted and unrestricted methods, singlet energies are inaccurate either by underestimating or by significantly overestimating the correlation energy. Therefore, it is generally not recommended to use standard single-reference schemes for radicals due to their imbalanced descriptions. As has been reported,69 SUHF gives reasonable ΔEST; the MAE is 6.42 kcal/mol, which is comparable to that with CCSD. It captures static correlation in singlet states, as shown in Table 3, enabling a qualitatively balanced description of ΔEST. With the dynamic correlation effects of ECID and ECISD, surprisingly good results are obtained: MAEs of 0.67 and 0.60 kcal/mol, respectively. It is found that ECID performs as well as ECISD, albeit the former is consistently at the upper bound in energy of the latter. Due to the size-consistency error, both ECID and ECISD correlation energies decrease with the number of electrons (from NH and OH+ to O2 and NF). Nonetheless, the percentage of correlation energy is similar between singlet and triplet states; therefore, their energy differences remain accurate. We note that for triplet states the ECISD correlation energies are similar to those of UCISD (Table 3), which is consistent with the fact that there is little static correlation in these high-spin triplet states. The Davidson correction successfully accounts for the residual correlation responsible for size-consistency. ECISD+Q reproduces the FCI results, except for O2, where the triplet correlation energy is slightly underestimated. The MAE is found to be only 0.30 kcal/mol, showing its potential usefulness for accurate estimation of ΔEST.

Table 1. Nonparallelity Errors for Different Correction Schemes method

HF

H2O

N2a

ECISD ΔEDC ΔERD ΔEDS ΔEPC ΔEMC CASPT2b

0.91 0.13 0.14 0.15 0.19 0.35 1.28

3.77 1.29 1.59 1.96 1.10 1.16 2.60

15.13 1.51 3.65 6.65 3.44 3.28 6.91

a

1s orbitals are frozen. bThe active spaces used are (2e, 2o) for HF, (6e, 5o) for H2O, and (6e, 6o) for N2.

since it gives an estimation of the stability and reactivity of the molecule. When the triplet state is energetically more stable, the molecule exhibits a radical character, and the lowest singlet state is typically open-shell, requiring two configurations. Rivero et al. showed the validity of SUHF and other PHF variants on computing qualitatively correct ΔEST.69 Their results were promising for a mean-field method, but there still remain appreciable errors. Here, we are interested in how the dynamic correlation effects of ECISD can improve the SUHF results. For this purpose, we consider four representative small molecules, NH, OH+, O2, and NF, all at experimental geometries with a 6-31G basis for a comparison with FCI results. For O2 and NF, 1s orbitals were frozen in all of the correlated methods. The calculations were performed with restricted methods (HF, CCSD, CCSD(T)), unrestricted methods (UHF, UCCSD, UCCSD(T), UCISD), and spin-projected methods (SUHF, ECID, ECISD, ECISD+Q). We use ms = 1 for triplet calculations, where a restricted open-shell formalism was employed for restricted methods.59,70,71 The singlet states of unrestricted methods are approximated by ms = 0 solutions. Table 2 shows ΔEST in kcal/mol and errors from FCI in parentheses with various methods. We observe a clear tendency: all spin-restricted methods overestimate ΔEST, whereas spin-unrestricted methods do the opposite. This implies that singlet states are not treated as accurately as triplet states in restricted methods due to missing static correlation. CCSD(T) improves the results of CCSD by almost a factor of 2, but it still yields a non-negligible amount of mean absolute error (MAE), 3.55 kcal/mol. On the other hand, when spins are unrestricted, UHF gives slightly better results (MAE = 17.29 kcal/mol) compared to those with restricted HF (20.03 kcal/mol). However, singlet states are now significantly mixed with a triplet state, and their energies are much lower than they should be. Spincontamination completely ruins ΔEST when dynamic correlation

Table 2. Singlet−Triplet Energy Gap, ΔEST (kcal/mol), for Representative Small Moleculesa molecule

HF

CCSD

CCSD(T)

UHF

UCCSD

UCCSD(T)

UCISD

SUHF

ECID

ECISD

ECISD+Q

FCI

NH

64.44 (18.94) 82.11 (23.77) 42.07 (16.54) 61.72 (20.85) 20.03 20.03

50.85 (5.34) 64.45 (6.11) 31.98 (6.44) 48.42 (7.55) 6.36 6.36

48.09 (2.58) 60.25 (1.91) 30.26 (4.72) 45.85 (4.99) 3.55 3.55

24.59 (−20.92) 30.99 (−27.36) 20.38 (−5.15) 25.13 (−15.73) −17.29 17.29

6.76 (−38.75) 6.78 (−51.56) 10.75 (−14.79) 11.87 (−29.00) −33.52 33.52

0.73 (−44.78) −0.31 (−58.65) 7.80 (−17.73) 6.27 (−34.60) −38.94 38.94

12.71 (−32.80) 13.93 (−44.41) 12.86 (−12.67) 18.20 (−22.66) −28.14 28.14

49.60 (4.09) 62.58 (4.24) 35.99 (10.46) 47.39 (6.52) 6.42 6.42

45.48 (−0.03) 58.39 (0.05) 27.04 (1.28) 39.54 (−1.33) 0.00 0.67

45.65 (0.15) 58.48 (0.14) 26.82 (1.55) 41.45 (0.58) 0.60 0.60

45.49 (−0.01) 58.32 (−0.03) 24.43 (−1.10) 40.81 (−0.06) −0.30 0.30

45.51

OH+ O2b NFb ME MAE a

58.34 25.53 40.87

Numbers in parentheses are the error from FCI. b1s orbitals are frozen. J


Article

Journal of Chemical Theory and Computation Table 3. Percentage of Correlation Energy Gained by Each Method with Respect to FCI (Shown in mEH) NH OH+ O2a NFa a

spin

CCSD

CCSD(T)

UCCSD

UCCSD(T)

UCISD

SUHF

ECID

ECISD

ECISD+Q

Ec (FCI)

ET ES ET ES ET ES ET ES

99.1 91.3 99.2 90.7 96.5 93.3 97.3 92.6

99.7 95.9 99.8 97.1 99.1 96.6 99.4 96.2

99.1 158.3 99.1 173.1 96.8 105.5 97.2 117.1

99.7 167.9 99.8 183.7 99.2 109.2 99.4 122.7

97.3 148.1 97.8 162.1 91.4 99.8 93.3 109.6

8.0 32.3 9.8 38.6 6.8 15.3 3.8 15.6

97.3 98.1 97.5 98.2 91.8 91.7 91.0 91.3

97.7 98.2 98.2 98.7 91.8 92.2 92.6 93.5

100.1 100.1 100.1 100.1 97.7 98.6 98.0 98.4

−70.1 −104.8 −68.4 −111.5 −241.0 −285.0 −197.8 −237.9

1s orbitals are frozen.

Finally, as was mentioned in Section 2.1, the choice of ms in the broken-symmetry HF determinant |Φ⟩ can influence the results of spin-projection for S > 1 . It is therefore interesting to see how 2 much singlet−triplet splitting energies can be affected by employing low-spin triplet states, |S = 1; ms = 0⟩, instead of high-spin states. We found that low-spin results improve ΔEST of SUHF, ECID, and ECISD, at least for the molecules considered in this section; MAEs with low-spin calculations are 3.10, 0.15, and 0.21 kcal/mol, respectively. This improvement can be rationalized from the viewpoint of the number of entangled electrons, Nentangled, which is the number of natural occupations that are significantly deviated from 0 and 2. This quantity is a rough measure of the degree of strong correlation present in a system. Given that Nentangled = 2 for both singlet and low-spin triplet states (in our test cases), it is expected that using a lowspin triple state will ensure a similar static correlation effect to the singlet state and thus provide an improved description of the energy dif ference. Note that high-spin triplet states have zero entanglements and Nentangled = 0; therefore, the correlation effect is considered to be all dynamic. On the other hand, low-spin ECISD+Q calculations slightly worsen the computed ΔEST; its MAE now increases to 0.51 kcal/mol. This indicates that there can be some error cancellation effect in the Davidson correction, as seen in Section 3.2, which now militates against the estimation of ΔEST. Nonetheless, the error is satisfactorily small. 3.4. BeH2. Next, we test our methods on the insertion of Be into H2, a well-known model for static correlation. As described in ref 75, coordinate x characterizes the vertical distance from Be to H2, whose bond length depends on x. The advent of static correlation is schematically understood to be a consequence of near degeneracy between the occupied 2s orbital of Be (a1) and the unoccupied σu orbital of H2 (b2), which is maximized when 2.5 ≤ x ≤ 3.0. Around this region, the potential energy curve of the first excited state crosses that of the ground state, posing a challenge to single-reference methods. Since ECISD offers not only ground state but also excited states as higher roots, we are interested in the performance of ECISD for both states of this strongly correlated system. Here, we report the results using the same basis set as that in ref 76, namely, Be(10s3p/3s2p) and H(4s/2s). Figures 6 and 7 display the potential energy curves of FCI and ECISD and the energy differences of ECISD and CCSD from FCI both as a function of x. The region where near degeneracy is manifest (2.5 ≤ x ≤ 3.0) is filled in red. For the excited state of CCSD, the equation of motion CCSD (EOMCCSD) is available. It is remarkable that ECISD is indistinguishable from FCI for the ground state, with the largest discrepancy being 0.5 mH at x = 4.0. On the other hand, for the first excited state, ECISD underestimates the correlation energy by 4−6 mH, but its potential energy curve is nicely parallel to that of FCI. Single-reference

Figure 6. Potential energy curves of the ground and first excited states. The degenerate region is filled in red.

Figure 7. Energy differences E − E(FCI) for BeH2. C0 coefficient is also shown. The degenerate region is filled in red.

CCSD also predicts a reasonable reaction pathway for the ground state, but there is a noticeable discontinuity at x ≃ 2.75, where the error is largest, 8 mH. This qualitatively incorrect character is attributed to two different reference HF determinants employed for different x values. The first excited state predicted by EOMCCSD also has a kink, but the error is ameliorated compared to that with CCSD, perhaps because the EOM equations flexibly adjust amplitudes in favor of the excited state. In Figure 7, we also show the ECISD+Q curve along with C0. The sudden drop in C0 implies the presence of the same transition as that seen in CCSD, although it is not discontinuous. As a result, the approximate quadruple correction is slightly K


Article

Journal of Chemical Theory and Computation Table 4. Error from the Accurate Correlation Energies (in mEH) for Be Isoelectronic Systems (cc-pCVQZ) Be B+ C2+ N3+ O4+ F5+ Ne6+ MAE NPE a

accuratea

SUHF

ECISD

ECISD+Q

CCSD

CCSD(T)

TPSS

−0.094 34 −0.111 34 −0.126 44 −0.140 53 −0.154 00 −0.167 08 −0.179 88

76.2 87.5 97.5 106.7 115.4 123.7 131.9 105.6 55.7

4.3 4.5 4.6 4.8 5.0 5.1 5.2 4.8 1.0

0.5 0.5 0.5 0.4 0.4 0.3 0.2 0.4 0.3

2.2 2.7 3.0 3.3 3.6 3.8 3.9 3.2 1.8

1.6 2.0 2.2 2.4 2.7 2.9 3.0 2.4 1.4

−4.3 2.0 8.4 15.0 21.6 28.2 34.9 16.3 39.2

Taken from ref 78.

Table 5. Spectroscopic Constants of the Singlet State of BNa method

Ee (a.u.)

Re (Å)

AP1roGc AP1roG-LCCSDc CCSD(T) RMR CCSD(T)d SUHF ECID ECISD ECISD+Q MRCI+Qe

−79.03693 −79.20986 −79.20598

1.310 (+0.012) 1.293 (−0.005) 1.285 (−0.013) 1.298 (±0.000) 1.278 (−0.020) 1.293 (−0.005) 1.291 (−0.007) 1.298 (±0.000) 1.298

AP1roGc AP1roG-LCCSDc CCSD(T) RMR CCSD(T)d SUHF ECID ECISD ECISD+Q MRCI+Qe

−79.07582 −79.28975 −79.26797

Deb (kcal/mol)

ωe (cm−1)

107.0 (−41.6) 169.0 (+20.4) 128.6 (−20.0)

1432 (−218) 1694 (+43) 1698 (+47) 1640 (−11) 1706 (+55) 1685 (+34) 1690 (+39) 1654 (+3) 1651

[cc-pVDZ]

−79.01784 −79.18929 −79.19018 −79.20807 −79.20682

112.7 (−35.9) 142.1 (−6.5) 142.2 (−6.4) 147.6 (−1.0) 148.6 [cc-pVTZ]

−79.03649 −79.24311 −79.24467 −79.26617 −79.26794

1.304 (+0.019) 1.275 (−0.010) 1.274 (−0.011) 1.284 (−0.001) 1.270 (−0.015) 1.278 (−0.007) 1.278 (−0.007) 1.285 (±0.000) 1.285

114.3 (−40.1) 178.6 (+24.2) 126.9 (−27.5) 112.4 (−42.0) 146.0 (−8.4) 146.6 (−7.8) 152.6 (−1.8) 154.4

1630 (−52) 1749 (+67) 1732 (+50) 1681 (+1) 1718 (+36) 1711 (+29) 1705 (+23) 1673 (−9) 1682

Numbers in parentheses are the difference from the MRCI+Q results. bE(50 bohr) − Ee. cTaken from ref 80. dTaken from ref 81. eFull valence active space. Taken from ref 82. a

overemphasized for x > 2.75. Note that excited states are not available with ECISD+Q except with a state-specific approach, which is beyond the scope of the current work. 3.5. Be Series. The Be isoelectronic series showcases another type of static correlation. The occupied 2s and three unoccupied 2p orbitals are close in energy and therefore are nearly degenerate. It is well-known that due to this type of static correlation the total correlation energy linearly increases with atomic number Z,77 as opposed to the He series, where correlation energy becomes constant at large Z. Table 4 tabulates the differences in correlation energy (in mH) from the accurate reference for Z = 4−10, which were computed using experimental data and ab initio calculations.78 The basis set used in this work is Cartesian cc-pCVQZ. The correct description of this series is obtained only by accounting for multiple doubly excited configurations, as is evident from the orbital energy diagram. For this reason, Kohn− Sham density functional theory (TPSS79) is not able to capture the static correlation in the Be series. SUHF cannot accurately predict the trend of correlation energy either because it relies on the pairs of broken-symmetry orbitals and thus takes into account only one doubly excited configuration. We should note that, as a result, the SUHF wave function breaks the spatial symmetry in

this case. The NPE of the SUHF correlation energy is as large as 55.7 mH. CCSD and CCSD(T) correctly describe the linear trend in correlation energy, but the error slightly accumulates as Z increases, resulting in NPEs of 1.8 and 1.4 mH, respectively. This is because static correlation increases with Z and CCSD cannot capture all of it. A similar but less pronounced tendency is found in the ECISD results. The computed energies are more accurate, especially when the size-consistency error is approximately removed in ECISD+Q. It should be noted, however, that they also break the spatial symmetry, as does SUHF, although we speculate that double excitations may ameliorate the symmetry breaking by optimizing the configuration coefficients. 3.6. Spectroscopic Constants of BN. The ground state of BN is a triplet state, but the lowest singlet state has attracted theoretical interest as it represents strong correlation. In the valence-bond picture, the 2pz orbitals of both atoms form a single bond. As a result, the singlet BN dissociates as 1Σ+ → B(2P) + N(2D), where N is “triradical” doublet, as opposed to the triplet dissociation 3Π → B(2P) + N(4S). Hence, the spins are locally frustrated in the N atom in the singlet state, and BN poses a theoretical challenge to single-reference methods. We have carried out calculations with cc-pVDZ and cc-pVTZ, for which the results of the recently developed strong correlation L


Article

Journal of Chemical Theory and Computation theory, AP1roG-LCCSD, are available for comparison. For spinprojected methods, we have used Cartesian basis functions. Table 5 shows the energy at equilibrium Ee, optimized bond length Re, dissociation energy De, and harmonic frequency ωe. The dissociation energies were computed by the supermolecule approach, De = E(50 bohr) − Ee. Since an analytic Hessian is not yet available for projected methods to compute ωe, we have performed numerical fitting to harmonic curves using 20 energy points around Re, each separated by 0.001 Å. The reference data was obtained by MRCI+Q with the full-valence active space (8e, 8o), taken from ref 82. In Table 5, CCSD(T), the gold standard for dynamic correlation, is not as accurate as its multireference version (RMR CCSD(T)), strongly indicating that static correlation plays an important role even at Re. AP1roG, also known as pair CCD, captures static correlation very well, in spite of its singlereference character, by including only the pair double excitations in the cluster amplitude. The method vastly neglects dynamic correlation, similar to SUHF, although it is size-consistent. While SUHF predicts Re to be too short and ωe to be slightly large, AP1roG behaves similarly but opposite. Both methods considerably underestimate Ee, meaning that the correct description of dynamic correlation is necessary at Re. Boguslawski and Ayers recently developed linearized CCSD on top of AP1roG to account for both dynamic and static correlations.80 While AP1roG-LCCSD improves the results over AP1roG, it tends to overestimate the dissociation energy by about 20 kcal/mol. Considering the accurate description of the method at Re, AP1roG-LCCSD appears to miss some correlation at the dissociation limit, where four electrons are experiencing strong correlation. The dynamic correlations gained by ECID and ECISD fix all of the SUHF results. In particular, De is substantially improved; this means that dynamic correlation is much more important at Re than at R → ∞, as expected. Again, ECID results are similar to those of ECISD. ECISD+Q further corrects the errors in ECISD, reproducing exactly the same Re as that with MRCI+Q, an approach that is quite expensive. It is impressive that, generally, the computed values are similar between ECISD+Q and MRCI +Q in Table 5. We see that although ECISD+Q is still not sizeconsistent (neither is MRCI+Q) the resultant size-consistency error is small enough for this system. Finally, it should be pointed out that these results are almost independent of the basis set size: ECISD+Q gives quite similar results to those with MRCI+Q in both cc-pVDZ and cc-pVTZ.

Table 6. Number of Redundant States and Single Excitations H6a HFa H2Oa N2b O2b F2b BNc a

ECID

ECISD

Nsingles

Ndoubles

2 7 7 6 2 10 5

53 59 79 109 119 125 191

54 60 80 110 120 126 192

945 1200 2161 4125 4950 5481 12 528

6-31G. b6-31G. 1s orbitals are frozen. ccc-pVDZ.

doubles space that is not yet clear, we will not investigate it here. What is more important is that this redundancy means single excitations should not contribute much to correlation energy; in fact, ECID correlation energy is very close to that of ECISD in most systems, as has been shown. However, this additional redundancy seems to bring a rather drastic stability to solving ECISD in practice. Occasionally, the convergence of the Davidson algorithm in ECID can become much more unstable than in ECISD. Figure 8 illustrates the convergence behavior of

Figure 8. Convergence behavior of ECID with each preconditioning for H6.

each preconditioning scheme applied to linear H6, which is the same as Figure 4a but with ECID. Clearly, all preconditioning schemes except for exact rotation suffer from a significantly slow convergence compared to that for the ECISD case. Although this does not always happen, including single excitations is likely to introduce a flexibility in the variational parameter space and is considered to be practically beneficial in spin-projected calculations. 4.2. PAV vs VAP. Here, we investigate the consequences of PAV, namely, spin-projected UCISD (PUCISD). There are two different schemes for PAV: in PUCISD(a), the UCISD energy is spin-projected, and in PUCISD(b), c is variationally optimized with eq 6 but with UHF orbitals. For this purpose, we use the F2 molecule, as it requires not only dynamic correlation but also static correlation near equilibrium, representing an ideal test case. As shown in Figure 9, UCISD noticeably underestimates the correlation energy compared to that with FCI. The improvement from spin-projection on top of a UCISD wave function is only marginal (see PUCISD(a)). A large amount of dynamic correlation could be recovered by the Davidson correction, but UCISD+Q produces an inaccurate potential curve and predicts Re to be too short. Hence, the dynamic and static correlation effects would be imbalanced in PUCISD(a). On the other hand, solving c in the correct spin space with eq 6 drastically recovers the correlation energy. Although it is slightly less accurate,

4. DISCUSSION 4.1. On the Importance of Single Excitations. In many cases, ECID energy is quite similar to that of ECISD. The difference between them is often much smaller than that in CC and CI. As in these traditional methods, the role of single excitations in ECISD might be regarded as an approximate orbital rotation on the reference broken-symmetry determinant |Φ⟩.56 However, this should have little effect since the doubles space largely overlaps with the singles one, with the latter being mathematically redundant. In Table 6, we tabulated the number of redundant states in ECID and ECISD as the number of zero eigenvalues in metric S, as well as the number of single and double substitutions, Nsingles and Ndoubles, respectively, for a set of singlet molecules. The redundancy in ECID is small compared to that in ECISD, as expected. Interestingly, we notice that the number of redundant states in ECISD is generically Nsingles − 1. While this implies some rigorous relation between singles and M


Article


Figure 9. Potential energy curves of F2 with a 6-31G basis.

Figure 10. Change in ΔEST with the number of the noninteracting He atoms n.

PUCISD(b) qualitatively reproduces the ECISD curve. This means that c in the UCISD wave function is qualitatively wrong because it is optimized in a biased manner under the influence of the triplet and higher spin-contaminants. It is therefore evident that optimizing c is much more important than optimizing ϕi. In other words, the choice of orbitals may not critically change the results, although using restricted orbitals will result in regular CISD before the Coulson−Fisher point and should be avoided. 4.3. Size-Intensivity in Excitation Energy: A Study of Noninteracting Systems. Because ECISD+Q is size-consistent only approximately, the sum of the energies of two noninteracting fragments A and B, E[A] + E[B], is not precisely equal to the energy of its supermolecule, E[A + B]. It is expected that the size-consistency error ΔSC = E[A + B] − (E[A] + E[B])

(52)

accumulates as the number of fragments becomes larger. As a result, the accuracy of size-extensive observables such as atomization energy is critically affected by this fact and deteriorates with system size. In particular, a projected energy is considered to become that of an unprojected method plus a nonzero constant at the thermodynamic limit.25 On the other hand, excitation energy is a size-intensive quantity and should not depend on the number of fragments. It is therefore plausible that such an “energy difference” is not as severely affected by the size-inconsistency of a method. To numerically validate this, we consider the NH molecule surrounded by n He or Be atoms, none of which is interacting. If the above conjecture is correct, then the singlet−triplet splitting energy ΔEST of a supermolecular system should stay more or less unaltered. We used the same basis set and geometries for NH as those in Section 3.3, and each fragment is separated by 10000 Å. For the He and Be fragments, symmetry-adapted orbitals are used. Figures 10 and 11 illustrate ΔEST computed as a function of n for He and Be, respectively. The FCI ΔEST is 45.51 kcal/mol, and we have already seen in Section 3.3 that the values computed with ECISD and ECISD+Q are in very good agreement with FCI at n = 0 (45.65 and 45.49 kcal/mol). When n is larger, ΔEST of ECISD becomes slightly larger in the cases of both He and Be. The computed ΔEST is more affected by n in the NH···Ben series than in the NH···Hen series because the double excitation effect in the Be atoms is more important than it is in the He atoms. ECISD+Q results show less dependence on n; the largest deviation from FCI results is 0.88 kcal/mol for NH···Be10.

Figure 11. Same as Figure 10 but with Be.

These results can be compared to the same calculations for unrestricted methods, UCISD and UCISD+Q. Initially, both methods predict unphysically small ΔEST at n = 0 (12.71 and 3.14 kcal/mol). Their ΔEST values increase considerably with n as opposed to spin-projected methods; however, they are still significantly underestimated at all n. It should be clear that if a method is size-consistent then ΔEST is constant for different n. This means that UCCSD, for example, could predict the wrong excitation energies regardless of the system size if spin-symmetry is broken; therefore, we suggest that the use of unrestricted methods should be avoided. Why does the increase in n not significantly affect ΔEST in ECISD and ECISD+Q? In Figure 12, we plot the dependence of ΔSC on n for both singlet and triplet states. While there seems to be no clear correlation between ΔSC of singlet and triplet states in unrestricted methods, the errors are almost the same between the two spin states in ECISD and ECISD+Q because of their balanced descriptions. This error cancellation is practically meaningful because ECISD and ECISD+Q can be potentially useful even for much larger systems. However, a general conclusion can be drawn only with more extensive computational analyses. 4.4. Equivalence between the Internally Contracted Scheme and Broken-Symmetry Picture. In Sections 2.1 and 2.2, we used the broken-symmetry picture for particle-hole excitations. In other words, the configurations are prepared by a “excite-then-project” scheme. Alternatively, one could formulate N


Article


Table 7. Total Energies of the HF and N2 Molecules at Re = 0.917 and 1.098 Å (Hartree)a spin-projection SCFb CISD CISD+Qc SCFb CISD CISD+Qc

MR

HF (2e, 2o) −100.008 818 −100.008 818 −100.112 355 −100.112 347 −100.114 834 −100.114 825 N2 (6e, 6o) −108.939 445 −109.015 511 −109.091 159 −109.100 088 −109.101 935 −109.103 122

difference ±0.000 −0.008 −0.009 +76.065 +8.929 +1.186

a c

description as in N2, this equivalence no longer holds; the performance of SUHF is not as good as that of CASSCF because the latter yields the exact electron correlation in the active space. Still, single and double excitations within the broken-symmetry orbitals followed by spin-projection in ECISD can bring a drastically improved description of the “SUHF active space” and partially retrieve the size-consistency.46 In particular, ECISD+Q gives quantitatively accurate results even for a system where a larger active space is required in the general CASSCF. It should be stressed that in SUHF one never specifies an active space and that all electrons (except core electrons if any) are correlated without increasing the computational effort.

Figure 12. Size-consistency error ΔSC in kcal/mol for singlet (S) and triplet (T) states.

ECISD by a “project-then-excite” scheme. Not surprisingly, this exactly corresponds to internally contracted spin-free CI based on a SUHF reference in the NO representation, in which a SUHF wave function can be expressed as MCSCF with spin-restricted determinants. In other words, one may write ⎛ |ΨIC⟩ = ⎜⎜1 + ⎝

p

∑ cq̃ pEq̂ pq

+

⎞

∑ crs̃ pqEr̂ Eŝ ⎟⎟P|̂ Φ⟩ p q

pqrs

⎠

The energy differences are shown in mH. bSUHF and CASSCF. Renormalized Davidson correction, ΔERD.

(53)

̂p where c̃pq and c̃pq rs are the CI coefficients in the NO basis and Eq are,

5. CONCLUDING REMARKS Wick’s theorem for nonorthogonal determinants in the operator basis can drastically simplify the evaluation of matrix elements as in the orthogonal case. In the present work, we applied the theorem to fully derive the working equations of ECISD. It turned out that the Hamiltonian and metric matrices are not diagonal dominant, hindering the straightforward application of diagonal preconditioning in the powerful Davidson technique for large matrix diagonalization. To improve the convergence behavior, we resorted to the corresponding pair orbital basis, which offers a special sparsity in spin-projected matrix elements. It was shown that approximately diagonal-dominant matrices are easily found by appropriate CI basis rotations, enabling a drastic reduction in the number of iterations in the Davidson diagonalization. We also proposed the widely used frozen-core approximation for post-SUHF methods, which enables computational savings as well as a direct comparison against other traditional wave function methods. All of the molecules tested in the current work require a balanced treatment of both dynamic and static correlations. The results obtained with ECISD and its approximately sizeconsistent version, ECISD+Q, are encouraging in many ways. We also studied various properties of ECISD. It was shown that solving the ECISD equation is much more important than simply spin-projecting a converged UCISD wave function. The systematic error cancellation in the size-consistency error between different states also ensures an approximate size-intensivity in excitation energies evaluated with ECISD. Our methods are viewed as an approximation to expensive MRCI in some cases; thus, it can be regarded as an intermediate between singlereference and multireference approaches, in terms of both accuracy and computational effort. In this respect, an extension to a priori size-consistent correction presents no difficulty, and we are currently working along this line. We should also note that the recently proposed multicomponent symmetry-projected

again, the unitary group generators. Since there is no clear partitioning between occupied and virtual orbitals in the NO basis (and in P̂|Φ⟩ itself), p, q, r, and s run over all orbitals. Therefore, it might appear that the NO representation potentially has a larger parameter space than does the brokensymmetry one. However, since Ê pq are spin-free, one may rewrite |ΨIC⟩ in the broken-symmetry representation by commuting Ê pq with P̂ and then transforming NO to an orbital basis that has a clear distinction between the occupied and virtual spaces, e.g., CO. It then becomes apparent that eq 53 contains many redundancies, e.g., de-excitations from the virtual space to the occupied space of |Φ⟩, which leave no net effect when c ̃ is optimized. The only meaningful configurations are those found in the excite-then-project scheme, which was our original proposal. Therefore, we conclude that |ΨIC⟩ gives exactly the same wave function as the broken-symmetry representation. Nonetheless, the equivalence between these two schemes sheds light on a different aspect of ECISD: it can be seen as a compromising alternative to expensive MRCI. Since SUHF mimics CASSCF with a minimal active space in some cases, one can expect ECISD to bring both qualitative and quantitative, but not exceeding, accuracies of the minimal MRCI level. Indeed, we should point out that if the spin-projection in SUHF is restricted to two-electron-two-orbitals (HOMO and LUMO) while other orbitals are symmetry-adapted using spin-constrained optimization then SUHF is identical to CASSCF with (2e, 2o) because they intrinsically have the same parameter space.60 Hence, ECISD on top of SUHF (2e, 2o) in principle gives MRCI with the same active space. This can be seen in Table 7, where we carried out ECISD and MRCI calculations for HF and N2 with active spaces of (2e, 2o) and (6e, 6o), respectively. As expected, the total energies of spin-projection methods and MR methods are numerically identical for (2e, 2o) in HF. On the other hand, when a larger active space is required for a correct zeroth-order O


Article

Journal of Chemical Theory and Computation approach,35,36 an analogue of NOCI, also tackles the residual dynamic correlation problem of PHF and produces accurate results for strongly correlated molecules. Therefore, a comparison of the performances of multicomponent SUHF and our particlehole CI scheme is interesting and will be reported in due time. Finally, many of the proposals in this article are applicable to other NOCI methods. We hope that this work stimulates further developments of post-NOCI.

■

it is shown that eq 56 holds even when virtual orbitals are present. Furthermore, although the number of a and b operators is the same in eq 56, this does not matter because one can expand one basis to another with eqs 13 and 14 to apply the theorem and then transform back to the original basis without loss of generality. We should note that a similar extension was given by Hara and Iwasaki for Bogoliubov-quasiparticles of nonorthogonal Hartree−Fock−Bogoliubov (HFB) states,83 although their results are not directly applicable to the HF case due to the use of Thouless parametrization for the HFB determinants with respect to the bare vacuum |−⟩, which immediately causes a divergence in HF.25 With generalized Löwdin’s results, we can now introduce the nonorthogonal Wick theorem in the operator basis. Recalling that both |Φ⟩ and |Φ̃⟩ are determinants, the concept of normal ordering is apparent only independently for each set of Fermions, ap and bp, as mentioned above. On the other hand, normal ordering of a string that is a mixture of different Fermions should become a linear combination of strings that are ordered differently rather than a simple rearrangement of the creation and annihilation operators in the string. Yet, similar to the extended Wick theorem for the multireference case,84 one can effectively define such normal ordering {ABC···} for two nonorthogonal determinants in the sense that the transition value ⟨{ABC···}⟩R is zero. In other words

APPENDICES

A. Nonorthogonal Wick Theorem

Let us define two different Fermions a†, a and b†, b, which are connected by a unitary rotation as eqs 13 and 14 but with some arbitrary rotation operator R̂ instead of that of spin-projection R̂ g. In standard single-reference approaches, one takes a HF determinant as a starting point Ne

|Φ⟩ =

∏ ai†|−⟩ = |ϕ1 ··· ϕN ⟩ e

(54)

i

where |−⟩ is the physical vacuum. In the regular Wick theorem, normal ordering relative to the Fermi vacuum |Φ⟩ is defined as a product where all quasiparticle annihilation operators a†i , aa lie to the right of all quasiparticle creation operators ai, a†a , which brings a drastic simplification to the evaluation of matrix elements. On the other hand, Fermions b† may constitute the basis for a different determinant |Φ̃⟩

{cp†cq} = cp†cq − ⟨cp†cq⟩R

where c = a, b. Here, the transition density matrix of the last term plays the role of contraction

Ne

|Φ̃⟩ =

∏ bi†|−⟩ = |ϕ1̃ ··· ϕÑ ⟩ = R̂|Φ⟩ e

i

(55)

which is not orthogonal to |Φ⟩. Then, there is no clear definition of normal ordering of a string with b† and b relative to |Φ⟩, and vice versa. In what follows, we show that there is an effective ordering such that a Wick-like theorem similar to the orthogonal case can be formulated. Löwdin showed that an n-body transition matrix (and thus the expectation value of an arbitrary operator) between nonorthogonal determinants can be written as a linear combination of the anti-symmetrized products of one-body transition matrix elements ⟨Φ|ai†1 ai†2⋯b j b j |Φ̃⟩ 2

1

⟨Φ|Φ̃⟩

(58)

What Löwdin’s result (and its generalization to include virtual orbitals) essentially tells us is that the sum of the full contractions is the transition value

where the same sign rule as the single-reference Wick theorem applies. The nonorthogonal Wick theorem then states that any operator can be written as a linear combination of normal products

= ⟨ai†1 ai†2⋯b j b j ⟩R 2

1

⎛⟨a †b ⟩ ⟨a †b ⟩ ⋯⎞ i1 j2 R ⎜ i1 j1 R ⎟ ⎜ ⎟ † † ≡ det ⟨a b ⟩ ⟨a b ⟩ ⎜ i2 j1 R ⎟ i 2 j2 R ⎜ ⎟ ⎝ ⋮ ⋱⎠

where the last term takes into account all single, double, ···, and full contractions. Similar to the generalized Wick theorem in the orthogonal framework,56 one can also write a product of normal-ordered strings as

(56)

where ⟨···⟩R represents the same notation as that in eq 11 but with R̂ . Löwdin’s derivation was based on the first quantization and was therefore restricted to occupied orbitals, as in eq 56. However, the generalization to the virtual space, e.g., excitation operators, is rather straightforward: considering the unitary rotation (eqs 13 and 14) and the anti-commutation relations [ap† , aq]+ = [bp† , bq]+ = δpq

(57a)

[ap† , bq]+

= R *pq

(57b)

[bp† , aq]+ = R qp

(57c)

where the prime on summation means that the contractions occur only between ABC ··· and XYZ ··· but not among themselves. B. Corresponding Pair Orbital

The corresponding spin-pair orbitals (CO) {ϕσp} are usually defined as ones with the following property

⟨ϕiα|ϕjβ ⟩ = σδ i ij P

(63) DOI: 10.1021/acs.jctc.6b00137 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation Extending this feature to virtual orbitals, it becomes clear that there are also one-to-one correspondences between ϕαi and ϕβa as well as between ϕαa and ϕβb , and one finds a unique set of orbitals that satisfy ⟨ϕiα|ϕaβ ⟩ = σi′δia

(64)

⟨ϕaα|ϕi β ⟩ = −σi′δai

(65)

⟨ϕaα|ϕbβ ⟩ = σaδab

(66)

1 − σi2

ACKNOWLEDGMENTS

■

REFERENCES

(1) Bartlett, R. J. Many-Body Perturbation Theory and Coupled Cluster Theory for Electron Correlation in Molecules. Annu. Rev. Phys. Chem. 1981, 32, 359−401. (2) Bartlett, R. J.; Musiał, M. Coupled-cluster theory in quantum chemistry. Rev. Mod. Phys. 2007, 79, 291−352. (3) Roos, B. O.; Taylor, P. R.; Siegbahn, P. E. A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach. Chem. Phys. 1980, 48, 157−173. (4) Siegbahn, P.; Heiberg, A.; Roos, B.; Levy, B. A Comparison of the Super-CI and the Newton-Raphson Scheme in the Complete Active Space SCF Method. Phys. Scr. 1980, 21, 323−327. (5) Siegbahn, P. E. M.; Almlöf, J.; Heiberg, A.; Roos, B. O. The complete active space SCF (CASSCF) method in a Newton-Raphson formulation with application to the HNO molecule. J. Chem. Phys. 1981, 74, 2384−2396. (6) Hirao, K. Multireference Møller-Plesset method. Chem. Phys. Lett. 1992, 190, 374−380. (7) Andersson, K.; Malmqvist, P.-Å; Roos, B. O.; Sadlej, A. J.; Wolinski, K. Second-order perturbation theory with a CASSCF reference function. J. Phys. Chem. 1990, 94, 5483−5488. (8) Angeli, C.; Cimiraglia, R.; Evangelisti, S.; Leininger, T.; Malrieu, J.P. Introduction of n-electron valence states for multireference perturbation theory. J. Chem. Phys. 2001, 114, 10252−10264. (9) Werner, H.-J.; Knowles, P. J. An efficient internally contracted multiconfiguration-reference configuration interaction method. J. Chem. Phys. 1988, 89, 5803−5814. (10) Knowles, P. J.; Werner, H.-J. An efficient method for the evaluation of coupling coefficients in configuration interaction calculations. Chem. Phys. Lett. 1988, 145, 514−522. (11) Mukherjee, D.; Moitra, R. K.; Mukhopadhyay, A. Correlation problem in open-shell atoms and molecules. Mol. Phys. 1975, 30, 1861− 1888. (12) Lindgren, I. A coupled-cluster approach to the many-body perturbation theory for open-shell systems. Int. J. Quantum Chem. 1978, 14, 33−58. (13) Jeziorski, B.; Monkhorst, H. J. Coupled-cluster method for multideterminantal reference states. Phys. Rev. A: At., Mol., Opt. Phys. 1981, 24, 1668−1681. (14) Mukherjee, D.; Pal, S. Use of Cluster Expansion Methods in the Open-Shell Correlation Problem. Adv. Quantum Chem. 1989, 20, 291− 373. (15) White, S. R. Density matrix formulation for quantum renormalization groups. Phys. Rev. Lett. 1992, 69, 2863−2866. (16) Chan, G. K.-L.; Head-Gordon, M. Highly correlated calculations with a polynomial cost algorithm: A study of the density matrix renormalization group. J. Chem. Phys. 2002, 116, 4462−4476. (17) Chan, G. K.-L.; Zgid, D. Chapter 7 The Density Matrix Renormalization Group in Quantum Chemistry. Annu. Rep. Comput. Chem. 2009, 5, 149−162. (18) Tsuchimochi, T.; Scuseria, G. E. Strong correlations via constrained-pairing mean-field theory. J. Chem. Phys. 2009, 131, 121102. (19) Scuseria, G. E.; Tsuchimochi, T. Constrained-pairing mean-field theory. II. Exact treatment of dissociations to nondegenerate orbitals. J. Chem. Phys. 2009, 131, 164119. (20) Gidofalvi, G.; Mazziotti, D. A. Active-space two-electron reduceddensity-matrix method: Complete active-space calculations without diagonalization of the N-electron Hamiltonian. J. Chem. Phys. 2008, 129, 134108. (21) Limacher, P. A.; Ayers, P. W.; Johnson, P. A.; De Baerdemacker, S.; Van Neck, D.; Bultinck, P. A New Mean-Field Method Suitable for Strongly Correlated Electrons: Computationally Facile Antisymmetric

(67)

Note that σa = σi if i and a constitute a corresponding pair. We describe this by a prime on an index, as a′ = i or, equivalently, i′ = a. For sufficiently large basis sets with M being the number of orbitals, there are Nβ occupied and virtual pairs and Nα − Nβ singly occupied α orbitals, which correspond to equivalent singly virtual β orbitals. There are also M − (Nα + Nβ) pure virtual orbitals, which are identical between α and β spaces and thus ⟨ϕαa |ϕβb ⟩ = δab. To better understand the character of this basis, the overlap matrix between α and β spatial orbitals is pictorially shown in Figure 13, where solid lines are added to help identify the occupied−virtual separation.

Figure 13. Orbital overlap matrix in the CO basis.

Working in the CO basis leads to a considerable simplification and reduction in the computational effort with ECISD due to the well-structured sparsity in the Hamiltonian and overlap matrix elements, as shown in Section 2.3.

■

■

The authors are grateful to Drs. Yuhki Ohtsuka, Yu-ya Ohnishi, and Motoyuki Uejima for valuable discussions and helpful comments.

with

σi′ =

Article

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (T.T.). *E-mail: [email protected] (S.T.). Funding

This work was supported in FLAGSHIP2020 by MEXT as priority issue 5 (development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use). We are also indebted to the HPCI System Research project for the use of computer resources (Project ID: hp150228, hp150278). Notes

The authors declare no competing financial interest. Q


Article

Journal of Chemical Theory and Computation Products of Nonorthogonal Geminals. J. Chem. Theory Comput. 2013, 9, 1394−1401. (22) Stein, T.; Henderson, T. M.; Scuseria, G. E. Seniority zero pair coupled cluster doubles theory. J. Chem. Phys. 2014, 140, 214113. (23) Henderson, T. M.; Bulik, I. W.; Stein, T.; Scuseria, G. E. Senioritybased coupled cluster theory. J. Chem. Phys. 2014, 141, 244104. (24) Scuseria, G. E.; Jiménez-Hoyos, C. A.; Henderson, T. M.; Samanta, K.; Ellis, J. K. Projected quasiparticle theory for molecular electronic structure. J. Chem. Phys. 2011, 135, 124108. (25) Jiménez-Hoyos, C. A.; Henderson, T. M.; Tsuchimochi, T.; Scuseria, G. E. Projected Hartree-Fock theory. J. Chem. Phys. 2012, 136, 164109. (26) Fukutome, H. Theory of Resonating Quantum Fluctuations in a Fermion System. Prog. Theor. Phys. 1988, 80, 417−432. (27) Tomita, N.; Ten-no, S.; Tanimura, Y. Ab initio molecular orbital calculations by the resonating Hartree-Fock approach: superposition of non-orthogonal Slater determinants. Chem. Phys. Lett. 1996, 263, 687− 690. (28) Koch, H.; Dalgaard, E. Linear superposition of optimized nonorthogonal Slater determinants for singlet states. Chem. Phys. Lett. 1993, 212, 193−200. (29) Thom, A. J.; Head-Gordon, M. Hartree-Fock solutions as a quasidiabatic basis for nonorthogonal configuration interaction. J. Chem. Phys. 2009, 131, 124113. (30) Burton, H. G. A.; Thom, A. J. W. Holomorphic Hartree-Fock theory: an inherently multireference approach. J. Chem. Theory Comput. 2016, 12, 167−173. (31) Olsen, J. Novel methods for configuration interaction and orbital optimization for wave functions containing non-orthogonal orbitals with applications to the chromium dimer and trimer. J. Chem. Phys. 2015, 143, 114102. (32) McClean, J. R.; Aspuru-Guzik, A. Compact wavefunctions from compressed imaginary time evolution. RSC Adv. 2015, 5, 102277− 102283. (33) Thouless, D. Stability conditions and nuclear rotations in the Hartree-Fock theory. Nucl. Phys. 1960, 21, 225−232. (34) Schmid, K. W.; Zheng, R.-R.; Grümmer, F.; Faessler, A. Beyond symmetry-projected quasi-particle mean fields: A new variational procedure for nuclear structure calculations. Nucl. Phys. A 1989, 499, 63−92. (35) Rodríguez-Guzmán, R.; Jiménez-Hoyos, C. A.; Schutski, R.; Scuseria, G. E. Multireference symmetry-projected variational approaches for ground and excited states of the one-dimensional Hubbard model. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 235129. (36) Jiménez-Hoyos, C. A.; Rodríguez-Guzmán, R.; Scuseria, G. E. Multi-component symmetry-projected approach for molecular ground state correlations. J. Chem. Phys. 2013, 139, 204102. (37) Garza, A. J.; Jiménez-Hoyos, C. A.; Scuseria, G. E. Capturing static and dynamic correlations by a combination of projected Hartree-Fock and density functional theories. J. Chem. Phys. 2013, 138, 134102. (38) Garza, A. J.; Jiménez-Hoyos, C. A.; Scuseria, G. E. Electronic correlation without double counting via a combination of spin projected Hartree-Fock and density functional theories. J. Chem. Phys. 2014, 140, 244102. (39) Tsuchimochi, T.; Van Voorhis, T. Extended Møller-Plesset perturbation theory for dynamical and static correlations. J. Chem. Phys. 2014, 141, 164117. (40) Tomita, N. Many-body wave functions approximated by the superposition of spin-projected nonorthogonal Slater determinants in the resonating Hartree-Fock method. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 045110. (41) Ten-no, S. Superposition of nonorthogonal Slater determinants towards electron correlation problems. Theor. Chem. Acc. 1998, 98, 182−191. (42) Yost, S. R.; Kowalczyk, T.; Van Voorhis, T. A multireference perturbation method using non-orthogonal Hartree-Fock determinants for ground and excited states. J. Chem. Phys. 2013, 139, 174104. (43) Tsuchimochi, T.; Van Voorhis, T. Time-dependent projected Hartree-Fock. J. Chem. Phys. 2015, 142, 124103.

(44) Tsuchimochi, T. Spin-flip configuration interaction singles with exact spin-projection: Theory and applications to strongly correlated systems. J. Chem. Phys. 2015, 143, 144114. (45) Wick, G. C. The Evaluation of the Collision Matrix. Phys. Rev. 1950, 80, 268−272. (46) Tsuchimochi, T.; Ten-no, S. Communication: Configuration interaction combined with spin-projection for strongly correlated molecular electronic structures. J. Chem. Phys. 2016, 144, 011101. Note the equations presented in the supplemental material of Paper I have some errors in factors. (47) Löwdin, P.-O. Quantum Theory of Many-Particle Systems. I. Physical Interpretations by Means of Density Matrices, Natural SpinOrbitals, and Convergence Problems in the Method of Configurational Interaction. Phys. Rev. 1955, 97, 1474−1489. (48) Löwdin, P.-O. Quantum Theory of Many-Particle Systems. III. Extension of the Hartree-Fock Scheme to Include Degenerate Systems and Correlation Effects. Phys. Rev. 1955, 97, 1509−1520. (49) Mayer, I. The Spin-Projected Extended Hartree-Fock Method. Adv. Quantum Chem. 1980, 12, 189−262. (50) Davidson, E. R. The Iterative Calculation of a Few of the Lowest Eigenvalues and Corresponding Eigenvectors of Large Real-Symmetric Matrices. J. Comput. Phys. 1975, 17, 87−94. (51) Löwdin, P.-O. Band Theory, Valence Bond, and Tight-Binding Calculations. J. Appl. Phys. 1962, 33, 251−280. (52) Amos, A. T.; Hall, G. G. Single Determinant Wave Functions. Proc. R. Soc. London, Ser. A 1961, 263, 483−493. (53) Karadakov, P. An extension of the pairing theorem. Int. J. Quantum Chem. 1985, 27, 699−707. (54) Zilberberg, I.; Ruzankin, S. P. Expansion of the unrestricted determinant in the basis of paired orbitals. Chem. Phys. Lett. 2004, 394, 165−170. (55) Ring, P.; Schuck, P. The Nuclear Many-Body Problem; SpringerVerlag: Berlin, 1980. (56) Shavitt, I.; Bartlett, R. J. Many-Body Methods in Chemistry and Physics; Cambridge University Press: Cambridge, 2009. (57) Koch, H.; Christiansen, O.; Kobayashi, R.; Jørgensen, P.; Helgaker, T. A direct atomic orbital driven implementation of the coupled cluster singles and doubles (CCSD) model. Chem. Phys. Lett. 1994, 228, 233−238. (58) Harriman, J. E. Natural Expansion of the First-Order Density Matrix for a Spin-Projected Single Determinant. J. Chem. Phys. 1964, 40, 2827−2839. (59) Tsuchimochi, T.; Scuseria, G. E. ROHF theory made simple. J. Chem. Phys. 2010, 133, 141102. (60) Tsuchimochi, T.; Scuseria, G. E. Constrained active space unrestricted mean-field methods for controlling spin-contamination. J. Chem. Phys. 2011, 134, 064101. (61) Langhoff, S. R.; Davidson, E. R. Configuration interaction calculations on the nitrogen molecule. Int. J. Quantum Chem. 1974, 8, 61−72. (62) Davidson, E. R. The World of Quantum Chemistry; Reidel: Dordrecht, The Netherlands, 1974. (63) Brueckner, K. A. Many-Body Problem for Strongly Interacting Particles. II. Linked Cluster Expansion. Phys. Rev. 1955, 100, 36−45. (64) Luken, W. L. Unlinked cluster corrections for configuration interaction calculations. Chem. Phys. Lett. 1978, 58, 421−424. (65) Davidson, E. R.; Silver, D. W. Size consistency in the dilute helium gas electronic structure. Chem. Phys. Lett. 1977, 52, 403−406. (66) Pople, J. A.; Seeger, R.; Krishnan, R. Variational Configuration Interaction Methods and Comparison with Perturbation Theory. Int. J. Quantum Chem. 1977, 12, 149−163. (67) Meissner, L. Size-consistency corrections for configuration interaction calculations. Chem. Phys. Lett. 1988, 146, 204−210. (68) Duch, W.; Diercksen, G. H. F. Size-extensivity corrections in configuration interaction methods. J. Chem. Phys. 1994, 101, 3018− 3030. (69) Rivero, P.; Jiménez-Hoyos, C. A.; Scuseria, G. E. Predicting Singlet-Triplet Energy Splittings with Projected Hartree-Fock Methods. J. Phys. Chem. A 2013, 117, 8073−8080. R


Article

Journal of Chemical Theory and Computation (70) Rittby, M.; Bartlett, R. J. An open-shell spin-restricted coupled cluster method: application to ionization potentials in nitrogen. J. Phys. Chem. 1988, 92, 3033−3036. (71) Scuseria, G. E. The open-shell restricted Hartree-Fock singles and doubles coupled-cluster method including triple excitations CCSD(T): application to C3+. Chem. Phys. Lett. 1991, 176, 27−35. (72) Schlegel, H. B. Møller-Plesset Perturbation Theory with Spin Projection. J. Phys. Chem. 1988, 92, 3075−3078. (73) He, Y.; Cremer, D. Spin-projected coupled-cluster theory with single and double excitations. Theor. Chem. Acc. 2000, 105, 132−144. (74) Yuan, H.; Cremer, D. The expectation value of the spin operator S2 as a diagnostic tool in couled cluster theory The advantages of using UHF-CCSD theory for the description of homolytic dissociation. Chem. Phys. Lett. 2000, 324, 389−402. (75) Purvis, G. D.; Shepard, R.; Brown, F. B.; Bartlett, R. J. C2V Insertion pathway for BeH2: A test problem for the coupled-cluster single and double excitation model. Int. J. Quantum Chem. 1983, 23, 835−845. (76) Evangelista, F. A.; Allen, W. D.; Schaefer, H. F., III High-order excitations in state-universal and state-specific multireference coupled cluster theories: Model systems. J. Chem. Phys. 2006, 125, 154113. (77) Linderberg, J.; Shull, H. Electronic Correlation Energy in 3- and 4Electron Atoms. J. Mol. Spectrosc. 1961, 5, 1−16. (78) Chakravorty, S. J.; Gwaltney, S. R.; Davidson, E. R.; Parpia, F. A.; Fischer, C. F. Ground-state correlation energies for atomic ions with 3 to 18 electrons. Phys. Rev. A: At., Mol., Opt. Phys. 1993, 47, 3649−3670. (79) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the density functional ladder: Nonempirical meta-generalized gradient approximation designed for molecules and solids. Phys. Rev. Lett. 2003, 91, 146401. (80) Boguslawski, K.; Ayers, P. W. Linearized Coupled Cluster Correction on the Antisymmetric Product of 1 Reference Orbital Geminals. J. Chem. Theory Comput. 2015, 11, 5252−5261. (81) Li, X.; Paldus, J. Singlet-triplet separation in BN and C2: Simple yet exceptional systems for advanced correlated methods. Chem. Phys. Lett. 2006, 431, 179−184. (82) Peterson, K. A. Accurate multireference configuration interaction calculations on the lowest 1 and 3 electronic states of C2, CN, BN, and BO. J. Chem. Phys. 1995, 102, 262−277. (83) Hara, K.; Iwasaki, S. On the quantum number projection: (I). General theory. Nucl. Phys. A 1979, 332, 61−68. (84) Kutzelnigg, W.; Mukherjee, D. Normal order and extended Wick theorem for a multiconfiguration reference wave function. J. Chem. Phys. 1997, 107, 432−449.

S


Black-Box Description of Electron Correlation with the Spin-Extended

Recommend Documents