Perturbation Expansion of Internally Contracted Coupled-Cluster

Subscriber access provided by WEBSTER UNIV

Quantum Electronic Structure

Perturbation expansion of internally contracted coupled-cluster theory up to third order Yuri Alexandre Aoto, Arne Bargholz, Daniel Kats, Hans-Joachim Werner, and Andreas Köhn J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b01301 • Publication Date (Web): 22 Feb 2019 Downloaded from http://pubs.acs.org on March 5, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Perturbation expansion of internally contracted coupled-cluster theory up to third order Yuri Alexandre Aoto,∗,†,‡ Arne Bargholz,∗,† Daniel Kats,∗,†,¶ Hans-Joachim Werner,∗,† and Andreas Köhn∗,† † Institut für Theoretische Chemie, Universität Stuttgart, Pfaffenwaldring 55, D-70569 Stuttgart, Germany ‡ Center for Mathematics Computing and Cognition, Federal University of ABC (UFABC), Avenida dos Estados 5001, Santo André, Brazil ¶ Max Planck Institute for Solid State Research, Heisenbergstraße 1, 70569 Stuttgart, Germany E-mail: [email protected]; [email protected]; [email protected]; [email protected]; [email protected] Abstract The internally contracted multireference coupled-cluster (icMRCC) method is analyzed through third order in perturbation theory. Up to second order, the icMRCC perturbation expansion is equivalent to that of standard Rayleigh-Schrödinger perturbation theory, which is based on a linear ansatz for the wave function, and the resulting theory is, depending on the employed zeroth-order Hamiltonian, equivalent to either second-order complete active space perturbation theory (CASPT2), N -electron valence perturbation theory (NEVPT2), or Fink’s retention of the excitation degree perturbation theory (REPT2). At third order, the icMRCC perturbation expansion features

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

additional terms in comparison to Rayleigh-Schrödinger perturbation theory, but these are shown to be nearly negligibly small by both analytic arguments and numerical examples. Considering these systematic cancellations, however, may be important in future work on approximations to icMRCC theory. In addition, we provide an extensive set of tests of the second and third-order perturbation theories based on three different zeroth-order Hamiltonians, namely the projected effective Fock operator as used for CASPT, the Dyall Hamiltonian as used for NEVPT, and the Fink Hamiltonian used for REPT. While the third-order variant of REPT often gives absolute energies that are rather close to values from higher level calculations, the results for relative energies and spectroscopic constants such as harmonic frequencies, give a less clear picture and a general conclusion about any best zeroth-order Hamiltonian does not emerge from our data. For small active spaces, REPT is rather prone to intruder state problems.

1

Introduction

Many-body perturbation theory (MBPT) is a very powerful approach to studying the electronic structure problem, due to two main reasons. First, because it gives rise to valuable and computationally very efficient procedures for calculating the correlation energy of the electronic system. And second, because it is an useful analysis tool in the development of other approximate electronic structure methods. In single-reference theory, the Møller-Plesset (MP) partitioning of the electronic Hamiltonian, using the Fock operator as zeroth-order Hamiltonian, has become the de facto standard and particularly the second-order theory, MP2, is now widely used as the most simple method to compute the correlation energy of fairly large systems from first principles. The fundamental importance of MBPT for the development and analysis of coupled-cluster (CC) theory is undeniable, 1,2 the most prominent example being the perturbative computation of connected triexcited clusters in the CCSD(T) 3,4 method (coupled-cluster theory with singles, doubles and perturbative triples clusters). Conversely, the exponential ansatz from coupled-cluster theory turned out to be

2


Page 2 of 47

Page 3 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


useful to analyze the structure of the MBPT wave function, 5,6 which originally was formulated by a linear parametrization, requiring additional algebraic effort to reveal factorizations (e.g. that of the second-order MBPT wave function, where the quadruple excitations can be expressed by products of lower-order double excitations 6 ). In multireference perturbation theory (MRPT), required whenever the zeroth-order description of the systems demands more than one configuration, the situation is less clear. One reason is the less obvious choice of the zeroth-order Hamiltonian. The effective Fock operator used in multireference Møller-Plesset (MRMP 7,8 ) theory and complete active space perturbation theory (CASPT 9–11 ) does not have the multiconfigurational reference wave functions as eigenfunctions, and requires additional projection operators to define a proper zeroth-order Hamiltonian. The Dyall Hamiltonian 12 (particularly used in N-electron valence perturbation theory, NEVPT 13,14 ) circumvents this problem, however at the cost of introducing a bielectronic part and neglecting couplings between subspaces. More recently, Fink has proposed the retaining the excitation degree perturbation theory (REPT), 15,16 which found promising application within a matrix product state (MPS) formalism, suggested by Sharma and Alavi as a simple yet accurate alternative to multireference coupled-cluster theory. 17,18 Another source of variety in MRPT approaches is the representation of the perturbed wave function. Although almost exclusively linear parametrizations are chosen, these differ by the degree of flexibility, either uncontracted, 7,8 partially contracted, 11 fully internally contracted, 10,19,20 or even strongly contracted, 13 which also has impact on their qualitative performance, e.g. concerning extensivity. 21 Like for MBPT, the multireference generalization of coupled-cluster theory is not entirely obvious and a multitude of approaches exist. 22–24 Analysis of these methods in terms of MRPT have been undertaken, 25,26 but the literature is more scarce and much less complete as compared to single-reference theory. Some of the present authors have in recent time been involved in the development of the internally contracted multireference coupled-cluster theory (icMRCC), which seems a very promising route to generalizing standard coupled-

3



cluster theory, albeit at the cost of both theoretical and computational complexity. 27–29 A perturbation theory based analysis of the method has previously been presented and was used to derive a multireference analog of CCSD(T). 30 In this work, we set out to extend the perturbation theory based analysis of icMRCC theory, also using different definitions of the zeroth-order Hamiltonian. In Section 2, we develop the perturbative expansion of the icMRCC energy, which we will call the icMRCC perturbation expansion, and establish its connection to the standard Rayleigh-Schrödinger formulation of MRPT based on linear expansions of the perturbed wave function. Non-trivial differences occur at third order already, which gives useful hints for developing approximations in icMRCC theory. In Section 3, we investigate the numerical performance of the second and third-order MRPT expansions based on different choices of the zeroth-order Hamiltonian and compare to icMRCCSD and icMRCCSD(T) results. We focus on both the accuracy of the methods and their stability, in particular with respect to intruder states. The sample calculations include a number of diatomics, ozone, heats of reaction and transition state energies of the F + H2 reaction, and the isomerization of a Cu2 O2 model complex.

2

Theory

We will divide the space of one-electron wave functions into three mutually orthogonal subspaces: closed, active, and virtual space. All reference determinants, |Φµ i, have doubly occupied orbitals in the closed space, no occupation in the virtual space, and varying occupation in the active space and they are assumed to span a complete active space (CAS). All other determinants belong to the complementary space, the orthogonal complement of the reference space. We will adopt the following convention: The indices {i, j, k, l} run only over spin orbitals of the closed space, {u, v, x, w} run only over spin orbitals of the active space, {a, b, c, d} run only over spin orbitals of the virtual space, and {p, q, r, s} run over the entire set of spin orbitals. One and two-electron excitation operators are denoted by

4


Page 4 of 47

