Subscriber access provided by - Access paid by the | UCSB Libraries
Quantum Electronic Structure
Correlation energy from the adiabatic connection formalism for complete active space wavefunctions Ewa Pastorczak, and Katarzyna Pernal J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00213 • Publication Date (Web): 22 May 2018 Downloaded from http://pubs.acs.org on May 29, 2018
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 34 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
Correlation energy from the adiabatic connection formalism for complete active space wavefunctions Ewa Pastorczak and Katarzyna Pernal∗ Institute of Physics, Lodz University of Technology, ul. Wolczanska 219, 90-924 Lodz, Poland E-mail:
[email protected]
Abstract Recently the adiabatic connection (AC) formula for the electron correlation energy has been proposed for a broad class of multireference wavefunctions. 1 We show that the AC formula used together with the extended random phase approximation (ERPA) can be efficiently applied to complete active space (CAS) wavefunctions to recover the remaining electron correlation. Unlike most of the perturbation theory approaches, the proposed AC-CAS method does not require construction of higher than two-electron reduced density matrices, which offers an immediate computational saving. In addition, we show that typically the AC-CAS systematically reduces the errors of both the absolute value of energy and of the energy differences (energy barrier) upon enlarging active spaces for electrons and orbitals. AC-CAS consistently gains in accuracy from including more active electrons. We also propose and study the performance of the AC0 approach resulting from the first-order expansion of the AC integrand at the coupling ∗
To whom correspondence should be addressed
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
constant equal to 0. AC0 avoids solving the full ERPA eigenequation, replacing it with small-dimension eigenproblems, while retaining good accuracy of the AC-CAS method. Low computational cost, compared to AC-CAS or perturbational approaches, makes AC0 the most efficient ab initio approach to capturing electron correlation for the CAS wavefunctions.
1
Introduction
Accounting for dynamic correlation in problems of multireference (MR) character has long been a challenge for quantum chemistry methods. The most prominent way of achieving this are the many-body perturbation theories (MBPT), in particular the complete active space perturbation theory (CASPT2) 2,3 and the second-order n-electron perturbation theory (NEVPT2). 4,5 CASPT2 (whose computational cost scales with the 8th power of the number of active orbitals) currently the most popular of those approaches, is plagued with two problems: the intruder states 6 and the lack of size-consistency. 7 The intruder-state problem is usually solved by level shifting 8 but this in turn often introduces arbitrariness to the solution. The size-inconsistency usually does not lead to serious problems, as for CASPT2 the error associated with this property is quite insignificant. NEVPT2 on the other hand is free of the mentioned problems thanks to the use of Dyall’s zeroth-order Hamiltonian 9 but the scaling remains similar (or higher) to that of CASPT2. Unfavorable scaling is a consequence of the presence of three- and four-particle reduced density matrices (RDMs) in the perturbation correction. The cost of their construction using the CAS expansion coefficients scales exponentially with the number of active orbitals. While it is possible to use cumulantbased expressions to approximate the RDMs, the consequences of it include loss of accuracy and appearance of the intruder states. 10 Recently, Chan and coworkers have introduced a time-dependent variant of NEVPT2, 11,12 which does not rely on 3- and 4-particle RDMs, and scales with the 6th power of the number of active orbitals. This is a significant improvement, but taking into account also the dependence of the computational cost on exponentially 2
ACS Paragon Plus Environment
Page 2 of 34
Page 3 of 34 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
growing number of determinants in CAS, the method remains expensive. Another option is DSRG-MRPT2, a perturbation theory applied to driven similarity renormalization group (DSRG) of Evangelista et al., 13 which uses only three-body density cumulant, but its scaling with the number of secondary orbitals is quite unfavorable. Alternatives to perturbationbased methods of describing the dynamic correlation in MR approaches are rather scarce. The existing MR coupled-cluster approaches are still far from settled and they require a lot of experience from their users. 14,15 The more pragmatic and much more computationally attractive approaches employ the density functional theory (DFT) to account for dynamic correlation. 16,17 Recently, one of us proposed a different path: an adiabatic connection (AC) between a system described exactly by the multireference wavefunction of choice (given in a form McWeeny’s 18 group product functions) and the system of interest. 1 The AC formalism has roots in the independent electron picture and it first appeared in the literature in the context of DFT. 19 The AC formalism together with random phase approximation and its variants has led to an outburst of a new family of orbital-dependent correlation functionals. 20–22 The AC combined with the extended random phase approximation 23 provides a way of obtaining dynamic correlation energy also for the multireference wavefunctions. 1 In this paper the AC is applied specifically to CAS wavefunction. The resulting method, AC-CAS, scales favorably with the number of active orbitals and is able to describe problems where both multireference effects and dynamic correlation play a role. Moreover, a computationally much more efficient version of the method, AC0-CAS (originating from the linear term in the expansion of the AC integrand for the coupling constant equal to zero), turns out to be as accurate as AC-CAS. In the following sections AC for CAS wavefunctions will be introduced and the explicit expression for AC0-CAS (previously only demostrated for a generalized valence bond wavefunction) will be presented. We will give the details of derivation and implementation of both methods and we will assess both their accuracy and cost.
3
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
2 2.1
Page 4 of 34
Theory AC-CAS: The Adiabatic Connection construction and Extended Random Phase Approximation for Complete Active Space Wavefunctions
Recently we have shown how to formulate the adiabatic connection correction to energy for multireference wavefunctions which can be presented as group product functions. 1 The complete active space (CAS) wavefunctions fall in this class since any CAS function can be written as a product of two functions ΨCAS = ψˆ1† ψˆ2† |vaci
,
(1)
where the operator ψˆ1† creates a single-reference Ninact - electron state with Minact = Ninact /2 inactive orbitals ψˆ1† = a ˆ†q1 a ˆ†q2 . . . a ˆ†qN
inact
,
(2)
whereas ψˆ2† creates an Nact -electron multireference state ψˆ2† =
active X
DQ a ˆ†q1 a ˆ†q2 . . . a ˆ†qNact
,
(3)
Q
where Q is a string of Nact indices corresponding to spinorbitals constructed from Mact active orbitals. A sum Ninact + Nact = N corresponds to a number of electrons in a system. If M is the number of basis set functions then Msec = M − Mact − Minact is equal to the number of secondary (unoccupied) orbitals. A deviation of the CAS energy from the exact value will be called a correlation energy, namely D
Ecorr = Eexact − Ψ
4
CAS
ˆ CAS |H|Ψ
ACS Paragon Plus Environment
E
,
(4)
Page 5 of 34 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
ˆ is a fully-interacting Hamiltonian where H
ˆ = H
X
hpq a ˆ†p a ˆq +
pq
1X † † a ˆa ˆa ˆq a ˆp hrs|pqi 2 pqrs r s
.
(5)
Both wavefunctions ψ1 = ψˆ1† |vaci and ψ2 = ψˆ2† |vaci resulting from the self-consistent field calculation satisfy a pertinent eigenequation
ˆ I ψI = EI ψI H
,
(6)
ˆ I involves an effective one-electron operator and a two-electron where a group Hamiltonian H part
ˆI = H
X
f † hef ˆp a ˆq + pq a
pq∈I
∀pq∈I
f hef pq = hpq +
XX
1 X † † a ˆa ˆa ˆq a ˆp hrs|pqi 2 pqrs∈I r s
nr [hpr|qri − hpr|rqi]
.
(7) (8)
J6=I r∈J
Indices I, J should be understood as sets of inactive or active spinorbitals, and it is assumed from now on that p, q, r, s are indices of natural spinorbitals pertaining to the CAS wavefunction, i.e. ∀pq where ∀p
γpq = ΨCAS |ˆ a†q a ˆp |ΨCAS = np δpq
,
(9)
0 ≤ np ≤ 1.
The adiabatic connection (AC) formalism proposed in Ref. 1 is based on the definition of the AC Hamiltonian
∀0≤α≤1
ˆα = H ˆ (0) + αH ˆ0 H ˆ0 = H ˆ −H ˆ (0) H
,
(10) (11)
for a coupling constant α. The zeroth-order Hamiltonian is a sum of group-Hamiltonians,
5
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
Page 6 of 34
Eq.(7), ˆ (0) = H
X
ˆI H
,
(12)
I
where in the case of the CAS wavefunction there are only two terms: one corresponding to the inactive orbital space, the other to the active one. Notice that the zeroth-order Hamiltonian, cf. Eq.(12) and Eq.(7), includes two-electron interactions within the active orbitals, similarly to the Dyall’s Hamiltonian employed as a zeroth-order Hamiltonian in the ˆ (0) assumed here treats inactive and active NEVPT2 perturbation method. 4,5 However, H spaces on equal footing so in general it is different from the Dyall’s Hamiltonian. Let Ψα be an eigenstate of the AC Hamiltonian
∀0≤α≤1
ˆ α Ψα = E α Ψα H
.
(13)
By varying the coupling parameter between 0 and 1 one smoothly switches between uncorrelatedgroup limit of α = 0 (no correlation among active, inactive, and secondary sets)
Ψ
α=0
CAS
=Ψ
,
E
α=0
D
CAS
= Ψ
E X (0) CAS ˆ |H |Ψ = EI
(14)
I
and fully-correlated-group limit of α = 1 (full correlation included)
Ψα=1 = Ψexact ,
E α=1 = Eexact
.
(15)
By exploiting the Hellman-Feynman theorem and the exact relation between the twoelectron reduced density matrix and one-electron reduced functions 24,25 a general AC correlation energy formula can be derived applicable to ground and excited states. The AC expression for the correlation energy pertaining to a ground state described by the CAS
6
ACS Paragon Plus Environment
Page 7 of 34 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
reference function Eq.(1) reads 1 ˆ AC Ecorr
1
(W α + ∆α ) dα ,
=
(16)
0
1 X Wα = 2 pqrs
! 0
X
α,0ν α,ν0 γpr γqs + (np − 1)nq δrq δps
hrs|pqi
,
(17)
ν
and
∆α =
XXXX
+
XXXX
I
I
α hpq γqp
p∈I J6=I q∈J α nr [hpr|qri − hpr|rqi] (δpq np − γqp ) ,
(18)
pq∈I J6=I r∈J
α ˆp |Ψα0 stands for one-electron reduced density matrix element correa†q a = Ψα0 |ˆ where γpq
α,0ν ˆp |Ψαν connect a†q a = Ψα0 |ˆ sponding to a ground state Ψα0 , transition density matrix elements γpq ˆ α . I,J denote sets of inactive, a ground and a νth excited state of the AC Hamiltonian H active, and secondary spinorbitals, and a prime sign in Eq.(17) indicates that terms for which the spinorbitals pqrs belong to the same set I are excluded. Therefore, if sets of inactive and secondary orbitals are empty and all orbitals are active (this would be a limiting case when the CAS wavefunction is equivalent to the full CI one) the correlation energy is equal to zero, as it should. In the first approximation one assumes that one-electron density matrix stays constant along the AC path ∀0≤α≤1
α=0 α γpq = γpq = δpq np
,
(19)
which results in vanishing of the ∆α term in Eq.(16) (its first line vanishes because indices p and q belong to different groups and therefore Kronecker delta in Eq.(19) would be equal to zero). Such approximation is justified if the reference wavefunction is of the multireference character since a major part of the static correlation would be already accounted for at α = 0. In this limit the structure of the occupation numbers vector is therefore correct
7
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
Page 8 of 34
and varying α from 0 to 1 would not change it significantly. In other words, for the multireference wavefunction the AC correction is responsible for accounting for mainly dynamic correlation effects, which do not alter one-electron density matrix, as discussed in Ref. 26. The second approximation to the exact AC correlation energy expression affects transition density matrices γ α,0ν . It has been proposed in Ref. 1 to employ Extended Random Phase Approximation (ERPA) equations 23,27 to approximate the latter. ERPA originates from ˆ† the Rowe’s equation of motion formalism 28 and it assumes that an excitation operator O ν creating an excited state |νi from the ground state |0i
ˆ † |0i = |νi O ν
(20)
includes only singly exciting particle-conserving operators, namely
ˆ ν† = O
X
Xpq a ˆ†p a ˆq + Ypq a ˆ†q a ˆp
,
(21)
p>q
where in the adopted notation a ˆ†p , a ˆp are, respectively, creation and annihilation operators pertaining to natural spinorbitals for the ground state. Notice that comparing with Ref. 23 ˆp are excluded from the excitation operator since it has been observed that diagonal terms a ˆ†p a their contribution to excited states is practically negligible. 27 By repeating derivation of the ˆ α given in Eq.(10) one arrives at the following ERPA equations for the AC Hamiltonian H generalized eigenproblem
α
α
Xαν
Xαν
−N 0 A B α = ων B α Aα Yνα 0 N Yνα
,
(22)
where ∀p>q
r>s
Npqrs = (np − nq )δpr δqs
8
ACS Paragon Plus Environment
.
(23)
Page 9 of 34 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
The α-dependent matrices Aα and B α read D E ˆ α, a [Aα ]pqrs = [B α ]pqsr = [ˆ a†p a ˆq , [H ˆ†s a ˆr ]]
,
(24)
so they are given solely in terms of the one- and two-electron reduced density matrices. For a given α the latter should correspond to Ψα0 but we propose to use in the α-ERPA equations (22) one- and two-electron density matrices obtained from the reference CAS wavefunction. This leads to the following explicit expressions for the elements of the α-ERPA matrices
[Aα ]pqrs = (hαqs δpr − hαpr δqs )(np − nq ) X X 1X α 1X α α α gurtp Γqtsu + gqtsu Γrupt + grqtu Γtups + g Γqrtu + 2 tu 2 tu tusp tu tu 1 X α 1 X α gtrwu Γwupt + δpr gtuws Γqwtu , + δqs 2 2 twu tuw
(25)
where
f hαpq = αhpq + δIp Iq (1 − α)hef , pq α gpqrs = α + δIr Is δIs Ip δIp Iq (1 − α) (hpq|rsi − hpq|sri)
(26) (27)
and two-electron reduced density matrix elements Γpqrs read
ˆq a ˆp |ΨCAS Γpqrs = ΨCAS |ˆ a†r a ˆ†s a
.
(28)
A symbol Ip denotes a set of inactive, active or secondary spinorbitals, which a spinorbital p belongs to. Eigenvectors corresponding to positive eigenvalues ων are normalized as follows [Yνα ]T N Yαν − [Xαν ]T N Xαν =
1 2
(29)
and they provide approximation to transition density matrix elements for the νth excitation,
9
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
Page 10 of 34
namely 29,30
∀p>q
γ α,0ν
∀q>p
γ α,0ν
qp
= (np − nq ) [Yνα ]pq
,
(30)
qp
= (np − nq ) [Xαν ]qp
.
(31)
Eigenvalues of the ERPA equations (22) come in pairs, i.e. for each positive eigenvalue ωνα there exists a negative counterpart equal to ωνα up to the absolute value. Since transition density matrices γ α,0ν pertain to excitations from a ground state then the only solutions of the ERPA equations corresponding to positive eigenvalues are exploited in Eqs.(30)-(31). Employing Eqs.(30)-(31) in the AC correlation energy expression, Eq.(16), and applying the approximation shown in Eq.(19) leads to following spin-summed approximate AC correlation energy expression
ˆ AC Ecorr
1
W α dα ,
=
(32)
0
where the AC integrand takes form
Wα = 2
X
0
p>q,r>s
{(np − nq )(nr − ns )
X
([Yνα ]pq − [Xαν ]pq )([Yνα ]rs − [Xαν ]rs )
ν
1 − [np (1 − nq ) + nq (1 − np )]δpr δqs } hpr|qsi 2
.
(33)
A prime indicates that all four indices pqrs do not belong at the same time to the set of active orbitals (i.e. terms for which p, q, r, s are active are not included in the sum). The same is true if pqrs are all inactive - such terms are simply vanishing because (np − nq )(nr − ns ) and AC np (1 − nq ) + nq (1 − np ) are all equal to 0. Thus, Ecorr introduces active-inactive-secondary
orbital correlation. Notice that if sets of inactive and secondary orbitals were empty (which corresponds to the case when CAS reference is equivalent to the full CI wavefunction) the approximate AC correlation energy would vanish by construction. Thus, one expects that AC when spaces of active orbitals and electrons are expanded, the value of Ecorr would decrease
and ultimately achieve 0. 10
ACS Paragon Plus Environment
Page 11 of 34 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
2.2
AC0-CAS: first-order expansion of the AC integrand at α = 0
As already noted in Ref. 1 multireference character of the wavefunction suggests that coupling among inactive, active and secondary orbitals will be weak which should be manifested in the nearly-linear behavior of the AC integrand given in Eq.(33). Were the AC integrand a perfectly linear function of α then its first-order expansion dW α α=0 ˜ α W =W + dα α=0
(34)
would obviously be equal to W α . As it will be discussed below, the linearization of the AC integrand allows one to reduce the computational cost of the AC method while retaining its good performance. This has been already observed for the GVB (generalized valence bond perfect-pairing) reference wavefunction in Ref. 1 and here we will demonstrate reliability of the linearization of W α for the CAS reference. To find an explicit expression of the linearized AC integrand as shown in Eq.(34) first of all, one observes that the zeroth-order term W α=0 vanishes for the CAS wavefunction. To show it one notices that at α = 0 the matrix Aα=0 acquires a block-diagonal structure with the following nonzero blocks
Aact :
α=0 Apqrs p>q,r>s
,
(35)
pqrs∈act
Aact−inact :
α=0 Apiqj pq∈act = δij Aα=0 piqi pq∈act ij∈inact
Asec−act :
α=0 Aapbq ab∈sec = δab Aα=0 apaq ab∈sec pq∈act
Asec−inact :
,
(36)
ij∈inact
,
(37)
pq∈act
α=0 Aaibj ab∈sec = δab δij Aα=0 aiai ab∈sec ij∈inact
,
(38)
ij∈inact
(where act, inact, and sec denote sets of active, inactive, and secondary orbitals, respectively) and in addition the elements of the B α matrix vanish
B α=0 = 0 ,
11
ACS Paragon Plus Environment
(39)
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
Page 12 of 34
with the exception of the act block. From now on assume that spinorbitals are in a descending order with respect to their occupancies, i.e. ∀p>q np − nq ≤ 0. The matrix Aα=0 is symmetric and the normalization condition (taking into account that Yνα=0 = 0) given in Eq.(29) yields
∀p>q
(np − nq )
X α=0 1 Xα=0 X = − δpr δqs ν ν pq rs 2 ν
.
(40)
Consequently, after employing Eq.(40) in Eq.(33) the AC integrand at α = 0 vanishes
W α=0 = 2
X
0
np (nq − 1) hpp|qqi = 0 ,
(41)
p>q
[it has been used that only terms corresponding to active-inactive (np < 1, ni = 1), secondaryactive (na = 0, nq < 1), and secondary-inactive (na = 0, ni = 1) terms enter the summation]. in the AC integrand expansion can be found after The first-order term W (1) = dW dα α=0 noticing that elements of the main matrix in the α-ERPA equation, cf. Eqs.(22)-(27), are linear in the coupling parameter α
Aα = Aα=0 + αA(1)
,
(42)
B α = B α=0 + αB (1)
,
(43)
and assuming that A(1) , B (1) are small perturbations. By applying the first-order perturbation theory to the α-ERPA problem one finds first-order α-expansion of the eigenvectors, namely
Xαν = Xα=0 + X(1) ν ν α ,
(44)
Yνα = Yνα=0 + Yν(1) α .
(45)
Using explicit expressions for the first-order corrections to the ERPA eigenvectors 29 it is straightforward to show that the first-order term in the expansion of the AC integrand with
12
ACS Paragon Plus Environment
Page 13 of 34 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
respect to α in Eq.(34) reads
W (1) = 4
X
0
(np − nq )(nr − ns )
p>q,r>s
(1) − (1) X Z+ + B (1) )Z+ − B (1) )Z− µ (A ν − Zµ (A ν
ωνα=0 + ωµα=0
ν,µ
− × Z− ν pq Zµ rs hpr|qsi
,
(46)
where
α=0 Z+ + Xα=0 ν = Yν ν
,
(47)
α=0 Z− − Xα=0 ν = Yν ν
,
(48)
and {ωνα=0 } are positive zeroth-order eigenvalues. The latter correspond to excitation energies for a group-uncorrelated, α = 0 limit, system. By confronting the block-diagonal structure of the main ERPA matrix eucidated in Eqs.(35)-(38) one notices that each of the eigenvalues is associated with one of the following types of excitations: active-active, activeinactive, secondary-active, and secondary-inactive ones. It is interesting to observe that even if a system is strongly correlated and some active orbitals become degenerate, which could result in zero-valued excitation energies in the active-active block, denominators in the expression for W (1) will in general stay different from zero. It is because ν and µ cannot both correspond to active-active excitations. Such terms are excluded due to a prime sign that indicates exclusion from the summation terms for which indices pqrs belong to a set of active orbitals. In other words, zero-valued excitation energies ωνα=0 do not lead to divergence in W (1) if ν’s correspond to the active-active block. Finally, the AC correlation energy in the linearized-AC-integrand approximation, named AC0, reads
ˆ AC0 Ecorr
1
(W α=0 + W (1) α) dα =
= 0
W (1) 2
,
(49)
where it has been used that W α=0 = 0. By inspection of W (1) it can be shown that the AC0 AC0 correlation correction Ecorr reduces to the well known MP2 second-order energy if the set of
13
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
active orbitals is empty, i.e. if the reference wavefunction is a single-determinant.
3
Implementation of AC-CAS and AC0-CAS
The exact formulation of the adiabatic connection correlation energy expression is given by AC Eqs.(16)-(18). It has led to obtaining two approximations for the correlation energy: the Ecorr AC0 and its less computationally demanding variant Ecorr , defined by Eqs.(32), (33) and Eqs.(46),
(49), respectively. Both of them rely on the extended random phase approximation and they involve solving either full ERPA problem (AC) or only smaller-sized ERPA equations for specific blocks, cf. Eqs.(35)-(38), at the coupling parameter α = 0 (AC0). Exploiting ERPA approximation implies that computation of the AC or AC0 correlation energy does not require construction of higher than two-electron reduced density matrices from a CAS wavefunction. This is a great advantage over the most often used implementations of the perturbation corrections to CAS which suffer from high computational cost due to the need of obtaining three- and four-electron reduced density matrices. It is worth noticing that if the active set is empty, then the reference wavefunction is a single determinant, the ERPA problem is equivalent to solving the time-dependent Hartree-Fock (TD-HF) equations, and the AC correlation correction is identical to one of the variants of the RPA correlation energy known as RPAx. 21 After taking into account the fact the matrices A(α) + B (α) and A(α) − B (α) are symmetric and positive-definite, if constructed from the optimized orbitals, it is convenient to turn the ERPA equation shown in Eq.(22) into a symmetric diagonalization problem the dimensionality of which is reduced by a factor of 2 with respect to the original ERPA equations. The matrices Aα and B α corresponding to a given CAS wavefunction are related, cf. the excitation operator Eq.(21), to excitations between inactive-active, active-active, inactive-secondary, and active-secondary orbitals so the dimension of the A/B matrix reads DimA = Minact Mact + Minact Msec + Mact Msec + Mact (Mact − 1)/2. The cost of constructing
14
ACS Paragon Plus Environment
Page 14 of 34
Page 15 of 34 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
2 the matrices Aα and B α scales as (DimA)2 Mact , whereas that of the diagonalization of the
ERPA matrix amounts to (DimA)3 . Thus, the overall scaling with respect to the number of 3 3 6 leading Msec . Notice that (DimA)3 involves a term reading Minact the active orbitals is Mact
to a similar overall scaling with respect to inactive and secondary orbitals as that of the RPA methods. To circumvent such high scaling, techniques like density fitting successfully employed in RPA could be also developed for AC. AC0 reduces the computational burden of AC, since finding the zeroth-order vectors , that enter Eq.(46) requires solving a set of symmetric eigenproblems for a sequence of Xα=0 ν small matrices A shown in Eqs.(35)-(38). For the largest matrix - the active block A = Aact 6 . Computation of the first-order term W (1) given in Eq.(46) - the scaling amounts to Mact 5 4 2 . and Msec Mact Mact involves steps the most expensive of which scale only as Msec
In Ref. 1 it has been anticipated that the AC integrand, Eq.(33), will show close to linear behavior for CAS wavefunction when α varies between 0 and 1 even for strongly correlated systems. Indeed, if the CAS wavefunction is employed, as in this work, the AC integrand stays nearly linear. However, for some systems we have observed that an instability may develop in the ERPA equations. When instability was occurring one of the eigenvalues was turning from a positive real number to an imaginary number. If α approaches the instability point such eigenvalue decreases toward zero and the AC integrand deviates from a linear behavior. While such instability is not possible to predict pre-calculation, it is easy to spot once the computation is done. In Fig.(1) and Fig.(6) we present illustrative behavior of the AC integrand for water and the nitrogen molecules in such a situation. When bonds dissociate the AC integrand changes abruptly close to α = 1 for H2 O (stretching of both OH bonds is considered) and N2 molecules (for N2 the instability shows up only if the CAS(10,10) reference is used). Since it only affects the AC integrand when α approaches 1, it has no impact on the potential energy curves, as it can be seen in Fig.(2) and Fig.(5). The AC0 approximation is also unaffected by an instability since it is based on expansion of W α around α = 0. On the other hand, the AC1 approximation proposed in Ref. 1 (shown to be equivalent
15
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
to the earlier introduced ERPA correlation correction 31,32 for strongly orthogonal geminal wavefunctions, and amounting to replacing W α with the linear function W α=1 α resulting in ´1 α=1 AC1 = 0 W α=1 α dα = W 2 ) should be the AC1 expression for the correlation energy: Ecorr used with caution with the CAS reference. Occurrence of instability would strongly affect the integrand at α = 1 and the AC1 correlation energy would show erratic behavior. Since for CAS reference AC0 is a more attractive variant of the adiabatic connection approach than the AC1 approximation the results presented in the next section focus on performance of AC and AC0 only.
4
Results
To assess the performance of the presented methods, let us look at several problems of multireference character. In order to evaluate the differences between the AC- and the perturbation-based approaches, we compare the AC-CAS and AC0-CAS results, obtained by adding, respectively, the AC and AC0 correction to the CAS energy, to two variants of NEVPT2 method - the partially-contracted (PC-) and the strongly-contracted (SC-) one 4 - employing the same CAS wavefunction. CAS wavefunctions used in AC-CAS and AC0CAS computations as well as all NEVPT2 results were obtained using DALTON software package. 33 MP2 natural orbitals were utilized as guess orbitals in CASSCF calculations carried out without imposing a system symmetry. For all systems a grid comprising 10 points was used for AC integration, except for 2, 6-pyridine, where only a 3-point quadrature has been employed. In fact, we have confirmed using the example of water molecule that a small grid with only 3 points is sufficient to achieve accuracy of 10−4 Hartree in the whole range of considered interatomic distances (also in the dissociation limit, where the instability on the AC integrand close to α = 1 occurs). To solve the problem of symmetric stretching of O-H bonds in a water molecule, an active space (4,4) [where the notation (e, o) stands for the number of active electrons e and
16
ACS Paragon Plus Environment
Page 16 of 34
Page 17 of 34 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
orbitals o] is required at minimum. As already noted, in Fig.(1) one can see, that in the equilibrium geometry the dependence of the AC integrand W α on the integration variable α is linear regardless of the CAS used. On the other hand, close to the dissociation limit (R = 7.0 a.u.) the curves abruptly deviate from linearity while α is nearing the value of 1. Despite this instability, the dissociation curve remains smooth (see Fig.(2)). In the case of the smallest employed CAS (4,4) the errors in the absolute energy are still fairly large for all methods varying from 12 to 22 mHartree, but the error in the dissociation energy exceeds one mHartree only for the AC0-CAS method (cf. Table 1). As the active space expands, the errors in the total energy diminish systematically for both AC methods. This is not the case, however, for the NEVPT2 approaches, as for the active space (8,8) the energy values are higher that these obtained for the smaller (6,10) space, in fact they are almost the same as for a small (4,4) space. When it comes to the dissociation energy: while the energy difference obtained by AC0-CAS is significantly improved by the expansion of CAS and the AC-CAS energy stays accurate, the PC-NEVPT2’s accuracy decreases upon the expansion, as it can be seen by comparing dissociating energies obtained in the (4,4) and (8,8) active spaces. Table 1: Total energies of H2 O molecule at equilibrium geometry and dissociation energies De computed in cc-pVDZ basis set. 34 Errors of De with respect to the FCI value are given in parentheses. CAS (4,4)
Energy at R= 1.8 a.u. (Hartree) AC-CAS AC0-CAS PC-NEVPT2 -76.22294 -76.22848 -76.22760
SC-NEVPT2 -76.22728
(6,10)
-76.23294
-76.23170
-76.23442
-76.23327
(8,8)
-76.23317
-76.23269
-76.22787
-76.22667
(8,12)
-76.23879
-76.23724
-76.23707
-76.23619
FCI
-76.24168
De (mHartree) AC-CAS AC0-CAS 332.1 336.1 (-0.8) (3.2) 332.0 332.6 (-0.9) (-0.2) 332.0 333.2 (-0.9) (0.3) 332.6 331.7 (-0.2) (-1.2)
PC-NEVPT2 332.4 (-0.4) 334.4 (1.5) 330.4 (-2.5) 331.8 (-1.1)
SC-NEVPT2 332.1 (-0.8) 334.5 (1.7) 330.6 (-2.3) 332.1 (-0.8)
To examine the dependence of the results on the chosen reference wavefunction in a more systematic manner, we chose the boron hydride molecule. With the number of active electrons in CAS fixed at 6 we were expanding the number of active orbitals from 6 to 16 17
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
0.00
(Hartree)
-0.05
-0.10
-0.15
-0.20
W
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
CAS(4,4) R=1.8 a.u.
-0.25
CAS(4,4) R=7.0 a.u. CAS(8,12) R=1.8 a.u.
-0.30
CAS(8,12) R=7.0 a.u.
0.0
0.2
0.4
0.6
0.8
1.0
Figure 1: AC integrand shown in Eq.(33) for the H2 O molecule at the equilibrium (R= 1.8 a.u.) geometry and near the dissociation limit (R= 7.0 a.u.) when both OH bonds are equally stretched. [CAS(6,19) would correspond to FCI in the cc-pVDZ basis set employed here]. As expected, the errors of all methods diminish with expansion of CAS (see Fig.(3) and Table 2). Note that the energies obtained with the NEVPT2 methods are higher than the ones from ACCAS and AC0-CAS irrespective of the CAS space. The dissociation energies provided by the AC-CAS and AC0-CAS approaches improve upon the expansion of the active space. In contrast, the NEVPT2’s behavior is less predictable, with the SC variant being especially erratic. Interestingly, for different active spaces and along the dissociation curves the total AC0-CAS energy parallels that of the PC-NEVPT2 one, but the AC0-CAS recovers more correlation energy for short BH distances, which results in better behavior of the dissociation energy. The problem of dissociation of a fluorine molecule requires an accurate description of dynamic correlation. Since even many-body perturbation methods with large active spaces often fail to correctly describe the dissociation curve, 36 it is a good case to test AC-CAS’s and AC0-CAS’s ability to capture the non-dynamic and dynamic correlation. Here, again
18
ACS Paragon Plus Environment
Page 18 of 34
Page 19 of 34
R (a.u.) 2
3
4
R (a.u.) 5
6
7
2
3
4
6
7 22
20
20
CAS=(8,8)
18
Error (mHartree)
5
22
18
16
16
14
14
12
12
10
10
8
8
6
6
CAS=(4,4)
4
4
2
2
22
22
20
20
CAS=(8,12)
CAS=(6,10)
18
18
16
16
AC-CAS
14
14
AC0-CAS
12
12
PC-NEVPT2 10
10
8
8
6
6
4
4
2
2 2
3
4
5
6
7
2
R (a.u.)
3
4
5
6
7
R (a.u.)
Figure 2: Errors of the total energy of H2 O molecule computed w.r.t. a FCI reference along the OH bond length in symmetric dissociation process. AC-CAS, AC0-CAS and PC-NEVPT2 energies were computed with a (4,4), (8,8), (6,10) and (8,12) CAS reference wavefunction in cc-pVDZ basis set.
19
ACS Paragon Plus Environment
Error (mHartree)
Error (mHartree)
Error (mHartree)
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
Journal of Chemical Theory and Computation
R (a.u.)
R (a.u.) 2
4
6
8
10
2
4
6
8
12
Error (mHartree)
10 14
14
CAS=(6,8)
CAS=(6,6)
12
10
10
8
8
6 6 4
Error (mHartree)
4 2 2 0
14
14
12
12
CAS=(6,12)
Error (mHartree)
CAS=(6,10) 10
10
8
8
6
AC-CAS AC0-CAS
4
6
4
PC-NEVPT2 2
2
0
0 2
4
6
8
10
2
R (a.u.)
4
6
8
10
R (a.u.)
Figure 3: Errors of the total energy of BH molecule computed w.r.t. a FCI reference along the BH bond length. AC-CAS, AC0-CAS and PC-NEVPT2 energies were computed with a (6,6), (6,8), (6,10) and (6,12) CAS reference wavefunction in cc-pVDZ basis set.
20
ACS Paragon Plus Environment
Error (mHartree)
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 20 of 34
Page 21 of 34
45
45 AC-CAS
Error (mHartree)
AC0-CAS 40
40
PC-NEVPT2
35
35
30
30
25
25
CAS=(10,10)
CAS=(2,2)
20
3
4
5
6
7
8
3
4
5
6
7
Error (mHartree)
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
20
8
R (a.u.)
R (a.u.)
Figure 4: Errors of the total energy of F2 molecule computed w.r.t. a CCSDTQ 35 reference along the F2 bond length. AC-CAS, AC0-CAS and PC-NEVPT2 energies were computed with a (2,2) and (10,10) CAS reference wavefunction in cc-pVDZ basis set.
21
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
Page 22 of 34
Table 2: Total energies of BH molecule at equilibrium geometry and dissociation energies De computed in cc-pVDZ basis set. 34 Errors of De with respect to the FCI value are given in parentheses. CAS (6,6)
Energy at R=2.3 a.u. (Hartree) AC-CAS AC0-CAS PC-NEVPT2 -25.20924 -25.20618 -25.20350
SC-NEVPT2 -25.19889
(6,8)
-25.21246
-25.21058
-25.20897
-25.20635
(6,10)
-25.21469
-25.21390
-25.21396
-25.21232
(6,12)
-25.21576
-25.21541
-25.21530
-25.21415
(6,14)
-25.21607
-25.21586
-25.21585
-25.21539
(6,16)
-25.21606
-25.21603
-25.21602
-25.21560
FCI
-25.21622
De (mHartree) AC-CAS AC0-CAS 120.5 122.9 (-5.8) (-3.4) 125.0 124.8 (-1.3) (-1.5) 125.3 125.4 (-1.0) (-0.9) 125.8 125.5 (-0.5) (-0.8) 126.8 126.3 (0.5) (0.0) 126.1 126.1 (-0.2) (-0.2)
PC-NEVPT2 122.0 (-4.3) 124.0 (-2.3) 125.9 (-0.4) 125.6 (-0.7) 126.1 (-0.2) 126.1 (-0.2)
SC-NEVPT2 120.9 (-5.4) 123.7 (-2.6) 126.1 (-0.2) 124.6 (-1.7) 125.7 (-0.6) 125.7 (-0.6)
(see Fig.4 and Table 3) all the methods predict correctly the dissociation energy. For the smaller active space, (2,2), the NEVPT2 and AC0-CAS energies are close to the benchmark energy, while the AC-CAS energy is much higher, revealing that the AC correction misses a large portion of dynamic correlation. This is related to the slight deviation of the W α integrand from linearity. Since the W α curve has a larger slope close to zero then near one, AC0-CAS recovers more correlation then AC-CAS, which in this case favors AC0. The corresponding AC-CAS dissociation energy is, however on par with the ones produced by other methods. In the larger, (10,10), active space all the methods produce similar total energies, with errors varying from 24 to 28 mHartree, and accurate dissociation energies. Again, the AC0-CAS total energies follow the NEVPT2 results rather than AC-CAS and both AC-CAS and AC0-CAS (and in this case also NEVPT2) dissociation energies improve upon the expansion of the active space. The dissociation of the nitrogen molecule is another multireference problem. It involves 6 electrons forming the triple bond and therefore the minimal CAS able to describe this process is (6,6). Indeed, from Table 4 it is clear the CAS(6,6) is sufficient to describe
22
ACS Paragon Plus Environment
Page 23 of 34 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
Table 3: Total energies of F2 molecule at equilibrium geometry and dissociation energies De computed in cc-pVDZ basis set. 34 Errors of De with respect to the CCSDTQ value 35 are given in parentheses. CAS (2,2)
Energy at R=2.8 a.u. (Hartree) AC-CAS AC0-CAS PC-NEVPT2 -199.05814 -199.08206 -199.08397
SC-NEVPT2 -199.08392
(10,10)
-199.07620
-199.07718
CCSDTQ 35
-199.07658
-199.07905
De (mHartree) AC-CAS AC0-CAS 41.5 48.2 (-3.5) (3.2) 43.9 44.5 (-1.1) (-0.5)
PC-NEVPT2 47.6 (2.6) 46.0 (1.0)
SC-NEVPT2 47.5 (2.5) 44.7 (-0.3)
-199.10390
the dissociation – the errors of the dissociation energy obtained by NEVPT2 methods are small, 0.3 and 1.4 mHartree for PC and SC variant, respectively. The AC-CAS and AC0CAS methods perform much worse, although the errors in the absolute energies are not dramatically higher. As CAS is expanded, the performance of AC0-CAS and AC-CAS improves both in terms of the dissociation energy and the absolute energy. This is not the case for NEVPT2 methods whose behavior is quite erratic upon the expansion – the errors in the dissociation energy grow and even the absolute PC-NEVPT2 energy for CAS(10,10) is for most geometries higher than for CAS(6,6) and CAS(6,12). This is a serious drawback of NEVPT2, as it prevents one from using the method in a “black box” manner. As already observed, the good behavior of AC-CAS is not hindered by the instability of the AC integrand (see Fig.(6)) appearing for the stretched structure when CAS(10,10) is employed. Table 4: Total energies of N2 molecule at equilibrium geometry and dissociation energies De computed in cc-pVDZ basis set. 34 Errors of De with respect to the MRCISDTQ value 37 are given in parentheses. CAS (6,6)
Energy at R=2.1 a.u. (Hartree) AC-CAS AC0-CAS PC-NEVPT2 -109.25025 -109.24645 -109.24849
SC-NEVPT2 -109.24734
(6,12)
-109.25405
-109.24677
-109.25153
-109.24796
(10,10)
-109.26208
-109.26140
-109.26043
-109.25794
MRCISDTQ 37
De (mHartree) AC-CAS AC0-CAS 325.6 327.7 (4.6) (6.7) 324.3 324.3 (3.3) (3.3) 321.7 325.8 (0.7) (4.8)
PC-NEVPT2 320.7 (-0.3) 316.4 (-4.6) 334.0 (13.0)
SC-NEVPT2 319.6 (-1.4) 314.1 (-6.9) 332.3 (11.3)
-109.27804
Finally, let us look at two forms of a 2,6-pyridine molecule. Coupled-cluster methods 38,39 23
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
Error (mHartree)
40
40
CAS=(6,12)
CAS=(10,10)
35
35
30
30
25
25
AC-CAS 20
20
AC0-CAS
CAS=(6,6)
PC-NEVPT2
15
15 2
4
6
8
10 2
4
R (a.u.)
6
8
10 2
R (a.u.)
4
6
8
R (a.u.)
Figure 5: Errors of the total energy of N2 molecule computed w.r.t. a MRCISDTQ 37 reference vs. the N2 bond length. AC-CAS, AC0-CAS and PC-NEVPT2 energies were computed with a (6,6) and (6,12) and (10,10) CAS reference wavefunction in cc-pVDZ 34 basis set. 0.00 CAS(6,6) R=2.1 a.u. CAS(6,6) R=10.0 a.u.
(Hartree)
-0.05
CAS(10,10) R=2.1 a.u. CAS(10,10) R=10.0 a.u.
-0.10
-0.15
W
-0.20
-0.25
-0.30
0.0
0.2
0.4
0.6
0.8
1.0
Figure 6: AC integrand for the N2 molecule at the equilibrium (R= 2.1 a.u.) geometry and near the dissociation limit (R= 10.0 a.u.). predict two stable structures of this molecule: monocyclic and bicyclic, but they disagree as to which of those is the global minimum. According to the state-of-art Mk-MRCCSDT 39 computations, the monocyclic structure is lower than the bicyclic by 5.7 kcal/mol, 39 while 24
ACS Paragon Plus Environment
10
Error (mHartree)
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 24 of 34
Page 25 of 34 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
Mk-MRCCSD(T) predicts a gap of 8.8 kcal/mol. 38 It has even been suggested that the bicyclic form might be an artifact of the theory. 39 While it seems impossible to answer that question without going to very large active spaces, pyridine is a good case to present the efficiency of both AC-CAS and the AC0-CAS method, as in the employed cc-pVTZ basis set 34 pyridine is described by 222 basis functions and we were able to use CAS as large as (14,14). The smallest CAS used was (8,8), suggested in, 40 while the largest was (14,14). CASSCF(14,14) clearly cannot describe the problem accurately without any correlation corrections, as the relative energy, ∆E of the bicyclic structure obtained with this method, 12.8 kcal/mol, is still over twice the energy obtained by Mk-MRCCSDT. Looking at the performance of the correlation corrections one observes that for the CAS(8,8) reference ACCAS produces accurate relative energy (see Table 5), differing by only 0.1 kcal/mol from the coupled cluster reference. The other methods - AC0-CAS and NEVPT2 yield similar ∆E of 7.3 − 7.4 kcal/mol. As a result of expanding the CAS space to (10,10), the ∆E obtained with NEVPT2 methods grows to 11.6 kcal/mol and 10.8 kcal/mol for PC-NEVPT2 and SCNEVPT2, respectively and to 9.1 kcal/mol and 10.4 kcal/mol for AC-CAS and AC0-CAS, respectively. Counterintuitively, the total NEVPT2 and AC0-CAS energies are higher for this larger space than for CAS(8,8). This deterioration of the performance of all the methods upon the expansion of CAS indicates that both (8,8) and (10,10) spaces are insufficient to capture the multireference character of the system. Namely, upon expanding CAS from (10,10) to (12,12) the composition of the wavefunction representing the monocyclic structure evolves significantly, with the dominant configuration remaining the same and changing contributions from configurations entering with coefficients of order of 0.1. In spaces (12,12) and (14,14) AC-CAS already predicts an accurate value of the relative energy of 5.0 kcal/mol. Similarly, the result produced by AC0-CAS remains the same for both those two spaces. In the future those conclusions should be verified by repeating the calculations in larger spaces. Unfortunately, the discrepancy between the results from Mk-MRCCSDT and MkMRCCSD(T) methods suggests that even our best benchmark might not be fully reliable.
25
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
Page 26 of 34
Nevertheless, AC-CAS and AC0-CAS have once again shown stability in large active spaces. In this case AC0-CAS was able to showcase its computational efficiency, beating for large active space the NEVPT2 method by an order of magnitude (even when confronting its plain implementation with an optimal one of NEVPT2 in Dalton) and AC-CAS by two orders in computational time (see Table 6). Table 5: Total energies of the monocyclic structure of pyridine and the relative energies of the bicyclic structure computed in cc-pVTZ basis set. 34 Geometries taken from Ref. 38. Errors with respect to a Mk-MRCCSDT/cc-pVTZ value 39 are given in parentheses. CAS (8,8)
Energy of the monocyclic structure (Hartree) AC-CAS AC0-CAS PC-NEVPT2 SC-NEVPT2 -246.43118 -246.48159 -246.48864 -246.48578
(10,10)
-246.44008
-246.46954
(12,12)
-246.44509
-246.46233
–
–
(14,14)
-246.45107
-246.46308
–
–
-246.47995
-246.47415
Mk-MRCCSDT 39
∆E (kcal/mol) AC-CAS AC0-CAS 5.6 7.3 (-0.1) (1.6) 9.1 10.4 (3.4) (4.7) 5.0 4.2 (-0.7) (-1.5) 5.0 4.3 (-0.7) (-1.4) 5.7
PC-NEVPT2 7.3 (1.6) 11.6 (5.9) – –
SC-NEVPT2 7.4 (1.7) 10.8 (5.1) – –
Table 6: Timings of the computation of the correlation corrections for a single structure of pyridine in seconds on a single CPU. NEVPT2 correction was computed in Dalton software package and AC-CAS and AC0-CAS in our in-house code.1 Estimate for single grid point. CAS (8,8) (10,10) (12,12) (14,14)
5
AC-CAS1 1 · 104 1 · 104 1 · 104 2 · 104
AC0-CAS 5 · 102 5 · 102 6 · 102 8 · 102
NEVPT2 2 · 102 3 · 103 -
Conclusions
We have used the idea presented in Ref. 1 of employing the adiabatic connection formalism to account for dynamic correlation energy for a multireference wavefunction to correct the complete active space model. Two variants have been implemented: the AC method based on Eqs.(32) and (33) and the AC0 approximation given in Eqs.(46), (49). Both methods 26
ACS Paragon Plus Environment
Page 27 of 34 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
are based upon the extended random phase approximation and consequently they are sizeconsistent. In the limit of no active orbitals the AC approach is equivalent to the RPAx variant, whereas the AC0 approximation turns into the well-known MP2 correction. When all orbitals belong to the active space (the full CI limit) both AC and AC0 corrections vanish by construction, as it should be. Using examples of molecules requiring multireference treatment we have shown that both AC and AC0 offer accuracy similar to or exceeding that of the NEVPT2 method for absolute values of the energy in the minimal active spaces. One would expect that employing a more and more elaborate CAS reference augmented with the correlation correction will lead to more accurate predictions. As it has been shown the AC correction meets this expectation unlike, in general, the NEVPT2 one. We have found that when spaces of active electrons and orbitals are expanded, errors for both energy values and energy barriers predicted by AC-CAS (note the exception of the pyridine molecule) systematically decrease. This desirable behavior has not been generally observed when the NEVPT2 method (in particular its strongly-contracted variant) was employed. Both AC and AC0 methods scale only with the 6th power with respect to the number of the active orbitals. This is a serious advantage over the perturbation corrections and it makes the adiabatic correction suited for treating large active spaces (with possible application in DMRG). The most demanding step of the AC approach is solving the ERPA equations. This is alleviated in the AC0 approximation, as only blocks of the full ERPA matrix need to be diagonalized. Overall the AC0 approach appears to be the most computationally efficient method for including the dynamic correlation energy in the CAS model. Taking into account that the extension of the AC and AC0 to open-shell systems is straightforward and the extension to excited states possible (work in progress) the ACderived approximations emerge as attractive alternatives to perturbation methods for CAS reference models.
27
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
6
Page 28 of 34
Acknowledgements
This work was supported by the National Science Center of Poland under grant no. 2016/23/B/ST4/02848.
References (1) Pernal, K. Electron Correlation from the Adiabatic Connection for Multireference Wave Functions. Phys. Rev. Lett. 2018, 120, 013001. (2) Andersson, K.; Malmqvist, P. A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. Secondorder perturbation theory with a CASSCF reference function. J. Phys. Chem. 1990, 94, 5483–5488. (3) 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. (4) 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. (5) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. n-electron valence state perturbation theory: A spinless formulation and an efficient implementation of the strongly contracted and of the partially contracted variants. J. Chem. Phys. 2002, 117, 9138. (6) Evangelisti, S.; Daudey, J.; Malrieu, J. Qualitative intruder-state problems in effective Hamiltonian theory and their solution through intermediate Hamiltonians. Phys. Rev. A 1987, 35, 4930. (7) 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. 28
ACS Paragon Plus Environment
Page 29 of 34 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
(8) Roos, B. O.; Andersson, K. Multiconfigurational perturbation theory with level shift ? the Cr2 potential revisited. Chem. Phys. Lett. 1995, 245, 215 – 223. (9) 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. (10) Zgid, D.; Ghosh, D.; Neuscamman, E.; Chan, G. K.-L. A study of cumulant approximations to n-electron valence multireference perturbation theory. J. Chem. Phys. 2009, 130, 194107. (11) Sokolov, A. Y.; Chan, G. K.-L. A time-dependent formulation of multi-reference perturbation theory. J. Chem. Phys. 2016, 144, 064102. (12) Sokolov, A. Y.; Guo, S.; Ronca, E.; Chan, G. K.-L. Time-dependent N-electron valence perturbation theory with matrix product state reference wavefunctions for large active spaces and basis sets: Applications to the chromium dimer and all-trans polyenes. J. Chem. Phys. 2017, 146, 244102. (13) Li, C.; Evangelista, F. A. Multireference Driven Similarity Renormalization Group: A Second-Order Perturbative Analysis. J. Chem. Theo. Comp. 2015, 11, 2097–2108. (14) Lyakh, D. I.; Musial, M.; Lotrich, V. F.; Bartlett, R. J. Multireference Nature of Chemistry: The Coupled-Cluster View. Chem. Rev. 2012, 112, 182–243. (15) Kohn, A.; Hanauer, M.; Muck, L. A.; Jagau, T.-C.; Gauss, J. State-specific multireference coupled-cluster theory. WIREs Comput. Mol. Sci. 2013, 3, 176–197. (16) Fromager, E.; Toulouse, J.; Jensen, H. J. A. On the universality of the long-/short-range separation in multiconfigurational density-functional theory. J. Chem. Phys. 2007, 126, 074111.
29
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
(17) Li Manni, G.; Carlson, R. K.; Luo, S.; Ma, D.; Olsen, J.; Truhlar, D. G.; Gagliardi, L. Multiconfiguration Pair-Density Functional Theory. J. Chem. Theo. Comp. 2014, 10, 3669–3680. (18) McWeeny, R. The density matrix in many-electron quantum mechanics I. Generalized product functions. Factorization and physical interpretation of the density matrices. Proc. Royal Soc. A 1959, 253, 242–259. (19) Langreth, D.; Perdew, J. Exchange-correlation energy of a metallic surface: Wavevector analysis. Phys. Rev. B 1977, 15, 2884. (20) Eshuis, H.; Bates, J.; Furche, F. Electron correlation methods based on the random phase approximation. Theor. Chem. Acc. 2012, 131, 1084. ´ (21) Toulouse, J.; Zhu, W.; Angy´ an, J. G.; Savin, A. Range-separated density-functional theory with the random-phase approximation: Detailed formalism and illustrative applications. Phys. Rev. A 2010, 82, 032502. (22) Ren, X.; Rinke, P.; Joas, C.; Scheffler, M. Random-phase approximation and its applications in computational chemistry and materials science. J. Mater. Sci. 2012, 47, 7447. (23) Chatterjee, K.; Pernal, K. Excitation energies from extended random phase approximation employed with approximate one-and two-electron reduced density matrices. J. Chem. Phys. 2012, 137, 204109. (24) McLachlan, A. D.; Ball, M. A. Time-Dependent Hartree-Fock Theory for Molecules. Rev. Mod. Phys. 1964, 36, 844. (25) Pernal, K. Correlation energy from random phase approximations: A reduced density matrices perspective. Int. J. Quantum Chem. 2018, 118, e25462–n/a, e25462.
30
ACS Paragon Plus Environment
Page 30 of 34
Page 31 of 34 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
(26) Ramos-Cordoba, E.; Salvador, P.; Matito, E. Separation of dynamic and nondynamic correlation. Phys. Chem. Chem. Phys. 2016, 18, 24015. (27) Pernal, K.; Chatterjee, K.; Kowalski, P. H. How accurate is the strongly orthogonal geminal theory in predicting excitation energies? Comparison of the extended random phase approximation and the linear response theory approaches. J. Chem. Phys. 2014, 140, 014101. (28) Rowe, D. J. Equations-of-motion method and the extended shell model. Rev. Mod. Phys. 1968, 40, 153. (29) Pernal, K. Intergeminal Correction to the Antisymmetrized Product of Strongly Orthogonal Geminals Derived from the Extended Random Phase Approximation. J. Chem. Theo. Comp. 2014, 10, 4332–4341. (30) Notice a sign error in Eq.(44) and Eq.(15) in Refs. 25 and 1, respectively. (31) Pastorczak, E.; Pernal, K. ERPA-APSG: a computationally efficient geminal-based method for accurate description of chemical systems. Phys. Chem. Chem. Phys. 2015, 17, 8622. (32) Chatterjee, K.; Pastorczak, E.; Jawulski, K.; Pernal, K. A minimalistic approach to static and dynamic electron correlations: Amending generalized valence bond method with extended random phase approximation correlation correction. J. Chem. Phys. 2016, 144, 244111. (33) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P.; Dalskov, E. K.; Ekstr¨om, U.; Enevoldsen, T.; Eriksen, J. J.; Ettenhuber, P.; Fern´andez, B.; Ferrighi, L.; Fliegl, H.; Frediani, L.; Hald, K.; Halkier, A.; H¨attig, C.; Heiberg, H.; Helgaker, T.; Hennum, A. C.; Hettema, H.; Hjertenæs, E.; Høst, S.; Høyvik, I.-M.; Iozzi, M. F.; Jans´ık, B.; Jensen, H.
31
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
J. Aa.; Jonsson, D.; Jørgensen, P.; Kauczor, J.; Kirpekar, S.; Kjærgaard, T.; Klopper, W.; Knecht, S.; Kobayashi, R.; Koch, H.; Kongsted, J.; Krapp, A.; Kristensen, K.; Ligabue, A.; Lutnæs, O. B.; Melo, J. I.; Mikkelsen, K. V.; Myhre, R. H.; Neiss, C.; Nielsen, C. B.; Norman, P.; Olsen, J.; Olsen, J. M. H.; Osted, A.; Packer, M. J.; Pawlowski, F.; Pedersen, T. B.; Provasi, P. F.; Reine, S.; Rinkevicius, Z.; Ruden, T. A.; Ruud, K.; Rybkin, V. V.; Salek, P.; Samson, C. C. M.; de Mer´as, A. S.; Saue, T.; Sauer, S. P. A.; Schimmelpfennig, B.; Sneskov, K.; Steindal, A. H.; SylvesterHvid, K. O.; Taylor, P. R.; Teale, A. M.; Tellgren, E. I.; Tew, D. P.; Thorvaldsen, A. J.; Thøgersen, L.; Vahtras, O.; Watson, M. A.; Wilson, D. J. D.; Ziolkowski, M.; ˚ Agren, H. The Dalton quantum chemistry program system. WIREs Comput. Mol. Sci. 2014, 4, 269–284. (34) Dunning Jr., T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007. (35) Musial, M.; Bartlett, R. J. Multi-reference Fock space coupled-cluster method in the intermediate Hamiltonian formulation for potential energy surfaces. J. Chem. Phys. 2011, 135, 044121. ¨ Reiher, M. Orbital Entanglement (36) Boguslawski, K.; Tecmer, P.; Barcza, G.; Legeza, O.; in Bond-Formation Processes. J. Chem. Theo. Comp. 2013, 9, 2959–2973. (37) Hanauer, M.; Kohn, A. Pilot applications of internally contracted multireference coupled cluster theory, and how to choose the cluster operator properly. J. Chem. Phys. 2011, 134, 204111. (38) Evangelista, F. A.; Prochnow, E.; Gauss, J.; Schaefer III, H. F. Perturbative triples corrections in state-specific multireference coupled cluster theory. J. Chem. Phys. 2010, 132, 074107.
32
ACS Paragon Plus Environment
Page 32 of 34
Page 33 of 34 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
(39) Prochnow, E.; Harding, M. E.; Gauss, J. Parallel calculation of CCSDT and MkMRCCSDT energies. J. Chem. Theo. Comp. 2010, 6, 2339–2347. (40) Cramer, C. J.; Debbert, S. Heteroatomic substitution in aromatic σ biradicals: the six pyridynes. Chem. Phys. Lett. 1998, 287, 320–326.
33
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
Graphical TOC Entry EAC corr =
34
Z
1
Wα dα 0
ACS Paragon Plus Environment
Page 34 of 34