Page 5 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a ˆqp = {ˆ a†q a ˆp } and a ˆqs a†q a ˆ†s a ˆr a ˆp }, where {·} indicates normal ordering with respect to the pr = {ˆ common closed-shell part of all reference determinants. The one-electron and the antisympr = (pq|rs) − (ps|rq), respectively, metrized two-electron integrals are denoted by hpq and gqs P P P ip EI = i hii + 21 ij gijij is the inactive energy, and (fI )pq = hpq + i giq are the elements of

the inactive Fock matrix. With these definitions, the usual non-relativistic clamped-nuclei electronic Hamiltonian can be written as ˆ = EI + H

X

ˆqp + (fI )pq a

pq

1 X pr qs ˆ . g a 4 pqrs qs pr

(1)

To perform a perturbative expansion, the Hamiltonian is divided into a zeroth-order Hamiltonian and the perturbation: ˆ =H ˆ (0) + W ˆ . H

(2)

ˆ (0) : The zeroth-order wave function, |Ψ0 i, is required to be an eigenfunction of H ˆ (0) |Ψ0 i = E (0) |Ψ0 i . H

(3)

For a multireference perturbative treatment, an appropriate choice for the zeroth-order wave ˆ projected to the reference function is a CASCI wave function, that is an eigenfunction of H space: ˆ |Ψα i = E α |Ψα i , Pˆ H 0 0 0

(4)

where Pˆ is the projector onto the reference space, Pˆ =

X 0 E D 0 Ψα0 , Ψα0

(5)

α0

and where |Ψα0 i are the normalized CASCI wave functions. The superscripts α or α0 are used to enumerate different CASCI wave functions, if required. If only a single specific reference 5



Page 6 of 47

state is referenced, we will omit the index α for the sake of a compact notation. Note that the zeroth-order energy, E (0) , is not necessarily the reference energy, E0 . These two energies match only for some choices of zeroth-order Hamiltonian, as will be discussed below. The total wave function can be decomposed into the reference wave function and a correlation wave function |Ψi = |Ψ0 i + |Ξi ,

(6)

ˆ . Intermediate normalization, where |Ξi will be expanded in orders of the perturbation W

Ψ0 Ψ0 = Ψ0 Ψ = 1, is assumed throughout this work. The theory outlined in this section employs the “diagonalize-then-perturb” approach that is appropriate for systems without (near-)degeneracies, only. This means that a variation of the wave function within the reference space is not considered and for appropriate choices ˆ (0) , the first-order wave function does not have any component in the reference space. of H Generalization to quasi-degenerate perturbation theory is possible, 31–35 but would add unnecessary complexity to this analysis. The main conclusions of this work will very likely directly transfer to such a generalization.

2.1

Zeroth-order Hamiltonian

The following three zeroth-order Hamiltonians will be considered:

Projected effective Fock operator The effective Fock operator is defined as Fˆ = EI,Fock +

X pq

(fI )pq +

X

X pu v ˆqp . γu a ˆqp = EI,Fock + f˜qp a gqv

(7)

pq

uv

where γuv = hΨ0 |ˆ avu | Ψ0 i is the reduced one-particle density of the CASCI reference state and P EI,Fock = i f˜ii is the energy contribution from the closed orbitals. However, we cannot use Fˆ directly as zeroth-order Hamiltonian, because it is a pure one-electron operator and Equation 3 does not hold in general. To circumvent this problem, the zeroth-order Hamiltonian is 6


Page 7 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


defined using projectors into different excitation spaces, making the CASCI wave function ˆ (0) by construction: 9–11 an eigenfunction of H ˆ (0) = H Fock

XD

E ˆ Fˆ Q ˆ. Ψα0 Fˆ Ψα0 |Ψα0 i hΨα0 | + Q

(8)

α

ˆ is the projector onto the complementary space: Here |Ψα0 i is the α-th CASCI root and Q ˆ = 1 − Pˆ . Q

(9)

D E ˆ (0) , with E (0) = Ψ0 Fˆ Ψ0 . The The CASCI wave function is thus an eigenfunction of H Fock ˆ (0) affects on the size-extensivity of the method, use of projectors in the definition of H Fock discussed in Section 2.4.3.

Dyall Hamiltonian The Dyall Hamiltonian is given by ˆ (0) = EI,Dyall + H Dyall

X

f˜ji a ˆji +

ij

X

f˜ba a ˆba +

X

(fI )uv a ˆvu +

uv

ab

1 X uv wx g a ˆ . 4 uvwx wx uv

(10)

It is still an effective Fock operator in the inactive (closed and virtual) subspaces, but for ˆ (0) does not contain any couplings the active subspace the full Hamiltonian is included. H Dyall between active and inactive spaces and thus equation 3 holds for each CASCI wave function. P ˜i ˆ (0) , hence the value of The closed-shell contribution is EI,Dyall = fi , the same as for H Fock E (0) differs by a constant from E0 . Fink Hamiltonian The partitioning of the Hamiltonian suggested by Fink 15,16 considers of zeroth-order all the parts of the full Hamiltonian that conserve the number of particles in

7



Page 8 of 47

each orbital subspace: ˆ (0) = EI + H Fink

X

ˆji + (fI )ij a

4

ijkl

1X 4

iajb

ˆvu + (fI )uv a

uv

ij

1X

X

X

ˆba + (fI )ab a

ab

1 X uv wx 1 X ab cd ij kl a îj + g a ˆ + gkl g a ˆ + 4 uvwx wx uv 4 abcd cd ab ia jb gjb a îa +

1 X iu jv 1 X ua vb g a ˆ + g a ˆ . 4 iujv jv iu 4 uavb vb ua

(11)

This zeroth-order Hamiltonian satisfies equation 3 as well, as it again contains the full Hamiltonian in the active space but no terms that induce couplings between the reference ˆ (0) indeterminants and determinants from the complementary space. Moreover, since H Fink cludes all two-electrons interactions within each orbital subspace, the use of an effective one-electron operator for the inactive spaces is no longer required. The closed-shell contribution equals the inactive energy EI of the full Hamiltonian, hence E (0) = E0 is fulfilled. ˆ (0) does not change the number of electrons in each subspace, the eight subspaces Since H Fink of the first-order interacting space (FOIS) that have different numbers of electrons in the closed/active/virtual subspaces do not couple through the zeroth-order Hamiltonian. This ˆ (0) , and this property is invoked in the implementations of NEVPT, that is also true for H Dyall treat each of these subspaces separately. 36

2.2

Standard Rayleigh-Schrödinger Perturbation-Theory

In this section we will recall a few well-known equations and properties of the RayleighSchrödinger perturbation theory, for comparison with the formulation derived in the next section. Starting from Equation 6, the complementary space component |Ξi can be parameterized in several ways. Most straightforward is an uncontracted expansion in terms of excited determinants or appropriately spin-adapted configuration state functions, as used for the implementation of MRMP2. 8 Here, we will focus on the internally contracted formulation, 37,38 where |Ξi is parametrized by a linear combination of excitations acting on top of

8


Page 9 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the reference wave function |Ψi = |Ψ0 i + Cˆ |Ψ0 i ,

(12)

with the operator Cˆ given by: Cˆ =

X

tρ τˆρ .

(13)

ρ

The symbol τˆρ denotes a general excitation operator and tρ its corresponding coefficient. The index ρ runs, in principle, over all possible excitations, thus generating a full configuration interaction expansion for |Ψi, and we assume that τˆρ is chosen such that intermediate ˆ (0) is chonormalization is not violated. In fact, if |Ψ0 i is a CASCI wave function and H sen as discussed in section 2.1, internal excitations (excitations to other CASCI states) do not contribute to the first order wave function 10 and we will usually exclude them in our discussion. The amplitudes tρ , and consequently the wave function, are expanded in perturbation orders: (1) tρ = t(0) ρ + tρ + . . .

|Ξi =

(14)

X X X X (i) Ψ(i) = tρ τˆρ |Ψ0 i = Cˆ (i) |Ψ0 i . i

i

ρ

(15)

i

The Schrödinger equation at each perturbation order is given by: ˆ (0) Cˆ (i) |Ψ0 i + W ˆ Cˆ (i−1) |Ψ0 i − H

i X

E (k) Cˆ (i−k) |Ψ0 i = 0 .

(16)

k=0

ˆ (0) , as Note that this equation is consistent at zeroth-order only if |Ψ0 i is eigenfunction of H we have required previously. Using the well known techniques of perturbation theory, the first-order wave function is obtained by solving the equation: D

τˆρ†

E ˆ (0) − E (0) Cˆ (1) + W ˆ H = 0.

(17)

In equation 17, we have introduced a shorthand notation that will be used throughout this 9



Page 10 of 47

article. We will write most expressions as expectation values of the reference wave function D E D E |Ψ0 i: Aˆ = Ψ0 Aˆ Ψ0 . The energy corrections, up to third order, are: D E ˆ (0) H D E ˆ = E0 − E (0) = W D E D E ˆ Cˆ (1) = H ˆ Cˆ (1) = W D E (1) † ˆ (1) ˆ (1) ˆ = (C ) W − E C .

E (0) =

(18)

E (1)

(19)

E (2) (3)

ERS

(20) (21)

We are using the subscript RS for the third-order energy to distinguish it from a slightly different expression for E (3) in the next section.

2.3

Multireference Coupled-Cluster Perturbation-Theory

We start with the following parametrization of the total wave function: ˆ

|Ψi = eT |Ψ0 i .

(22)

ˆ Equation 13, but of course with different amplitudes, The operator Tˆ is formally similar to C, ˆ Equation 22 is the icMRCC ansatz and the energy expression and the as compared to C. amplitude equations of the icMRCC method can be fused into a single energy functional: 29 D

ˆ −Tˆ He ˆ Tˆ L = (1 + Λ)e

E

.

(23)

ˆ = P λρ (ˆ The Lagrange multipliers, λρ , are implicitly given by Λ τρ )† . Differentiation of ρ equation 23 with respect to λρ recovers the amplitude equations of icMRCC theory and the Lagrangian equals the energy expression for optimized amplitudes. In contrast to previous work on icMRCC, we do not consider a relaxation of the reference state |Ψ0 i. In fact, ˆ (0) given in section 2.1, such a change of the reference state under the assumptions for H will not occur for the first-order wave function. 30 A change of |Ψ0 i can be important in 10


Page 11 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


near-degeneracy cases, which we exclude here from our discussion (see above). We expand the amplitudes and the Lagrange multipliers in orders of the perturbation ˆ , and by collecting terms of the same total order, we obtain expressions for the perturbed W wave functions and the energy corrections. 1,5 To derive them, we have to assume that the Tˆ operator does not include any excitations purely within the active space. Such excitations are undesirable, as intermediate normalization cannot be fulfilled as soon as quadratic terms in Tˆ occur. The assumption is not problematic for our analysis, however, as long as the zeroth-order wave function fulfills equation 4, because corrections within the reference space will occur only for the second-order wave function, as already pointed out before. This absence of pure active-active excitations implies that: (a) Tˆ(i) |Ψ0 i belongs to the orthogonal ˆ (i) complement of the reference space; and (b) hΨ0 | Tˆ(i) = 0. Analogous properties hold for Λ and thus, together with equation 3, a number of terms have vanishing contributions in the following. The zeroth-order Lagrangian is then

ˆ (0) )e−Tˆ(0) H ˆ (0) eTˆ(0) L(0) = (1 + Λ

ˆ (0) , Tˆ(0) ], Tˆ(0) ] + . . . ˆ (0) + [H ˆ (0) , Tˆ(0) ] + 1 [[H ˆ (0) ) H = (1 + Λ 2

(24) (25)

ˆ (0) = 0 is the only meaningful choice and it becomes immediately clear that Tˆ(0) = 0 and Λ for the zeroth-order amplitudes and multipliers. Thus, |Ψ0 i is the zeroth-order wave function ˆ (0) is the zeroth-order energy: and the expectation value of H

L

(0)

D

E (0) ˆ = H = E (0) .

At first-order, the Lagrangian is:

ˆ + [H ˆ (0) , Tˆ(1) ] + Λ ˆ (1) H ˆ (0) . L(1) = W

11


(26)


Page 12 of 47

Again, using Equation 3, some terms can be eliminated and the first-order energy expression is obtained:

ˆ = E (1) . L(1) = W

(27)

The second-order Lagrangian becomes:

ˆ , Tˆ(1) ] + Λ ˆ (1) W ˆ +Λ ˆ (1) [H ˆ (0) , Tˆ(1) ] . L(2) = [W

(28)

The first-order amplitudes can be determined from the condition that L(2) should be sta(1)

tionary with respect to λρ : ∂L(2)

0=

(1)

∂λρ

D E ˆ + [H ˆ (0) , Tˆ(1) ] = τˆρ† W .

(29)

Using equation 3, this can be turned into D

τˆρ†

E (0) (0) ˆ (1) ˆ ˆ H −E T +W = 0,

(30)

exactly like in standard Rayleigh-Schrödinger perturbation theory, see Equation 17. For optimized Tˆ(1) , the second-order Lagrangian corresponds to the second-order correction to ˆ , Tˆ(1) ]i = hW ˆ Tˆ(1) i, as in equation 20. the energy, E (2) = h[W The first-order Lagrangian multipliers are obtained from L(2) , as well. Using the station(1)

ary condition with respect to tρ , one arrives at

0=

∂L(2) (1)

∂tρ

=

D

E ˆ , τˆρ ] + Λ ˆ (1) [H ˆ (0) , τˆρ ] [W

=

D

E ˆ + [Λ ˆ (1) , H ˆ (0) ])ˆ (W τρ .

(31)

D E ˆ does not contribute, because τˆρ anniThe last equality uses the fact that the term τˆρ W ˆ (0) , thus hilates the bra reference state, and that the reference function is eigenfunction of H

12


Page 13 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the negative term of the commutator can be recast as: D

E D E D E ˆ (1) τˆρ H ˆ (0) = E (0) Λ ˆ (1) τˆρ = H ˆ (0) Λ ˆ (1) τˆρ . Λ

(32)

ˆ (1) = (Tˆ(1) )† . This equation is just the adjoint of Equation 29 and thus Λ At third-order, the Lagrangian is:

L(3) =

ˆ , Tˆ(1) ], Tˆ(1) ] + ˆ , Tˆ(2) ] + 1 [[W [W 2 ˆ (0) , Tˆ(1) ], Tˆ(1) ] + ˆ (1) [H ˆ (0) , Tˆ(2) ] + Λ ˆ (1) [W ˆ , Tˆ(1) ] + Λ ˆ (1) 1 [[H Λ 2 ˆ (2) W ˆ +Λ ˆ (2) [H ˆ (0) , Tˆ(1) ] . Λ

(33)

ˆ (2) and Tˆ(2) contain the equations 29 and 31, which vanish for optimized Tˆ(1) The terms with Λ ˆ (1) = (Tˆ(1) )† (independently of the values of Λ ˆ (2) and Tˆ(2) ). Thus, the third-order and for Λ correction to the energy becomes: (3)

ECC =

1 ˆ , Tˆ(1) ], Tˆ(1) ] + (Tˆ(1) )† [W ˆ , Tˆ(1) ] + 1 (Tˆ(1) )† [[H ˆ (0) , Tˆ(1) ], Tˆ(1) ] . [[W 2 2

(34)

This result clearly differs from the expression in Rayleigh-Schrödinger perturbation theory, equation 21, and we will analyze this difference further in section 2.4.2.

2.4

Analysis of both approaches

The equations from sections 2.2 and 2.3 are independent of the choice of the zeroth-order ˆ formally includes Hamiltonian (as long as equation 3 holds), and the operator Tˆ (or C) excitations from all ranks. Also, the equations hold for the single reference case, where |Ψ0 i equals the Hartree-Fock determinant, although some terms of equation 34 vanish and it becomes equivalent to equation 21. In this subsection we will clarify what happens to these terms in the multireference case.

13



2.4.1

Page 14 of 47

Contributions from different excitation ranks

Both approaches to the multireference perturbation theory, described in Sections 2.2 and 2.3, lead to the same set of equations to determine the first-order wave function, XD

E E D ˆ ˆ (0) − E (0) )ˆ Ψ0 τˆσ† (H τρ Ψ0 tρ = − Ψ0 τˆσ† W Ψ0 .

(35)

ρ

The set of excited wave functions {ˆ τρ |Ψ0 i} is here chosen in terms of internally contracted configurations. These functions are not orthogonal and contain, in general, linear dependencies, but standard techniques exist to transform them into a proper non-singular basis. 9,20,29 For this discussion, we also admit purely internal excitations, such that the expansion space for equation 35 is still complete up to the full configuration interaction limit. A truncation of this space can be introduced based on the first-order interacting space (FOIS), which is the subspace of the complementary space that interacts with the reference wave function |Ψ0 i through the Hamiltonian. Because the Hamiltonian is a two-body operator, the FOIS can be spanned by only single and double excitations on top of |Ψ0 i. 37 Higher rank excitation operators can be chosen such that the configurations created by them are orthogonal to the FOIS, 29,39,40 and the right hand side of Equation 35 is zero for excitations of rank three or E D †ˆ higher. If |Ψ0 i is a state-specific optimized CASSCF wave function, Ψ0 τˆσ W Ψ0 is also zero for single excitations, due to the generalized Brillouin theorem. However, this fact will not be further exploited, and the discussion is also applicable if the orbitals are obtained by some other procedure, e.g. state-averaged CASSCF. Furthermore, as |Ψ0 i is by construction ˆ (0) , equation 3, and also an eigenstate of the full Hamiltonian projected an eigenstate of H to the CASCI space, equation 4, the right-hand side of equation 35 also vanishes for any internal excitation. Based on the same argument, internal, as well as triple and higher rank excitations do not give any (direct) contribution to the second-order energy, Equation 20. However, the left ˆ (0) discussed in section 2.1, hand side of Equation 35 contains, for all the three choices of H

14


Page 15 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


couplings between excitations of different ranks (doubles and triples, for instance). Thus, triple and higher rank excitations have a non-vanishing contribution to the first-order wave function. The coupling will also slightly influence the single and double amplitudes, resulting in an indirect contribution from triple and higher rank excitations to the second-order energy. This effect, however, is generally expected to be small and higher rank excitations are neglected in any practical implementation. 10,11,13 In CASPT, for instance, Equation 35 is explicitly made block diagonal, by projecting out the coupling between the FOIS and the rest of the complementary space. 10,11 In NEVPT, the perturber functions are also chosen to be only within the FOIS. 13 In the present work we will restrict the first-order wave function to only single and double excitations, as well. Internal and external excitations do not couple ˆ (0) and hence the neglect of internal excitations, as already indicated in sections through H 2.2 and 2.3, is justified. With respect to higher rank excitations, we note that they can directly contribute to the third-order energy, and we will discuss this effect in section 2.4.2. Interestingly, the MRPT formulation based on the icMRCC ansatz , developed in Section 2.3, will give an extra numerical argument for considering only single and double excitations.

2.4.2

Differences between Rayleigh-Schrödinger and coupled-cluster based perturbation theories

The expression for the third-order energy obtained from the multireference coupled-cluster approach, equation 34, has a number of extra terms that do not appear in the expression for (3)

ERS , equation 21, obtained using the linear expansion of the wave function. The difference between the two expressions is explicitly: (3)

(3)

ECC − ERS =

1 ˆ ˆ(1) ˆ(1) 1 ˆ(1) † ˆ (0) ˆ(1) ˆ(1) [[W , T ], T ] + (T ) [[H , T ], T ] + 2 2

(1) † (1) (1) † (1) ˆ (Tˆ ) Tˆ ˆ W − (Tˆ ) Tˆ W

15


(36)


Page 16 of 47

In the single-reference case the above equation yields zero. This is because Tˆ(1) contains only ˆ cannot give a non-vanishing matrix double excitations and thus the first term vanishes as W element between the reference determinant and fourfold excited determinants generated by ˆ (0) is a single-particle operator, which cannot connect (Tˆ(1) )2 , whereas in the second term H the doubly excited determinants generated by (Tˆ(1) )† and the fourfold excited determinants generated by (Tˆ(1) )2 . The last two terms simply cancel each other. In the multireference case, each of the terms above has a non-negligible contribution and the last two are not equal, due to the presence of single excitations and excitations to and from the active space. However, in our calculations we numerically observed that equations 21 and 34 give very similar results, although it is not immediately clear that equation 36 should vanish or be small in the multireference case. To see what happens, some algebraic manipulation is applied to Equation 36 to give: (3)

(3)

ECC − ERS =

E D E 1 D ˆ ˆ (0) ] Tˆ(1) Tˆ(1) − (Tˆ(1) )† Tˆ(1) W ˆ + [H ˆ (0) , Tˆ(1) ] W + [(Tˆ(1) )† , H 2

(1) † (1) ˆ (Tˆ ) Tˆ (37) + W

The first two terms of the equation above are very similar to the equations to determine Tˆ(1) ˆ (1) (Equations 17 and 31). The first term, with the factor 1 , can be written as: and Λ 2 E 1 X ˜ D ˆ (1) † ˆ (0) ˆ tσ W + [(T ) , H ] τˆσ , 2 σ

(38)

with τˆσ and t˜σ being products of excitation operators and amplitudes, respectively. If the full set of excitations (including triples and higher rank excitations) were used to construct Tˆ(1) , the excitations τˆσ in the summation above would be part of the projection manifold in Equation 31 and thus this term would vanish for optimal Tˆ(1) . However, as discussed in section 2.4.1, we restrict Tˆ(1) to the FOIS, and thus this term is not exactly zero, because several of the excitations τˆσ are triple and quadruple excitations, not included in the restricted excitation manifold. 16


Page 17 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The analysis of the second term of equation 37 is more intricate. Using the resolution of the identity, we obtain: D

E D E ˆ W ˆ + [H ˆ (0) , Tˆ(1) ] ˆ + (Tˆ(1) )† Tˆ(1) (Pˆ + Q) = (Tˆ(1) )† Tˆ(1) Pˆ W D E ˆ (0) , Tˆ(1) ] + (Tˆ(1) )† Tˆ(1) Pˆ [H D E ˆ W ˆ + [H ˆ (0) , Tˆ(1) ] (Tˆ(1) )† Tˆ(1) Q . (39)

The first term of the right hand side is equal to the last term of equation 37, canceling it: D

ˆ (Tˆ(1) )† Tˆ(1) Pˆ W

E

=

XD

ED E ˆ(1) † ˆ(1) α (0) α ˆ ˆ Ψ0 (T ) T Ψ0 Ψ0 H − H Ψ0

Dα ED E ˆ , = (Tˆ(1) )† Tˆ(1) W

(40)

where the projector Pˆ has been expanded on the basis of the CASCI wave functions (Equation 5), being the reference wave function in question, |Ψ0 i, one of them (say α = 1). The last equality holds after using Equations 3 and 4: D

Ψα0

E

α

α D E ˆ α (0) (0) ˆ δα,1 ˆ Ψ0 Ψ0 = W H − H Ψ 0 = E0 Ψ 0 Ψ 0 − E

ˆ (0) , Tˆ(1) ] |Ψ0 i is in the complementary The second term of Equation 39 is zero, because [H space and the projector Pˆ annihilates it. Finally, the third term in the right hand side ˆ of Equation 39 has an interpretation similar to the term in Equation 38: hΨ0 | (Tˆ(1) )† Tˆ(1) Q is in the complementary space and the whole term is similar to the first-order amplitude ˆ were completely within the manifold where this equation is equation. If hΨ0 | (Tˆ(1) )† Tˆ(1) Q solved, this term would be exactly zero. It contains, however, terms outside the FOIS, such as hΨ0 | a ˆvi âx ua a zv , that is an effective triple excitation not included in the excitation manifold of only singles and doubles excitations. The perturbative analysis of the icMRCC ansatz carried out in Section 2.3 shows that some non-linear terms in Tˆ(1) appear at third-order, that are not present in the standard 17



Page 18 of 47

perturbative expansion in terms of a linear ansatz for the wave function. Based on the discussion of this section, we see that these terms should cancel systematically. This only does not happen exactly due to the rank restriction made in the excitation manifold, that neglects third and higher rank excitations. We numerically observed that the magnitude of equation 37 is indeed very small. For the calculations we carried out (discussed in Section 3), ˆ (0) and H ˆ (0) , but slightly larger the value of the terms are in between 10 to 50 µEh for H Dyall Fink ˆ (0) , in the order of 102 µEh . See the Supporting Information for detailed graphs, e.g. for H Fock Fig. SI-10. These values are quite small and certainly much smaller than the errors of the MRPT3 calculations. However, the individual terms of Equation 36 can be much larger, and thus it is important that all of them are treated together (all retained or all neglected) in any approximation to the icMRCC ansatz . In the following of the present article, all the results for third-order energies are obtained using the traditional and simplest form, Equation 21.

2.4.3

ˆ (0) and size-extensivity Some observations regarding H Fock (0)

ˆ The projectors used in the definition of H Fock destroy the one-particle character of the zerothorder Hamiltonian, as they act in the many-particle space (their second quantization expressions actually depend on the number of particles). This destroys the size-extensivity of the ˆ (0) and Tˆ do not imply connected CASPT2 method, 11,39,40 because commutators between H Fock ˆ (0) into the amplitude equation 30 gives: quantities. Inserting the definition of H Fock D

E E D D E ˆ (τρ )† Fˆ − Fˆ Tˆ(1) = − Ψ0 τˆσ† H Ψ0

where we explicitly used E (0) =

(41)

E E D E D D ˆ ˆ Fˆ and that Ψ0 τˆσ† W Ψ0 = Ψ0 τˆσ† H Ψ0 . The

left-hand side can be rewritten as a connected part and a remainder D

D E E D E D

E ˆ Tˆ(1) ] + (τρ )† Tˆ(1) Fˆ − Fˆ (τρ )† Fˆ − Fˆ Tˆ(1) = (τρ )† [F,

18


(42)

Page 19 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ˆ and it becomes clear that the remainder would vanish, if |Ψ0 i were an eigenfunction of F. ˆ (0) and H ˆ (0) the remainder cancels out and a strictly size-consistent theory Indeed, for H Dyall Fink ˆ this is not the case and small violations of size-consistency (whenever is obtained. For F, the active space is split among subsystems) are reported in the literature. 11,21,39,40 Disconnected terms can also be found in the third-order energy expression, which, if ˆ (0) , can be rewritten as evaluated for H Fock

E

(3) ˆ (0) CC,H Fock

=

1 ˆ Tˆ(1) ], Tˆ(1) ] + (Tˆ(1) )† [H ˆ − F, ˆ Tˆ(1) ] + 1 (Tˆ(1) )† [[F, ˆ Tˆ(1) ], Tˆ(1) ] [[H, 2 2 D E 1D

E (1) † ˆ (1) ˆ (1) † ˆ (1) ˆ (1) ˆ ˆ ˆ ˆ − (T ) T F− F (T ) T T F − Fˆ − 2 D E + (Tˆ(1) )† Tˆ(1) Pˆ Fˆ Tˆ(1) . (43)

Again, the first term is strictly connected but the remainder is not and does not obviously cancel. One could consider a less strict formulation of a connected MRPT variant based on a truly one-particle zeroth-order Hamiltonian, discarding these disconnected terms. However, ill-defined energy denominators appear, leading to intruder states. This is the limiting case (s → ∞) of the MRPT method based on the driven similarity renormalization group (DSRGMRPT) of Li and Evangelista, 41,42 where s parametrizes a similarity transformation of the original Hamiltonian, which regularizes the energy denominators and leads to a method free ˆ (0) operator as is done of intruder states. 41 In the present article, we will simply use the H Fock in the CASPT method.

3

Numerical Examples

In this section we present results of MRPT calculations at second and third order in energy, carried out with the three zeroth-order Hamiltonians described in section 2.1. The objective is to analyze in particular the third-order methods as possible candidates for approximating

19



the computationally more demanding icMRCCSD method. We investigate the general trend in accuracy one can obtain with these zeroth-order Hamiltonians, and how it depends on the choice of the active space. While for most of the considered cases full valence or even doubleshell active spaces are affordable, we judiciously chose a number of smaller active spaces, including ones which are not optimal for the entire range of probed molecular structures. The rationale behind choosing small and non-optimal active spaces is to investigate the robustness of the methods. A comprehensive collection of graphs and tables can be found in the Supporting Information, and only a selection of representative results is discussed here. For most of the present calculations, the cc-pVDZ basis set 43,44 has been used, but calculations with larger basis sets (cc-pVTZ and cc-pVQZ) also have been carried out for some examples, to check if the conclusions drawn from the results with the smaller basis remain. We compared the results to full configuration interaction (FCI) calculations, internally contracted multireference coupled-cluster theory (icMRCC 27–30 ), and/or multireference configuration interaction, as described in Ref. 45,46, with the relaxed version of the Davidson correction 47 (MRCI+Q), depending on the system size. Core orbitals are always excluded from the correlation treatment. Moreover, they are always kept fixed to the full-valence CASSCF orbitals, when comparison is made among calculations with different reference wave functions. All FCI and MRCI calculations were carried out with Molpro, 48 whereas MRPT and icMRCC calculations were carried out in the General Contraction Code (GeCCo), developed in one of our groups, that allows a straightforward implementation and automated evaluation of the working equations. 29,30 CASSCF orbitals and transformed integrals were imported from a development version of Molpro. The present implementation of the MRPT ˆ (0) is equivalent to CASPT, 9–11 with small deviations to what is obtained equations using H Fock from Molpro calculations due to the different contraction scheme. 11 A recent implementation of CASPT2 in Molpro also offers the full contraction scheme, 49 in this case identical results could be obtained from both GeCCo and the new Molpro code. Our calculations

20


Page 20 of 47

Page 21 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ˆ (0) are equivalent to the partially contracted NEVPT method, which also has been using H Dyall confirmed numerically by comparison with calculations from the original version of Angeli et ˆ (0) is used, the retaining the excitation degree al., 13,14,36 as interfaced to Molpro. 48 When H Fink perturbation theory (REPT) method of Fink is obtained, 15,16 and the corresponding MRPT3 results are very close to those of the MPS-LCC method of Sharma and Alavi, 17,18 where their implementation uses an uncontracted basis represented by matrix product states (MPS).

3.1

Diatomic molecules

We calculated the potential energy curves (PEC) for the ground state of a sequence of 3 − 1 + diatomic molecules, namely, CH (2 Π), NH (3 Σ− ), OH (2 Π), C2 (1 Σ+ g ), N2 ( Σg ), O2 ( Σg ),

CN (2 Σ+ ), CO (1 Σ+ ), and NO (2 Π). We used the cc-pVDZ basis set and the quality of the results is compared against FCI, MRCI+Q, and icMRCC calculations. For most of the cases, FCI PECs are taken as reference for the accuracy of the results, except for the O2 and NO molecules, for which the reference PECs are calculated with the icMRCCSD(T) method and full valence active space. We also carried out calculations with the cc-pVTZ and cc-pVQZ basis sets, for the full valence active spaces. FCI calculations are not affordable for these basis sets, but the MRPT results can be compared to icMRCCSD(T), and the trend observed for the cc-pVDZ basis set remains. We will first analyze the general behaviour of the potential energy curves for larger active spaces (full valence, for instance). For an explicit example, see Figure 1 showing the results of full valence active space MRPT calculations on the CN molecule in the X 2 Σ+ electronic ground state. Similar figures are given in the Supporting Information for all diatomic molecules studied here, and the same behavior as for this particular case is found for all calculations with full valence active space. In all these cases, the most accurate MRPT results (in terms of the root-mean-square deviation from the FCI reference curve) are obˆ (0) . At second order, the method usually overshoots the external correlation tained with H Fink energy, defined as the difference between FCI or icMRCCSD(T) and the CASSCF energies. 21



25

(0) ) MRPT2 (HDyall

[2.29 mEh]

E EFCI (mEh)

20 (0) ) MRPT2 (HFock [0.94 mEh]

15 10 5 0

(0) ) MRPT3 (HDyall [0.37 mEh]

(0) ) MRPT3 (HFock

[0.77 mEh]

(0) ) MRPT3 (HFink

[1.11 mEh]

5 92.40 1.8 Energy (Eh)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 47

(0) ) MRPT2 (HFink [1.33 mEh] 2.6 2.8 3.0

2.0

2.2

2.4 R (a0)

2.0

2.2

2.4 2.6 R (a0)

92.45 92.50

1.8

2.8

3.0

Figure 1: Potential energy curves of the X 2 Σ+ state of CN obtained from the MRPT calculations and full valence active space, with the cc-pVDZ basis set. Upper graph shows the difference to the FCI potential energy curve and the lower graph shows total energies. Non-parallelity errors, NPE, calculated as NPE = max(E − EFCI ) − min(E − EFCI ) for the range shown in the graph are reported inside brackets.

22


Page 23 of 47

The deviation to the reference PEC is largely decreased at third order, and the MRPT3 cal(0)

ˆ culations with H Fink using the full valence active spaces present very good results, as already (0)

(0)

ˆ ˆ pointed out by Sharma and Alavi. 17,18 Potential energy curves with H Fock and HDyall are less accurate and tend to underestimate the correlation energy. For the molecules considered (0)

ˆ here, PECs calculated with H Dyall show the largest deviation to the reference PEC, although very often with rather small non-parallelity error (NPE).

75.3 75.4

(0) , CAS (4e, 4o) HDyall

75.5 Energy (Eh)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(0) , HDyall CAS (2e, 2o)

75.6

FCI

75.7 75.8 (0) , CAS (4e, 4o) HFink (0) , CAS (2e, 2o) HFink 75.9 2.0 2.5 3.0 3.5 R (a0)

4.0

ˆ (0) for potential energy curves of C2 Figure 2: Intruder state problems of MRPT using H Fink (state X 1 Σ+ g ) when using small active spaces. Shown are results at the MRPT2 level with ˆ (0) and H ˆ (0) for two active spaces (see text). H Fink Dyall When smaller active spaces are used, the trend of the results from the three zeroth-order Hamiltonians is much less evident. For some cases, the order of performance observed for ˆ (0) the larger active spaces remains, but for others the good performance obtained with H Fink 23



is lost. For example, we observed strong intruder state problems for the C2 molecule for the active spaces CAS (4e, 4o) and CAS (2e, 2o). The active space CAS (2e, 2o) is chosen to describe the two main configurations (2σu )2 (1πu )4 and (3σg )2 (1πu )4 close to equilibrium, while CAS (4e, 4o) includes all four π orbitals and can describe the homolytic cleavage of ˆ (0) close to the molecule. In Figure 2 we see that CAS (2e, 2o) gives good results for H Fink the equilibrium, but severely breaks down beyond 3 a0 , where other configurations start to become important and a number of intruder states appear. Likewise, CAS (4e, 4o) works well for large separations but shows an intruder state at the equilibrium distance, where the (3σg )2 (1πu )4 configuration is missing in the active space. Although it is clear that these small ˆ (0) is surprising, active spaces are not optimally chosen, the severity of the problems for H Fink ˆ (0) . For O and small given the non-problematic behavior found for the other choices of H 2 ˆ (0) are much poorer than active spaces, we also observed that the PECs obtained with H Fink ˆ (0) or H ˆ (0) . In particular for CAS(2e, 2o), corresponding the PECs calculated with H Dyall Fock to a single reference, but still reasonable zeroth-order description of the triplet state close ˆ (0) leads to a repulsive potential at second order and strongly to equilibrium geometry, H Fink overbinds at third order, see the Supporting Information Fig. SI-81. The other two zerothorder Hamiltonian lead to a slightly more regular result, but also with strong deviations from the reference curve. The quality of the MRPT calculations has also been measured by the results for spectroscopic constants. The trend is in accordance with what is observed for the potential energy curves. For the full valence active space calculations, the root mean square deviation to the reference results is reported in Table 1. Among the MRPT methods, the smallest errors for ˆ (0) , but H ˆ (0) leads equilibrium distances and harmonic frequencies are obtained using H Fink Dyall to smaller NPEs on the average. For reduced active spaces, however, it is rather difficult to establish a general rule for the quality of the results, even for this restricted subset of first-group diatomic molecules. To show the spread of the deviations to the reference value for smaller active spaces, we give in Figure 3 a scatter plot, with the absolute value of the

24


Page 24 of 47

Page 25 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Square root mean deviation to the reference PEC (FCI or icMRCCSD(T) with full valence active space, with corresponging basis set) over the nine diatomic molecules, for the full valence active spaces, for equilibrium distance (Re ), harmonic frequency (ωe ), anharmonic constant (ωe xe ), and non-parallelity error (NPE). Method Basis set: MRCI MRCI+Q icMRCCSD(T) ˆ (0) MRPT2 H Fock ˆ (0) MRPT2 H Dyall ˆ (0) MRPT2 H Fink ˆ (0) MRPT3 H Fock ˆ (0) MRPT3 H Dyall

Re (pm) DZ TZ QZ 0.02 0.04 0.05 0.01 0.01 0.01 0.00 0.00 0.00 0.04 0.04 0.04 0.07 0.09 0.08

ωe DZ 1.20 0.48 0.15 4.15 5.61

(cm−1 ) TZ QZ 2.54 2.76 0.84 0.61 0.00 0.00 4.04 3.96 5.56 5.03

ωe xe (cm−1 ) DZ TZ QZ 0.12 0.66 0.75 0.02 0.51 0.60 0.07 0.00 0.00 0.21 0.74 0.86 0.47 0.26 0.37

NPE (kJ/mol) DZ TZ QZ 1.60 2.16 2.29 0.56 0.40 0.47 0.15 0.00 0.00 3.02 4.10 4.67 4.23 3.37 3.17

0.02 0.05 0.06 2.21 3.11 3.42 0.61 0.28 0.34 2.77 4.93 5.36 0.03 0.06 0.07 2.32 4.62 5.95 0.08 0.59 0.67 2.02 3.22 4.02 0.03 0.03 0.03 2.37 2.11 1.97 0.09 0.44 0.56 1.12 1.67 1.79

ˆ (0) 0.01 0.01 0.01 1.05 0.84 1.29 0.14 0.54 0.64 1.52 1.59 1.50 MRPT3 H Fink

deviation to the reference spectroscopic constants calculated with MRPT3 plotted against the number of configuration state functions (CSF) in the reference space, normalized such that full valence active spaces correspond to one and references with only one CSF correspond to zero. For comparison, also the results from icMRCCSD calculations are shown in this figure and it is evident that the perturbative methods clearly loose in robustness. Among the third-order methods, no clear picture emerges concerning the optimal choice for the zeroth-order Hamiltonian. The average errors (µ, shown in the legend of Figure 3) over ˆ (0) . However, all molecules and active spaces are quite similar for all the three choices for H ˆ (0) is much larger (except for ωe xe ). This comes from the the standard deviation (σ) for H Fink ˆ (0) gives the best results for most of the cases, it might give excepfact that, although H Fink tionally large deviation for some, here mainly for the small active spaces of O2 . The PECs of C2 that show intruder states are not included in this analysis, since the NPE and the ˆ (0) tends spectroscopic constants are not meaningful for this case. We also observe that H Fock ˆ (0) with respect to the reduction of the active space, in agreement to be more stable than H Dyall with the conclusions of Havenith and coworkers. 50 However, strong exceptions like O2 exist, ˆ (0) . leading to larger standard deviation of the results for H Fock

25



ref. PEC Re|, in pm

4

|

5

2

Re, in pm Average ( ) and standard deviation ( ) (0) ) MRPT3 (HFock = 0.70; = 0.94 (0) ) MRPT3 (HDyall = 0.67; = 0.75 (0) ) MRPT3 (HFink = 0.71; = 1.32 icMRCCSD = 0.24; = 0.32

ref. PEC

|

400 300 200 100

1 0

0 one CSF Full Valence 0.0 1.0 0.2 0.4 0.6 0.8 Number of reference CSF, normalized to the F. V. active space

8

exe, in cm 1 Average ( ) and standard deviation ( )

NPE, in kJ/mol Average ( ) and standard deviation ( )

600

(0) ) MRPT3 (HFock = 2.37; = 3.09 (0) ) MRPT3 (HDyall = 2.65; = 3.88 (0) ) MRPT3 (HFink = 1.41; = 2.38 icMRCCSD = 1.06; = 1.06


500 NPE|, in kJ/mol

10

one CSF Full Valence 0.0 1.0 0.2 0.4 0.6 0.8 Number of reference CSF, normalized to the F. V. active space

ref. PEC

6 4

400 300 200

|

ref. PEC

exe|, in cm 1

12


500

3

14

e, in cm 1 Average ( ) and standard deviation ( )

600

e|, in cm 1

6

|

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 47

100

2 0

0 one CSF Full Valence 0.0 0.2 0.4 0.6 0.8 1.0 Number of reference CSF, normalized to the F. V. active space

one CSF Full Valence 0.0 0.2 0.4 0.6 0.8 1.0 Number of reference CSF, normalized to the F. V. active space

Figure 3: Absolute value of the deviation to the reference value (FCI or icMRCCSD(T) with full valence active space) of main spectroscopic constants and the non-parallelity error (NPE) against the normalized number of configuration state functions (CSF) in the reference space, calculated as (nCSF − 1)/(nCSF, F.V. − 1), where nCSF is the number of CSFs in the reference space and nCSF, F.V. is the number of CSFs in the reference space for a full valence active space. The cc-pVDZ basis set was used.

26


Page 27 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.2

Ozone

We determined the geometry and harmonic frequencies of O3 , and calculated the potential energy curves for its symmetric and asymmetric dissociation. Three active spaces have been used: CAS (12e, 9o), CAS (8e, 7o), and CAS(2e, 2o). The selected orbitals are shown in the Supporting Information (page SI-163). Convergence of the harmonic frequencies with the basis set is known to be slow, 30 but MRCI+Q results with the cc-pVQZ basis set are very close to the experimental values. 51,52 We thus take the experimental spectroscopic constants as reference for calculations with the cc-pVQZ basis set. (We note that the experimental spectroscopic constants of Refs. 51,52 are corrected for anharmonic effects such that a comparison to computed equilibrium structures and harmonic frequencies is fully valid.) For the smaller active spaces, for which results are given in the Supporting Information, cc-pVDZ calculations were performed, and compared to MRCI+Q with CAS (12e, 9o) results for the same basis. For this small basis set, direct comparison to the experimental values is not appropriate. Table 2 shows the geometry and harmonic frequencies calculated with the cc-pVQZ basis set and the largest active space considered in this work, CAS (12e, 9o). For this case, the ˆ (0) accuracy in the geometry is similar for all variants, with slightly larger errors for H Fink at the second-order level. For harmonic frequencies, very good results are obtained with ˆ (0) , whereas H ˆ (0) and H ˆ (0) show stronger deviations, especially for the MRPT3 using H Dyall Fock Fink ω2 harmonic frequency corresponding to the antisymmetric stretching mode, which is known to be very challenging. 53 The results for smaller active spaces (see Supporting Information for details) can be summarized as follows: For CAS (8e, 7o), the geometry is well described by the three zeroth-order Hamiltonians, but for CAS (2e, 2o), the bond distance is strongly underestimated, with deviˆ (0) ). For harmonic frequencies, the results with these ations reaching 4 pm (MRPT3 with H Fock smaller active spaces become very poor, especially for CAS (2e, 2o) that can be considered unsuitable for the MRPT calculations for this system. We recall that icMRCCSD(T) and 27



icMRCCSDT calculations with CAS (2e, 2o), on the other hand, give accurate frequencies. 30 Table 2: Equilibrium geometries and harmoic frequencies for O3 , with basis cc-pVQZ and CAS (12e, 9o), given as the difference to the experimental values. Bond lenghts (R) in pm, bond angle in degree, and harmonic frequencies in cm−1 . R α 0.7 -0.1 0.2 -0.1 1.0 0.4 -0.1 -0.0 0.3 -0.1 0.5 -0.1 -0.2 0.1 0.2 -0.0 127.3 116.8

method ˆ (0) MRPT2 H Fock ˆ (0) MRPT2 H Dyall ˆ (0) MRPT2 H Fink ˆ (0) MRPT3 H Fock ˆ (0) MRPT3 H Dyall (0) ˆ MRPT3 HFink MRCI MRCI+Q Experimental

ω1 (A1 ) ω2 (B1 ) ω3 (A1 ) -13 -60 -26 -3 -18 11 -19 24 64 4 24 10 -3 9 -5 -3 59 -13 6 23 11 -2 6 -1 715 1087 1133

The potential energy curves for the symmetric dissociation of ozone, with the bond angle fixed at 142.76◦ , is shown in Figure 4. This is not the equilibrium bond angle, but the motivation for this choice is that the second-order approximate coupled-cluster (CC2) model fails to describe the equilibrium structure of ozone, and this failure is more evident in the curve for the symmetric dissociation at this angle. 54 Calculations with the icMRCCSD(T) method and CAS (12e, 9o) have also been carried out, and are used as reference to analyze the MRPT results. The MRPT energies behave very similarly to the PECs of the diatomic (0)

ˆ molecules: For the larger active spaces, MRPT2 with H Fink overshoots the correlation energy, ˆ (0) largely underestimates it. The curve closest to the icMRCCSD(T) reference whereas H Dyall (0)

ˆ is obtained by MRPT3 using H Fink . For the smallest active space, CAS (2e, 2o), an intruder ˆ (0) at a bond length close to 1.8 Å. Although the quality of the CAS (2e, state occurs for H Fink ˆ (0) or H ˆ (0) is also rather poor, they are smooth and do not 2o) curves calculated with H Fock Dyall show intruder state problems. We found a very similar behavior for the potential energy curves for the asymmetric dissociation of ozone, that is, the extension of one O–O bond, while keeping the bond angle and the other ROO distance fixed. Graphs are given in the Supporting Information, see in

28


Page 28 of 47

Page 29 of 47

icMRCCSD(T) CAS (12e, 9o) (0) ) (0) ) MRPT2 (HFock MRPT3 (HFock (0) ) (0) ) MRPT2 (HDyall MRPT3 (HDyall (0) ) (0) ) MRPT2 (HFink MRPT3 (HFink

224.70 224.75 Energy (Eh)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


224.80 224.85 CAS (12e, 9o)

224.90 1.2

1.4 1.6 ROO (Å)

1.81.2

CAS (8e, 7o) 1.4 1.6 ROO (Å)

1.81.2

CAS (2e, 2o) 1.4 1.6 ROO (Å)

1.8

Figure 4: Symmetric dissociation of Ozone, at bond angle of 142.76◦ . Calculations carried out with the cc-pVDZ basis set.

29



particular Fig. SI-136. We note that the correct description of this curve is a challenge for ab initio methods, as they tend to predict an unphysical barrier at about 1.8 Å. 55 It has been shown that only with large basis set and appropriate multireference treatment the correct behavior is obtained. 55 For the cc-pVDZ basis set, used for this case, even the icMRCCSD(T) calculations with CAS (8e, 7o) show a small barrier, and the MRPT calculations perform similar for the largest active space.

3.3

F + H2 reaction

The F + H2 −−→ FH + H reaction is an important prototype for chemical dynamics. Along with its isotopic variations, it has been widely studied by theoretical 56–62 and experimental 63–65 methods, such as state-to-state quantum dynamics over full dimensional ab initio potential energy surfaces and crossed molecular beam experiments, with very good agreement between the results provided from both sides. This system is small enough to allow ab initio calculations with highly correlated wave functions, like the CCSDTQ method and MRCI+Q with large active spaces, as employed in the benchmark studies of Werner, Kállay and Gauss. 47 These calculations give accurate predictions of the potential energy surface at the transition state region, and predict that this reaction has a low barrier (∼ 1.7 kcal mol−1 ), associated with a transition state with bent geometry, but the transition state restricted to the linear geometry has only a slightly higher energy (∼ 2 kcal mol−1 ). Although these values are not directly comparable to experiments, because they do not consider spin-orbit coupling and core-valence correlation effects, they can be used as reference to check the accuracy of MRPT. We made calculations with four active spaces, and used the basis set (aug-cc-pVTZ basis set 66 for the fluorine atom and the cc-pVTZ for the hydrogen atoms) and the “benchmark geometry” of Reference 47 (see Table I of that publication). The calculated barrier heights and the exothermicity of the reaction are given in Table 3. For the largest active spaces, CAS (7e, 7o) and CAS (9e, 10o), the barrier heights are overestimated by about 1 kcal mol−1 for all MRPT3 variants. This deviation is probably 30


Page 30 of 47

Page 31 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the best accuracy one can expect from a MRPT3 calculation for such a small barrier height. Also, the linear transition state is correctly predicted to have slightly higher energy than the bent transition state. For the smaller active spaces, CAS (3e, 3o) and the full valence active space CAS (9e, 6o), the MRPT3 methods overestimate the barrier height by 3 to 4 kcal mol−1 and the energetic order of the transition states is inverted. This behavior can also be observed by the potential energy curves of the bending potential, shown in the Supporting Information Figs. SI-141 to SI-149. Whereas the relative energies of the transition states calculated with the three zerothorder Hamiltonians agree very well among each other, larger differences are obtained for the (0)

ˆ exothermicity. For instance, the value from MRPT2 calculations with H Dyall and the largest active space overestimate the CCSDTQ and MRCI+Q predictions by about 3.5 kcal mol−1 . ˆ (0) and H ˆ (0) are much closer. For the full valence active space, CAS(9e, Results with H Fock Fink ˆ (0) and H ˆ (0) overestimate the reaction energy by 9 and 5.5 kcal mol−1 , respectively, 6o), H Dyall Fink ˆ (0) ) overshoots by 2 kcal mol−1 only. at second order, whereas the CASPT2 value (using H Fock ˆ (0) underestimate the CCSDTQ and At third order, on the other hand, calculations with H Fink MRCI+Q predictions by almost 2 kcal mol−1 even for the large CAS(9e, 10o) active space, ˆ (0) and H ˆ (0) practically match the reference values. This agreement is somewhat whereas H Fock Dyall fortuitous, since the accuracy of the MRPT calculations is at most 1 kcal mol−1 , based on what has been observed for the barrier heights. But it becomes clear that the performance of the three different zeroth-order Hamiltonians is not easy to rank. For comparison, we also carried out icMRCCSD and icMRCCSD(T) computations for all active spaces. The results, also shown in Table 3, again demonstrate the robustness of these methods with respect to active space choice. In particular the icMRCCSD(T) results are very accurate even for the smallest active space, CAS(3e,3o).

31



Page 32 of 47

Table 3: Relative energies, in kcal/mol, for the F + H2 reaction. Calculations carried out with the aug-cc-pVTZ basis set for Fluorine and cc-pVTZ basis set for Hydrogen and the “benchmark geometry” of Ref. 47. Method ˆ (0) H Fock ˆ (0) H Dyall ˆ (0) MRPT2 H Fink ˆ (0) MRPT3 H Fock ˆ (0) MRPT3 H Dyall ˆ (0) MRPT3 H Fink MRCI+Q icMRCCSD icMRCCSD(T) CCSDTQa MRCI+Q (9e,10o)a,b a Ref. 47. b The small difference MRPT2 MRPT2

Linear barrier Bent barrier Exothermicity (3e, 3o) (9e, 6o) (7e, 7o) (9e, 10o) (3e, 3o) (9e, 6o) (7e, 7o) (9e, 10o) (3e, 3o) (9e, 6o) (7e, 7o) (9e, 10o) 3.30 3.74

3.13 2.72

2.60 3.01

2.59 2.66

3.52 3.96

3.18 2.71

2.57 2.91

2.39 2.48

32.45 38.28

32.69 39.51

31.10 30.82

31.30 34.07

3.46 4.68 4.67

3.19 4.49 4.32

2.61 2.81 2.71

2.86 2.79 2.70

3.29 5.41 5.39

3.24 5.04 4.86

2.40 2.74 2.56

2.80 2.70 2.59

32.31 29.27 27.66

36.01 29.26 29.23

26.44 30.84 31.58

30.01 30.77 30.59

4.56 3.16 2.70 2.04

4.14 3.04 2.33 2.13

2.69 2.23 1.94 1.91 2.052 2.208

2.86 2.19 2.08 2.05

5.38 3.28 2.67 1.73

4.61 3.09 2.25 1.89

2.50 1.97 1.65 1.61 1.767 1.960

2.88 1.93 1.83 1.77

27.14 29.94 29.91 30.57

27.96 31.05 30.12 30.61 30.04 30.59 30.32 30.48 30.58 30.45

28.87 30.54 30.56 30.54

to the values in the present work comes from the fact that we are not using SA-CASSCF.

Figure 5: Isomerization pathway for the complex [Cu2 O2 ]2+ (NH3 )6 . Calculations have been carried out with the CAS (8e, 10o) active space. The MRCI+Q curve of Rode et al. has been calculated with a perturbative estimation of the correlation of some of the internal orbitals (see Ref. 67 for details), whereas the present MRCI+Q calculation fully includes the correlation of all orbitals.

32


Page 33 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.4

Copper-oxygen complex

Our last example is the isomerization of the dicopper oxygen complex [Cu2 O2 ]2+ (NH3 )6 , from the µ-η 2 :η 2 -peroxo to the bis(µ-oxo) structures. This is a model system to study potential catalysts, in particular for selective methane oxidation. 68 A number of catalysts for these processes have their active sites based on copper-oxygen complexes, like zeolites, 69,70 metal organic frameworks, 71 and enzymes. 72,73 Previous computational studies on this system 67,74 have shown that MRCI+Q predicts the µ-η 2 :η 2 -peroxo structure to be more stable than the bis(µ-oxo) isomer, whereas CASPT2 predicts the opposite. The relative energies obtained from density functional theory using the B3LYP functional qualitatively agree with the MRCI+Q results, but this was shown to strongly depend on the functional, in particular on the amount of exact exchange. 67 Given the higher accuracy of the MRCI+Q theory, it has been concluded that CASPT2 is not appropriate for this system. 67 In the MRCI calculations previously carried out by Rode and Werner, 67 the number of correlating orbitals was restricted to 32, a limitation in the MRCI implementation of Molpro at that time. Thus, their MRCI+Q results did not include the correlation of the NH orbitals, and the hybrid CIPT2 approach 75 was used to estimate the missing correlation perturbatively. We carried out further MRCI+Q calculations, now using the most recent implementation that allows the correlation of all orbitals. 76 These calculations were done with the same basis set as used before, 67,74 namely cc-pVDZ for Cu and H, aug-cc-pVDZ for N and O, and a 10-electron relativistic pseudopotential for the Cu atoms. 77 The active space has eight electrons distributed in ten orbitals, CAS (8e, 10o). For the isomerization pathway, we used optimized structures at the restricted Kohn-Sham level of theory employing the B3LYP functional. Although the isomerization pathway has been reoptimized at the CASPT2 level of theory in the previous works, 67,74 this does not lead to major differences in the isomerization curves. 74 The present MRCI+Q results, correlating all orbitals, show a slight increase of the stabilization of the µ-η 2 :η 2 -peroxo structure relative to the bis(µ-oxo) structure. However, the difference between the isomerization curves is rather small and the 33



present MRCI+Q results agree well with the curves reported by Rode and Werner. Calculations at the icMRCCSD level were also performed, but since the system size makes these calculations very demanding, only few points with the CAS (8e, 10o) have been obtained. Their relative energies agree surprisingly well with the original MRCI+Q curve of Rode and Werner, see the Supporting Information Fig. SI-150 and Tab. SI-113. Further icMRCCSD calculations using a much smaller active space, CAS (2e, 2o), have also been attempted, but they strongly underestimate the energy of the µ-η 2 :η 2 -peroxo structure and are not further considered. Based on these facts, we conclude that the pair of MRCI+Q curves gives references accurate enough to analyze the MRPT results. ˆ (0) , the MRPT2 The MRPT isomerization curves are shown in Figure 5. When using H Fock calculations predict the µ-η 2 :η 2 -peroxo structure to be less stable than the bis(µ-oxo) structure by almost 10 kcal mol−1 , in agreement with the previous CASPT2 calculations, 67,74 but ˆ (0) , on clearly in contrast with MRCI+Q computations. The MRPT3 result based on H Fock the other hand, strongly overstabilises the µ-η 2 :η 2 -peroxo isomer, still giving results differing from the MRCI+Q computations by 10 kcal mol−1 . The quality of the MRPT curves is ˆ (0) or H ˆ (0) are used. At second order, the two methods give largely improved when H Dyall Fink very similar results, but still overestimate the relative energy of the µ-η 2 :η 2 -peroxo structure compared to the MRCI+Q results by nearly 10 kcal mol−1 . Including third-order correcˆ (0) and H ˆ (0) give very good results, which are close to the MRCI+Q curves. tions, both H Dyall Fink These findings corroborate previous conclusions that second-order perturbation theory (in particular CASPT2) is not sufficiently accurate for describing this isomerization reaction. At third order, qualitatively better results are obtained, but the results still depend strongly ˆ (0) . on the choice of H

34


Page 34 of 47

Page 35 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4

Conclusions

Based on a partitioning of the electronic Hamiltonian into an appropriate zeroth-order part ˆ (0) and a first-order part W ˆ =H ˆ −H ˆ (0) giving rise to the major part of dynamic correH lation, we have analyzed the internally contracted multireference coupled-cluster (icMRCC) theory 27–29 through third order, resulting in an icMRCC perturbation expansion. This analysis is analogous to similar undertakings in the single-reference framework. 1,2,6 In the present work, we have not accounted for relaxation of the reference coefficients and we only consider non-degenerate perturbation theory. Up to second order, our analysis results in expressions that are equivalent to those from Rayleigh-Schrödinger perturbation theory based on a linear parametrization of the perturbed wave function. Depending on the choice for the zeroth-order Hamiltonian, the resulting icMRCC perturbation expansion to second order is equivalent to second-order complete active space perturbation theory (CASPT2, 9–11 using a projected Fock operator), N -electron valence perturbation theory (NEVPT2, 13,14 using Dyall’s Hamiltonian, 12 ) or retention of the excitation degree perturbation theory (REPT2, 15,16 using the zeroth-order Hamiltonian proposed by Fink 15 ). At third order, the formal result for the icMRCC perturbation expansion differs from that of Rayleigh-Schrödinger perturbation theory, as a number of extra terms appear. These terms have been analyzed in detail and it has been shown based on both analytical and numerical results that the extra contributions are very small (typically smaller than 0.1 mEh for the investigated cases and thus much smaller than the intrinsic error of the third-order energy). The cancelling terms come from different parts of the icMRCC equations (that is from energy and residual equations) and it will be important to consider them together in any approximate version of icMRCC theory. This study also provides a number of numerical examples, where second and third-order multireference perturbation theories (MRPT) have been compared using three different zeroth-order Hamiltonians, as mentioned above. The examples comprise a number of di35



Page 36 of 47

atomic molecules, the ozone moleclue, stationary points on the F + H2 −−→ FH + H reaction energy surface, and the isomerization pathway of a copper-oxygen complex. The MRPT results have been benchmarked against Full CI, icMRCCSD(T), and MRCI+Q results. In cases where large active space can be used, third-order MRPT using the Fink Hamiltonian often gives results that are rather close to reference values, as already observed by Sharma and Alavi. 17,18 Sometimes, however, despite the closeness to the reference energies, the resulting potential energy surface does not run sufficiently parallel to the reference and (0)

ˆ larger deviations are found for spectroscopic constants. For instance, using H Fink the MRPT3 value for the O3 antisymmetric stretch frequency deviates by nearly 60 cm−1 from the refer(0)

(0)

ˆ ˆ ence result, while MRPT3 using either H Dyall or HFock performs significantly better. Results for the F + H2 barrier heights and reaction energy show that there are limits to the expected accuracy of MRPT3 in general and deviations of 1 or 2 kcal mol−1 relative to high-accuracy reference calculations remain even for large double-shell active spaces. For smaller active spaces, the MRPT methods clearly deteriorate in comparison to the full icMRCC methods, as demonstrated by the results for F + H2 barrier heights and reaction energies, but also for ˆ (0) can be the other examples considered in this work. We also found that MRPT using H Fink rather prone to intruder state problems when the active space becomes too small, as exemplified for C2 and O3 . While these examples may seem artificial, it is clear that optimally sized active spaces are hard to achieve for large systems and a certain robustness against intruder state problems is thus important. We also point out that the single-reference limit ˆ (0) is the linearized coupled-cluster doubles (L-CCD) method, which of MRPT3 using H Fink clearly limits the accuracy of the correlation treatment of the inactive part of the system. ˆ (0) . Overall, no clear picture emerges concerning the accuracy of the three choices for H While the Fink Hamiltonian provides highly accurate results in many cases, there are clear outliers and other examples exist where either the projected Fock operator or the Dyall Hamiltonian perform better. In comparison to full icMRCCSD the accuracy of all thirdorder approaches remains limited and particularly the stability of coupled-cluster methods

36


Page 37 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with respect to active space choice is not reached. Nonetheless, we believe that the present work will be helpful in creating alternative approximate methods that can avoid a huge part of the complexity inherent to full-blown icMRCC theory.

Acknowledgement This work has been supported by the European Research Council (ERC) Advanced Grant 320723 (ASES). Y.A.A also thanks the grants #2017/21199-0 and #2018/04617-6, São Paulo Research Foundation (FAPESP), and to Márcio Fabiano da Silva for his kind support.

Supporting Information Available A comprehensive collection of graphs and tables with the results for the systems described in this article are given in the Supporting Information. This information is available free of charge via the Internet at http://pubs.acs.org.

References (1) Bartlett, R. J. Many-body perturbation theory and coupled cluster theory for electron correlation in molecules. Ann. Rev. Phys. Chem. 1981, 32, 359–401. (2) Bartlett, R. J. The coupled-cluster revolution. Mol. Phys. 2010, 108, 2905–2920. (3) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A fifth-order perturbation comparison of electron correlation theories. Chem. Phys. Lett. 1989, 157, 479–483. (4) Stanton, J. F. Why CCSD(T) works: a different perspective. Chem. Phys. Lett. 1997, 281, 130–134.

37



(5) Koch, H.; Christiansen, O.; Jørgensen, P.; de Merás, A. M. S.; Helgaker, T. The CC3 model: An iterative coupled cluster approach including connected triples. J. Chem. Phys. 1997, 106, 1808–1818. (6) Helgaker, T.; Jørgensen, P.; Olsen, J. Molecular electronic structure theory; Wiley, 2000. (7) Wolinski, K.; Sellers, H. L.; Pulay, P. Consistent generalization of the Møller-Plesset partitionining to open-shell and multiconfigurational SCF reference states in manybody perturbation theory. Chem. Phys. Lett. 1987, 140, 225–231. (8) Hirao, K. Multireference Møller-Plesset method. Chem. Phys. Lett. 1992, 190, 374–380. (9) Andersson, K.; Malmqvist, P.-A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. Second-order perturbation theory with a CASSCF reference function. J. Phys. Chem. 1990, 94, 5483–5488. (10) Andersson, K.; Malmqvist, P.-A.; Roos, B. O. Second-order perturbation theory with a complete active space self-consistent field reference function. J. Chem. Phys. 1992, 96, 1218–1226. (11) Werner, H.-J. Third-order multireference perturbation theory. The CASPT3 method. Mol. Phys. 1996, 89, 645–661. (12) Dyall, K. G. The choice of a zeroth-order Hamiltonian for second-order perturbation theory with a complete active space self-consistent-field reference function. J. Chem. Phys. 1995, 102, 4909–4918. (13) 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. (14) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. n-electron valence state perturbation theory:

38


Page 38 of 47

Page 39 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A spinless formulation and an efficient implementation of the strongly contracted and of the partially contracted variants. J. Chem. Phys. 2002, 117, 9138–9153. (15) Fink, R. F. Two new unitary-invariant and size-consistent perturbation theoretical approaches to the electron correlation energy. Chem. Phys. Lett. 2006, 428, 461–466. (16) Fink, R. F. The multi-reference retaining the excitation degree perturbation theory: A size-consistent, unitary invariant, and rapidly convergent wavefunction based ab initio approach. Chem. Phys. 2009, 356, 39–46. (17) Sharma, S.; Alavi, A. Multireference linearized coupled cluster theory for strongly correlated systems using matrix product states. J. Chem. Phys. 2015, 143, 102815. (18) Sharma, S.; Knizia, G.; Guo, S.; Alavi, A. Combining Internally Contracted States and Matrix Product States To Perform Multireference Perturbation Theory. J. Chem. Theory Comput. 2017, 13, 488–498. (19) Werner, H.-J.; Reinsch, E. A. An iterative CI method with configurations generated from several reference determinants. Proceedings of the fifth seminar on computational methods in quantum chemistry, Groningen, Holland 1981, (20) Werner, H.-J.; Reinsch, E. A. The self-consistent electron pairs method for multiconfiguration reference state functions. J. Chem. Phys. 1982, 76, 3144–3156. (21) Rintelman, J. M.; Adamovic, I.; Varganov, S.; Gordon, M. S. Multireference secondorder perturbation theory: How size consistent is “almost size consistent”. J. Chem. Phys. 2005, 122, 044105. (22) Lyakh, D. I.; Musiał, M.; Lotrich, V. F.; Bartlett, R. J. Multireference Nature of Chemistry: The Coupled-Cluster View. Chem. Rev. 2012, 112, 182–243. (23) Köhn, A.; Hanauer, M.; Mück, L. A.; Jagau, T.-C.; Gauss, J. State-specific multireference coupled-cluster theory. WIREs Comput. Mol. Sci. 2012, 10.1002/wcms.1120. 39



(24) Evangelista, F. A. Perspective: Multireference coupled cluster theories of dynamical electron correlation. J. Chem. Phys. 2018, 149, 030901. (25) Evangelista, F. A.; Simmonett, A. C.; Schaefer III, H. F.; Mukherjee, D.; Allen, W. D. A companion perturbation theory for state-specific multirefernece coupled cluster methods. Phys. Chem. Chem. Phys. 2009, 11, 4728–4741. (26) Giner, E.; Angeli, C.; Garniron, Y.; Scemama, A.; Malrieu, J.-P. A Jeziorski-Monkhorst fully uncontracted multi-reference perturbative treatment. I. Principles, second-order versions, and tests on ground state potential energy curves. J. Chem. Phys. 2017, 146, 224108. (27) Banerjee, A.; Simons, J. The coupled-cluster method with a multiconfiguration reference state. Int. J. Quantum Chem. 1981, 19, 207–216. (28) Evangelista, F.; Gauss, J. An orbital-invariant internally contracted multireference coupled cluster approach. J. Chem. Phys. 2011, 134, 114102. (29) Hanauer, M.; Köhn, A. Pilot applications of internally contracted multireference coupled cluster theory, and how to choose the cluster operator properly. J. Chem. Phys. 2011, 134, 204111. (30) Hanauer, M.; Köhn, A. Perturbative treatment of triple excitations in internally contracted multireference coupled cluster theory. J. Chem. Phys. 2012, 136, 204107. (31) Finley, J.; Malmqvist, P.-Å.; Roos, B.; Serrano-Andrés, L. The multi-state CASPT2 method. Chem. Phys. Lett. 1998, 288, 299–306. (32) Angeli, C.; Borini, S.; Cestari, M.; Cimiraglia, R. A quasidegenerate formulation of the second order n-electron valence state perturbation theory approach. J. Chem. Phys. 2004, 121, 4043–4049.

40


Page 40 of 47

Page 41 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(33) Angeli, C.; Borini, S.; Cavallini, A.; Cestari, M.; Cimiraglia, R.; Ferrighi, L.; Sparta, M. Developments in the n-electron valence state perturbation theory. Int. J. Quantum Chem. 2005, 106, 686–691. (34) Granovsky, A. A. Extended multi-configuration quasi-degenerate perturbation theory: The new approach to multi-state multi-reference perturbation theory. J. Chem. Phys. 2011, 134, 214113. (35) Shiozaki, T.; Győrffy, W.; Celani, P.; Werner, H. Communication: Extended multi-state complete active space second-order perturbation theory: Energy and nuclear gradients. J. Chem. Phys. 2011, 135, 081106. (36) Angeli, C.; Pastore, M.; Cimiraglia, R. New perspectives in multireference perturbation theory: the n-electron valence state approach. Theor. Chem. Acc. 2007, 117, 743–754. (37) Meyer, W. In Methods of Electronic Structure Theory; Schaefer III, H. F., Ed.; Plenum, 1977. (38) Siegbahn, P. E. M. Direct configuration interaction with a reference state composed of many reference configurations. Int. J. Quantum Chem. 1980, 18, 1229–1242. (39) van Dam, H. J. J.; van Lenthe, J. H.; Pulay, P. The size consistency of multi-reference Møller Plesset perturbation theory. Mol. Phys. 1998, 93, 431–439. (40) van Dam, H. J. J.; van Lenthe, J. H.; Ruttink, P. J. A. Exact size consistency of multireference Møller–Plesset perturbation theory. Int. J. Quantum Chem. 1999, 72, 549–558. (41) Li, C.; Evangelista, F. A. Multireference Driven Similarity Renormalization Group: A Second-Order Perturbative Analysis. J. Chem. Theory Comput. 2015, 11, 2097–2108. (42) Li, C.; Evangelista, F. A. Driven similarity renormalization group: Third-order multireference perturbation theory. J. Chem. Phys. 2017, 146, 124132. 41



(43) Dunning Jr., T. H. Gaussian basis set for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (44) Woon, D. E.; Dunning Jr., T. H. Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. J. Chem. Phys. 1993, 98, 1358– 1371. (45) Werner, H.-J.; Knowles, P. J. An efficient internally contracted multiconfigurationreference configuration interaction method. J. Chem. Phys. 1988, 89, 5803–5814. (46) 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. (47) Werner, H.-J.; Kállay, M.; Gauss, J. The barrier height of the F+H2 reaction revisited: Coupled-cluster and multireference configuration-interaction benchmarck calculations. J. Chem. Phys. 2008, 128, 034305. (48) MOLPRO, version 2015 a package of ab initio programs. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, and others, see http://www.molpro.net. (49) Menezes, F.; Kats, D.; Werner, H.-J. Local complete active space second-order perturbation theory using pair natural orbitals (PNO-CASPT2). J. Chem. Phys. 2016, 145, 124115. (50) Havenith, R. W. A.; Taylor, P. R.; Angeli, C.; Cimiraglia, R.; Ruud, K. Calibration of the n-electron valence state perturbation theory approach. J. Chem. Phys. 2004, 120, 4619–4625. (51) Tyuterev, V. G.; Kochanov, R. V.; Tashkun, S. A.; Holka, F.; Szalay, P. G. New analytical model for the ozone electronic ground state potential surface and accurate ab initio vibrational predictions at high energy range. J. Chem. Phys. 2013, 139, 134307.

42


Page 42 of 47

Page 43 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(52) Barbe, A.; Chichery, A.; Cours, T.; Tyuterev, V. G.; Plateaux, J. J. Update of the anharmonic force field parameters of the ozone molecule. J. Mol. Struct. 2002, 616, 55–65. (53) Stanton, J. F.; Lipscomb, W. N.; Magers, D. H.; Bartlett, R. J. Highly correlated single-reference studies of the O3 potential surface. I. Effects of high order excitations on the equilibrium structure and harmonic force field of ozone. J. Chem. Phys. 1989, 90, 1077–1082. (54) Pabst, M.; Köhn, A.; Gauss, J.; Stanton, J. F. A worrisome failure of the CC2 coupledcluster method when applied to ozone. Chem. Phys. Lett. 2010, 495, 135–140. (55) Dawes, R.; Lolur, P.; Li, A.; Jiang, B.; Guo, H. Communication: An accurate global potential energy surface for the ground electronic state of ozone. J. Chem. Phys. 2013, 139, 201103. (56) Manolopoulos, D. E.; Stark, K.; Werner, H.-J.; Arnold, D. W.; Neumark, D. M. The transition state of the F + H2 reaction. Science 1993, 262, 1852–1855. (57) Castillo, J. F.; Hartke, B.; Werner, H. J.; Aoiz, F. J.; Bañares, L.; Martínez-Haya, B. Quantum mechanical and quasiclassical simulations of molecular beam experiments for the F + H2 → HF + H reaction on two ab initio potential energy surfaces. J. Chem. Phys. 1998, 109, 7224–7237. (58) Stark, K.; Werner, H.-J. An accurate multireference configuration interaction calculation of the potential energy surface for the F + H2 → HF + H reaction. J. Chem. Phys. 1996, 104, 6515–6530. (59) Castillo, J. F.; Manolopoulos, D. E. Quantum mechanical angular distributions for the F + HD reaction. Faraday Discuss. 1998, 110, 119–138.

43



(60) Castillo, J. F.; Manolopoulos, D. E.; Stark, K.; Werner, H.-J. Quantum mechanical angular distributions for the F + H2 reaction. J. Chem. Phys. 1996, 104, 6531–6546. (61) Che, L.; Ren, Z.; Xingan,; Dong, W.; Dai, D.; Wang, X.; Zhang, D. H.; Yang, X.; Sheng, L.; Li, G.; Werner, H.-J.; Lique, F.; Alexander, M. H. Breakdown of the BornOppenheimer Approximation in the F + o-D2 → DF + D Reaction. Science 2007, 317, 1061–1064. (62) Alexander, M. H.; Manolopoulos, D. E.; Werner, H.-J. An investigation of the F + H2 reaction based on a full ab initio description of the open-shell character of the F(2 P) atom. J. Chem. Phys. 2000, 113, 11084–11100. (63) Neumark, D. M.; Wodtke, A. M.; Robinson, G. N.; Hayden, C. C.; Shobatake, K.; Sparks, R. K.; Schafer, T. P.; Lee, Y. T. Molecular beam studies of the F + D2 and F + HD reactions. J. Chem. Phys. 1985, 82, 3067–3077. (64) Neumark, D. M.; Wodtke, A. M.; Robinson, G. N.; Hayden, C. C.; Lee, Y. T. Molecular beam studies of the F + H2 reaction. J. Chem. Phys 1985, 82, 3045–3066. (65) Qiu, M.; Ren, Z.; Che, L.; Dai, D.; Harich, S. A.; Wang, X.; Yang, X.; Xu, C.; Xie, D.; Gustafsson, M.; Skodje, R. T.; Sun, Z.; Zhang, D. H. Observation of Feshbach Resonances in the F + H2 → HF + H Reaction. Science 2006, 311, 1440–1443. (66) Kendal, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of first-row atoms revisited. Systematic basis sets and wave functions. J. Chem. Phys. 1992, 96, 6796. (67) Rode, M. F.; Werner, H.-J. Ab initio study of the O2 binding in dicopper complexes. Theor. Chem. Acc. 2005, 114, 309–317. (68) Olivos-Suarez, A. I.; Szécsényi, À.; Hensen, E. J. M.; Ruiz-Martinez, J.; Pidko, E. A.; Gascon, J. Strategies for the Direct Catalytic Valorization of Methane Using Heterogeneous Catalysis: Challenges and Opportunities. ACS Catal. 2016, 6, 2965–2981. 44


Page 44 of 47

Page 45 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(69) Groothaert, M. H.; Smeets, P. J.; Sels, B. F.; Jacobs, P. A.; Schoonheydt, R. A. Selective Oxidation of Methane by the Bis(µ-oxo)dicopper Core Stabilized on ZSM-5 and Mordenite Zeolites. J. Am. Chem. Soc. 2005, 127, 1394–1395. (70) Alayon, E. M. C.; Nachtegaal, M.; Bodi, A.; Ranocchiari, M.; van Bokhoven, J. A. Bis(µ-oxo) versus mono(µ-oxo)dicopper cores in a zeolite for converting methane to methanol: an in situ XAS and DFT investigation. Phys. Chem. Chem. Phys. 2015, 17, 7681–7693. (71) Ikuno, T.; Zheng, J.; Vjunov, A.; Sanchez-Sanchez, M.; Ortunño, M. A.; Pahls, D. R.; Fulton, J. L.; Camaioni, D. M.; Li, Z.; Ray, D.; Mehdi, B. L.; Browning, N. D.; Farha, O. K.; Hupp, J. T.; Cramer, C. J.; Gagliardi, L.; Lercher, J. A. Methane Oxidation to Methanol Catalyzed by Cu-Oxo Clusters Stabilized in NU-1000 Metal-Organic Framework. J. Am. Chem. Soc. 2017, 139, 10294–10301. (72) Rosenzweig, A. C. The metal centres of particulate methane mono-oxygenase. Biochem. Soc. Trans. 2008, 36, 1134–1137. (73) Chan, S. I.; Yu, S. S.-F. Controlled Oxidation of Hydrocarbons by the MembraneBound Methane Monooxygenase: The Case for a Tricopper Cluster. Acc. Chem. Res. 2008, 41, 969–979. (74) Flock, M.; Pierloot, K. Theoretical Study of the Interconversion of O2 -Binding Dicopper Complexes. J. Phys. Chem. A 1999, 103, 95–102. (75) Celani, P.; Stoll, H.; Werner, H.-J.; Knowles, P. The CIPT2 method: Coupling of multireference configuration interaction and multi-reference perturbation theory. Application to the chromium dimer. Mol. Phys. 2004, 102, 2369–2379. (76) Shamasundar, K. R.; Knizia, G.; Werner, H.-J. A new internally contracted multireference configuration interaction method. J. Chem. Phys. 2011, 135, 054101.

45



(77) Figgen, D.; Rauhut, G.; Dolg, M.; Stoll, H. Energy-consistent pseudopotentials for group 11 and 12 atoms: adjustment to multi-configuration Dirac-Hartree-Fock data. Chem. Phys. 2005, 311, 227–244.

46


Page 46 of 47

Page 47 of 47

Perturbation expansion of internally contracted coupled-cluster theory up to third order Graphical abstract

Graphical TOC Entry

Wˆ + 0)

H ˆ( =

ˆ (0)

H

 (0) ˆ   HFock ˆ (0) = H Dyall   ˆ (0) HFink

Relative Energy (kcal/mol)

ˆ

|ΨicMRCC i = eT |ΨCASCI i

Hˆ

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


E (2)

10

E (2)

0

E (3)

10

E (3)

20 1.4

E = E (0) +E (1) +E (2) +E (3) +. . .

E (2)

1.6

E (3) 1.8 2.0 ROO (Å) (0) MRPT2 HFock (0) MRPT2 HFink (0) MRPT2 HDyall

1

47


2.2

2.4

(0) MRPT3 HFock (0) MRPT3 HFink (0) MRPT3 HDyall

Perturbation Expansion of Internally Contracted Coupled-Cluster

Recommend Documents