Table Of ContentBenchmark of dynamic electron correlation models for seniority-zero wavefunctions
and their application to thermochemistry
Katharina Boguslawski∗
Institute of Physics, Faculty of Physics, Astronomy and Informatics,
Nicolaus Copernicus University in Torun, Grudziadzka 5, 87-100 Torun, Poland and
Departement of Chemistry, Nicolaus Copernicus University in Torun, Gagarina 7, 87-100 Torun, Poland
Pawel(cid:32) Tecmer†
Institute of Physics, Faculty of Physics, Astronomy and Informatics,
Nicolaus Copernicus University in Torun, Grudziadzka 5, 87-100 Torun, Poland
7
1 Wavefunctions restricted to electron-pair states are promising models to describe
0 static/nondynamic electron correlation effects encountered, for instance, in bond-dissociation
2 processes and transition-metal and actinide chemistry. To reach spectroscopic accuracy, however,
the missing dynamic electron correlation effects that cannot be described by electron-pair states
n
need to be included a posteriori. In this article, we extend the previously presented perturbation
a
J theory models with an Antisymmetric Product of 1-reference orbital Geminal (AP1roG) reference
function that allow us to describe both static/nondynamic and dynamic electron correlation
7
effects. Specifically,ourperturbationtheorymodelscombineadiagonalandoff-diagonalzero-order
1
Hamiltonian, a single-reference and multi-reference dual state, and different excitation operators
] used to construct the projection manifold. We benchmark all proposed models as well as an a
h posteriori linearized coupled cluster correction on top of AP1roG against CR-CCSD(T) reference
p data for reaction energies of several closed-shell molecules that are extrapolated to the basis set
- limit. Moreover,wetesttheperformanceofournewmethodsformultiplebondbreakingprocesses
m
in the N , C , and BN dimers against MRCI-SD and MRCI-SD+Q reference data. Our numerical
2 2
e results indicate that the best performance is obtained from a linearized coupled cluster correction
h as well as second-order perturbation theory corrections employing a diagonal and off-diagonal
c zero-order Hamiltonian and a single-determinant dual state. These dynamic corrections on top of
.
s AP1roG allow us to reliably model molecular systems dominated by static/nondynamic as well as
c dynamic electron correlation.
i
s
y
h I. INTRODUCTION tum chemistry do not a priori account for the dynamic
p part of the correlation energy. Exceptions are multi-
[
Accurate and reliable theoretical predictions of ther- reference coupled cluster models, some variants of the
1 modynamic properties of chemical reactions remain an multi-reference configuration interaction method, and
v active field of research [1]. Computational studies are the perturb-then-diagonalize multi-reference perturba-
3 particularly desirable especially when experimental ma- tion theory approaches that use effective Hamiltonians.
6 nipulations are limited or impossible. Conventional To remedy this problem, a posteriori corrections have
5
“ab initio” quantum chemistry tools used in thermo- been developed that are usually based on perturbation
4
0 chemistry are commonly based on different flavours of theory [11–15], canonical transformation theory [16] or
. many-body perturbation theory and coupled cluster ap- coupled cluster methods [17].
1
proaches [2–9]. Over the years, the CCSD(T) (Cou-
0 Although multi-reference methods are frequently ap-
7 pled Cluster Singles, Doubles and perturbative Triples)
plied to model strongly-correlated systems, like bond-
1 method has unfolded as one of the most reliable and ac-
breaking processes and heavy-element chemistry, some
: curate models for systems that are well-represented by
v flavours of single-reference coupled cluster theory pose
a single electron configuration. Therefore, it is often
i
X referred to as the “gold standard” of quantum chem- promising alternatives to describe static/nondynamic
electron correlation. Among the most representa-
r istry. Although CCSD(T) performs extraordinary well
a tive examples are the CR-CC (Completely Renor-
forsystemsdominatedbydynamiccorrelation,themodel
malized Coupled Cluster) and CR-EOMCC (Com-
often fails for multi-reference systems. For molecu-
pletely Renormalized Equation of Motion Coupled
lar systems with a substantial multi-reference character,
Cluster) approaches [18–22], the active-space CC and
multi-reference methods [10] are commonly used, which
EOMCC methods [23–26], the EA/IP- (Electron Affin-
are computationally very expensive. Furthermore, the
ity/IonizationPotential)andDEA/DIP-(DoubleEA/IP)
most popular multi-reference approaches used in quan-
EOMCC models [27], and the spin-flip CC/EOMCC
formalism [28–30]. A different, computationally fea-
sible approach suitable for strongly-correlated sys-
∗ k.boguslawski@fizyka.umk.pl tems uses seniority-zero wavefunctions to describe the
† ptecmer@fizyka.umk.pl static/nondynamic part of the electron correlation en-
2
ergy. [31–44] The missing dynamic electron correlation tion, two different PT models have been proposed that
effects are included a posteriori in these ans¨atze using, allow us to describe electron correlation effects beyond
for instance, many-body perturbation theory [45–47], electron pairs. [46] While these approaches provide reli-
coupled-cluster theory [48–52], extended random phase able spectroscopic constants for some first-row diatomic
approximation[53],anddensityfunctionaltheory(DFT) molecules, [32, 51] their performance deteriorates when
corrections [54, 55]. moving to heavy-element containing compounds like ac-
In this article, we seek computationally cheap and tinide species. [51] To remedy this problem, a linearized
reliable dynamic electron correlation models on top of coupledcluster(LCC)correctionsontopofAP1roGcan
the Antisymmetric Product of 1-reference orbital Gemi- beusedthatallowsustoaccuratelydescribeclosed-shell
nal (AP1roG) of different levels of approximation. The moleculeswith(quasi-)degeneratef-,d-,p-,ands-shells,
AP1roG wavefunction can be written using an exponen- asencounteredinactinide-containingcompounds.[51]In
tial ansatz, the following, we will briefly summarize some theoretical
models for dynamic correlation with an AP1roG refer-
(cid:32)(cid:88)P (cid:88)K (cid:33) ence function. Specifically, we will extend the existing
|AP1roG(cid:105)=exp caia†aa†a¯a¯iai |0(cid:105), (1) PT models and relate the proposed models to the re-
i=1a=P+1 cently presented PTa and PTb methods. [46]
wherethesumrunsoverallelectronpairsP (equaltothe
number of occupiedorbitals) and virtual orbitals K−P,
A. 2nd-order Perturbation Theory
while|0(cid:105)issomereferencedeterminant. a† anda arethe
p p
fermioniccreationandannihilationoperatorsforspin-up
One drawback of multi-reference perturbation theory
p and spin-down p¯ electrons, while ca are the AP1roG
i is the arbitrariness of the theoretical model. For exam-
amplitudes, also called geminal coefficients. The gemi-
ple, there are different choices for the zero-order Hamil-
nal coefficients ca are obtained by solving the projected
i tonian Hˆ , the dual state (cid:104)Ψ˜| in the projector, and the
Schr¨odinger equation, as in coupled cluster theory. The 0
choice of the projection space. [14] Poor choices may
AP1roG model allows us to describe static/nondynamic
lead to technical difficulties and unphysical solutions. In
correlation and to a minor extent some dynamic correla-
this work, we will investigate different selections for the
tion(seealsoref.[56]). Toincludethemissingpartofthe
zero-order Hamiltonian Hˆ , the dual state (cid:104)Ψ˜|, and the
dynamic electron correlation effects a posteriori that go 0
projection space, while the zero-order wavefunction is
beyondtheAP1roGansatz,thatis,beyondelectronpair
restricted to the AP1roG reference function of eq. (1),
states, we will use perturbation theory as well as a lin-
|Ψ(0)(cid:105)=|AP1roG(cid:105).
earized coupled cluster corrections. [51] Specifically, our
Furthermore, in the derivation of all PT models, we
perturbation theory models will combine a diagonal and
willusethequantumchemicalHamiltonianinitsnormal
off-diagonal zero-order Hamiltonian, a single- and multi-
productfromshifted bythecorrelationenergyofAP1roG
determinant wavefunction as dual, and different excita-
tionoperatorsusedtoconstructtheprojectionmanifold. Ec(o0r)r (the shift in energy is indicated by “(cid:48)”),
The performance of all perturbation theory models as
Hˆ(cid:48) =Hˆ −(cid:104)0|Hˆ|0(cid:105)−E(0)
well as the linearized coupled cluster correction, as pre- N corr
sented in ref. 51, will be assessed against spectroscopic =(Hˆ −(cid:104)0|Hˆ |0(cid:105))+(Vˆ −(cid:104)0|Vˆ|0(cid:105)−E(0) )
0 0 corr
constants of homo- and heteronuclear dimers and reac-
=(Hˆ ) +Vˆ(cid:48) , (2)
tion energies of closed shell systems. 0 N N
This article is organized as follows. In section II,
where the quantum chemical Hamiltonian Hˆ is divided
we present different theoretical models to correct for
into a zero-order contribution Hˆ and a perturbation Vˆ.
the missing dynamic electron correlation in the AP1roG 0
It is convenient to rewrite Hˆ into a sum of a one- (Hˆ )
ansatz. Computational details are described in sec- 1
and a two-electron (Hˆ , here indicated by Wˆ) part,
tion III. Numerical results and comparison to standard 2
electroncorrelationmethodsformolecularsystemsdom-
inated by static/nondynamic (bond dissociation pro- Hˆ =Hˆ +Wˆ =(cid:88)h a†a +1 (cid:88)(cid:104)pq|rs(cid:105)a†a†a a . (3)
1 pq p q 2 p q s r
cesses) and dynamic electron correlation (reaction en- pq pqrs
ergies of main group compounds) are presented in sec-
In the above equation, h and (cid:104)pq|rs(cid:105) are the one- and
tion IV. Finally, we conclude in section V. pq
two-electron (written in physicists’ notation) integrals,
respectively, determined for the one-particle basis func-
tions p,q,r,s. We should note that, in this work, we will
II. THEORETICAL MODELS FOR DYNAMIC
CORRELATION restrictHˆ0 tobeaone-bodyoperatorsothatthenormal-
product form of the perturbation Vˆ = Hˆ −Hˆ can be
0
written as an operator shifted by the AP1roG correla-
A common way to account for dynamic correlation ef-
fects in quantum many-body systems is to use perturba- tionenergy,VˆN(cid:48) =VˆN−Ec(o0r)r (againindicatedby“(cid:48)”),as
tion theory (PT). Specifically for the AP1roG wavefunc- we have E(0) = (cid:104)0|Wˆ |AP1roG(cid:105). Introducing a shifted
corr N
3
perturbationoperatorVˆ(cid:48) isequivalenttoneglectingcon- conventional CC theory, but adjusted during the opti-
N
tractions (or diagrams) that correspond to the AP1roG mization of the AP1roG wavefunction. Choosing |0(cid:105) as
correlation energy in the PT equations, which will be reference determinant, the first-order correction can be
indicated by the (cid:48) in the sum of the Wˆ(cid:48) operator. written as
N
As in conventional Rayleigh–Schr¨odinger PT, the ex-
actwavefunctioncanbewrittenasanorder-by-orderex- |Ψ(1)(cid:105)=Tˆ|0(cid:105), (9)
pansion, |Ψ(cid:105)=|AP1roG(cid:105)+λ|Ψ(1)(cid:105)+..., where λ is the
order parameter. The first-order correction to the wave- whereTˆissomeexcitationoperatorthatsubstituteselec-
function is expanded in a set of Slater determinants Φ , tronsfromtheoccupiedtothevirtualspacewithrespect
p
to |0(cid:105). Furthermore, we will restrict Tˆ to contain dou-
|Ψ(1)(cid:105)=(cid:88)t |Φ (cid:105), (4) ble excitations without electron pairs, Tˆ = Tˆ(cid:48), as well
p p 2
p as single and double excitations, Tˆ = Tˆ1 +Tˆ2(cid:48). If only
double excitations are included, the excitation operator
and forced to be orthogonal to the zero-order wavefunc-
is specified as
tion, here, |AP1roG(cid:105),
occ virt
(cid:104)Ψ(1)|AP1roG(cid:105)=0. (5) Tˆ(cid:48) = 1(cid:88)(cid:88)(cid:48)tabEˆ Eˆ , (10)
2 2 ij ai bj
ij ab
This orthogonality constraint restricts the choice of the
projection space used for the expansion of |Ψ(1)(cid:105) (and where Eˆai = a†aai +a†a¯a¯i is the singlet excitation oper-
higherorders). Byconstruction,allpair-excitationswith ator and the perturbation amplitudes tab are symmetric
ij
respect to |0(cid:105) have to be excluded as they don’t satisfy
with respect to pair-exchange, i.e., tab = tba. Similar
(cid:104)Φ |AP1roG(cid:105) = 0. Similarly, introducing an order pa- ij ji
p to our notation of the shifted normal-product Hamilto-
rameterλintheHamiltonianHˆ =Hˆ +λVˆ andequating
0 nian, the prime in the above summation indicates that
coefficientsofpowersofλ,weobtainthezero-,first-,and
pair-excited determinants are excluded in the excitation
higher-orderPTequations. Specificallyforthefirst-order operator, i.e., taa¯ = 0. Exclusion of pair-excitations ful-
correction to the wavefunction, we have to solve i¯i
fils the orthogonality condition eq. (5). If the excitation
operatorcontainsbothsingleanddoubleexcitations,the
(Hˆ ) |Ψ(1)(cid:105)+Vˆ(cid:48) |AP1roG(cid:105)=0
0 N N single excitations can be accounted for by adding
(cid:88)t (Hˆ ) |Φ (cid:105)+Vˆ(cid:48) |AP1roG(cid:105)=0. (6)
p 0 N p N
occ virt
p Tˆ =(cid:88)(cid:88)taEˆ (11)
1 i ai
Since we have introduced a shifted normal-product i a
Hamiltonian, the zero- and first-order energy corrections
to the double excitation operator. Furthermore, as basis
vanish,
for the bra states of the projection manifold, we follow
1 refs. 46, 51 and use the convenient choice
E(0)+E(1) = (cid:104)Ψ˜|Vˆ(cid:48) |AP1roG(cid:105)=0, (7)
(cid:104)Ψ˜|AP1roG(cid:105) N 1 1
(cid:104)ab|= (cid:104)ab|+ (cid:104)ab|, (12)
ij 3 ij 6 ji
where(cid:104)Ψ˜|isthedualoftheunperturbedstate|AP1roG(cid:105).
Specificchoicesfor(cid:104)Ψ˜|willbeconsideredbelow. Thefirst for doubly excited determinants where (cid:104)ab| = (cid:104)0|Eˆ Eˆ
ij jb ia
non-zero correction to the energy is of second-order and and
can be calculated from the first-order wavefunction and
1
the shifted normal-product perturbation Hamiltonian, (cid:104)a|= (cid:104)a|, (13)
i 2 i
(cid:104)Ψ˜|Vˆ(cid:48) |Ψ(1)(cid:105)
E(2) = N . (8) for singly excited determinants with (cid:104)a| = (cid:104)0|Eˆ . The
(cid:104)Ψ˜|AP1roG(cid:105) i ia
bra basis then forms a biorthogonal basis that satisfies
the normalization conditions
Before we focus on possible choices of (Hˆ ) and Vˆ(cid:48)
0 N N
as well as (cid:104)Ψ˜|, we will define our projection space used (cid:104)ab|cd(cid:105)=δ +δ (14)
in the expansion of |Ψ(1)(cid:105) in eq. (4). Following previ- ij kl iajb,kcld jbia,kcld
ous PT models as well as linear CC corrections with an for doubly-excited determinants, while for the singly-
AP1roGreferencefunction,theprojectionspacewillcon- excited determinants we have
tainallpossibleexcitationswithrespecttoareferencede-
terminant. This reference determinant is, however, not (cid:104)a|b(cid:105)=δ . (15)
i j ia,jb
arbitrary, but restricted to the reference determined of
AP1roG, |0(cid:105) of eq. (1). We should emphasize that |0(cid:105) Choosingeq.(4)asthefirst-ordercorrectiontothewave-
is not equivalent to the Hartree–Fock determinant as in functionwith excitationoperatorsas defined ineqs. (10)
4
and(11), thecorrespondingperturbationamplitudesare while for single and double excitations we have
then obtained by solving
E(2) =2(cid:88)f ta+(cid:88)tab((cid:104)ij||ab(cid:105)+(cid:104)ij|ab(cid:105)). (22)
sd ia i ij
(cid:88)t (cid:104)Φ |(Hˆ ) |Φ (cid:105)+(cid:104)Φ |Vˆ(cid:48) |AP1roG(cid:105)=0, (16) ia iajb
p q 0 N p q N
p We should emphasize that this PT model is equivalent
to the PTa model presented in ref. 46. Note, however,
Note that the above equations depend on the partition
that pair excitations are not excluded in ref. 46 and
oftheHamiltonianHˆ intothezero-orderandthepertur-
that the full Hamiltonian Hˆ is taken as the perturba-
bation part as well as the projection manifold, but not
tion Hamiltonian so that E(2) = (cid:104)0|Hˆ|Ψ(1)(cid:105). Since the
on the choice of the dual state (cid:104)Ψ˜|. In the following, we
PTa amplitude equations for the pair excitations van-
will consider different partitionings of Hˆ as well as two
ish(astheyequaltheAP1roGamplitudeequations),the
choices for (cid:104)Ψ˜|.
PTa energy corrections is nonetheless determined from
eq. (20). For reasons of consistency, we will abbreviate
thesePTmodelsusingasingledeterminant(SD)asdual
1. Diagonal one-body zero-order Hamiltonian
and a diagonal (d) zero-order Hamiltonian as PT2SDd.
Thechoiceoftheexcitationoperatorwillbeindicatedin
First, we will restrict (Hˆ0)N to be a diagonal one- parentheses, that is, PT2SDd(d) for double excitations
electron operator. Specifically, we will choose (Hˆ ) to and PT2SDd(sd) for both single and double excitations.
0 N
be the diagonal of the inactive Fock operator b. Choosing (cid:104)AP1roG| as dual (cid:104)Ψ˜|. If (cid:104)AP1roG| is
(cid:32) occ (cid:33) takenasthedualstate(cid:104)Ψ˜|, wehavetoevaluatetermsas
Fˆ =(cid:88) h +(cid:88)((cid:104)pi||qi(cid:105)+(cid:104)pi|qi(cid:105)) a†a (cid:104)AP1roG|AP1roG(cid:105) in the energy expression, which be-
pq p q comes computationally intractable for large systems. In
pq i
order to arrive at a computationally feasible model, we
(cid:88)
= fpqa†paq, (17) willfollowref.46to,atleastpartially,eliminatetheover-
pq lap (cid:104)AP1roG|AP1roG(cid:105) in the PT equations and energy
expression. For that purpose, we redefine the zero-order
where(cid:104)pi||qi(cid:105)arethetwo-electronintegralsinphysicists’
Hamiltonian of eq. (18) by introducing the inverse of the
notation containing the Coulomb (cid:104)pi|qi(cid:105) and exchange
overlap (cid:104)AP1roG|AP1roG(cid:105) as a scaling factor,
(cid:104)pi|iq(cid:105) terms. Using normal-product operators, the zero-
order Hamiltonian reads (Hˆ ) =F¯d =(cid:88) fpp {a†a }
0 N N (cid:104)AP1roG|AP1roG(cid:105) p p
(Hˆ ) =Fˆd =(cid:88)f {a†a }, (18) p
0 N N p pp p p =(cid:88)f¯pp{a†pap}. (23)
p
while the perturbation becomes (cf. eq. (2))
By changing the zero-order Hamiltonian of eq. (18), we
Vˆ(cid:48) =Fˆo +Wˆ(cid:48) also have to adjust the corresponding perturbation part
N N N
given in eq. (19),
(cid:88) 1 (cid:88)
= f {a†a }+ (cid:48)(cid:104)pq|rs(cid:105){a†a†a a } (19)
pq p q 2 p q s r Vˆ(cid:48) =Fˆ −F¯d +Wˆ(cid:48)
p(cid:54)=q pqrs N N N N
=(cid:88)(f −f¯ δ ){a†a }+ 1 (cid:88)(cid:48)(cid:104)pq|rs(cid:105){a†a†a a }.
Eqs. (18) and (19) are then substituted into eq. (16) to pq pp pq p q 2 p q s r
solve for the PT amplitudes. Note that, in the case of a p,q pqrs
(24)
diagonal zero-order Hamiltonian, the PT amplitudes are
obtained from a set of uncoupled equations. To fully avoid the evaluation of the overlap
a. Restricting the dual (cid:104)Ψ˜| to (cid:104)0|. If (cid:104)Ψ˜| is re- (cid:104)AP1roG|AP1roG(cid:105) in the PT amplitude equations
stricted to the AP1roG reference determinant (cid:104)0|, we eq.(16),theinverseoftheAP1roGwavefunctionoverlap
canstraightforwardlyevaluatetheoverlap(cid:104)Φ0|AP1roG(cid:105), will be absorbed in the PT amplitudes. Thus, the
which equals 1 due to intermediate normalization of the first-order wavefunction contains scaled amplitudes,
AP1roG wavefunction. The expression for the second-
order energy correction E(2) given in eq. (8) thus simpli- |Ψ(1)(cid:105)=(cid:88) tp |Φ (cid:105)=(cid:88)t¯ |Φ (cid:105). (25)
fies to (cid:104)AP1roG|AP1roG(cid:105) p p p
p p
E(2) =(cid:104)0|Vˆ(cid:48) |Ψ(1)(cid:105). (20) By substituting eq. (23) into eq. (16) and introducing
N
the scaled PT amplitudes from eq. (25), we can elimi-
Specifically, the energy correction for Tˆ =Tˆ(cid:48) is given as natethewavefunctionoverlapfromthePTworkingequa-
2
tions. Furthermore,scalingthePTamplitudesbythein-
E(2) =(cid:88)tab((cid:104)ij||ab(cid:105)+(cid:104)ij|ab(cid:105)), (21) verseoftheoverlap(cid:104)AP1roG|AP1roG(cid:105)restoresthezero-
d ij order Hamiltonian of eq. (18) and we get (Hˆ ) = Fˆd.
iajb 0 N N
5
Note that the wavefunction overlap is still present in canbeevaluatedfromeq.(20). Note,however,thatonly
the perturbation part Vˆ(cid:48) . Due to the structure of the the doubly-excited determinants directly contribute to
N
zero-order wavefunction and the choice of the projection the energy correction. Since the perturbation Hamilto-
manifold, the diagonal part of the modified Fock oper- nian is a two-electron operator, the second-order energy
ator in eq. (24) does not contribute to the PT ampli- correction of the single excitations vanishes. Single exci-
tudeequationsandtheresultingperturbationreducesto tationscontributeindirectlythroughcouplingtothedou-
Vˆ(cid:48) = Fˆo + Wˆ(cid:48) . Since we have chosen (cid:104)AP1roG| as ble excitation manifold in the PT amplitude equations.
N N N
dual, the second-order energy correction is determined For both including and excluding the singles projection
from E(2) = (cid:104)AP1roG|Vˆ(cid:48) |Ψ(1)(cid:105) (cf. eq. (8)), where the manifold,thesecond-orderenergyE(2) isthuscalculated
N
first-order correction to the wavefunction is calculated from eq. (21). We will abbreviate the PT corrections
from the scaled PT amplitudes (eq. (25)), using an off-diagonal (o) zero-order Hamiltonian and a
single determinant for its dual state as PT2SDo, while
E(2) =(cid:88)t¯ (cid:104)AP1roG|Vˆ(cid:48) |Φ (cid:105). (26) the projection manifold will be indicated in parentheses
p N p
(d for doubles, sd for singles and doubles, respectively).
p
b. Choosing (cid:104)AP1roG| as dual (cid:104)Ψ˜|. Similar to the
Thesumintheaboveequationrunsoveralldeterminants PT methods with a diagonal Hˆ Hamiltonian, choos-
0
intheprojectionmanifold(doublyexcitedorsingly-and ing (cid:104)AP1roG| as dual forces us to evaluate the wave-
doubly-excited determinants). We should note that al- functionoverlap(cid:104)AP1roG|AP1roG(cid:105),whichbecomespro-
thoughwecanexactlyevaluatetheenergycorrectionand hibitive for large systems. In order to (partially) remove
thePTequations,thezero-andfirst-orderenergycorrec- thewavefunctionoverlapfromtheworkingequations,we
tions eq. (7) do not vanish. However, we can neglect the follow the procedure from above and introduce a scaled
weights of the (AP1roG) amplitudes beyond single pair zero-order Hamiltonian, where we have to modify both
excitations and assume that E(0) + E(1) ≈ 0. [46] We the diagonal Fd and off-diagonal Fo Fock operator,
N N
will label this PT model as PT2MDd as it uses a multi-
determinant (MD) wavefunction as dual and a diagonal (Hˆ0)N =F¯Nd +F¯No
(indd)iczaetreos-otrhdaetrtHheamexilctiotantiaionn. oFpuerrtahtoerrmcoonret,aiPnsTo2nMlyDddo(du)- =(cid:88)(cid:104)AP1roGfp|AqP1roG(cid:105){a†paq}=(cid:88)f¯pq{a†paq}.
ble excitations (without pairs), while in PT2MDd(sd) p,q p,q
both single and double excitations are included in Tˆ. (29)
Since we use a modified Fock operator as zero-order
Hamiltonian, we have to account for it in the definition
2. Off-diagonal one-body zero-order Hamiltonian of the perturbation part,
Vˆ(cid:48) =Fˆ −F¯d −F¯o +Wˆ(cid:48)
Next, we will consider an off-diagonal one-electron op- N N N N N
erator as the zero-order Hamiltonian. Similar to the PT =(cid:88)(f −f¯ ){a†a }+ 1 (cid:88)(cid:48)(cid:104)pq|rs(cid:105){a†a†a a }.
modelsemployingadiagonalzero-orderHamiltonian,we pq pq p q 2 p q s r
choose the Fock operator of eq. (17) as our Hˆ Hamilto- p,q pqrs
0 (30)
nian,
Analogous to PT2MDd-type methods, the scaling factor
(Hˆ0)N =FˆNd +FˆNo =(cid:88)fpq{a†paq}. (27) in the (Hˆ0)N Hamiltonian is absorbed in the PT ampli-
tudes (cf.eq. (25)) so that the zero-order Hamiltonian in
p,q
the PT amplitude equations eq. (16) contains only the
Then, the perturbation part Vˆ =Hˆ −Hˆ in its (shifted) unscaled Fock operator, (Hˆ ) = Fˆd +Fˆo, while the
0 0 N N N
normal-product form reads PT amplitudes t are replaced by the scaled amplitudes
p
t¯ . Although we eliminated the wavefunction overlap in
p
VˆN(cid:48) = 12 (cid:88)(cid:48)(cid:104)pq|rs(cid:105){a†pa†qasar}. (28) tohfethseecoPnTd-aomrdpelriteundeergeyqeuxaptiroensssi,o(cid:104)nAaPn1drionGt|hAeP(1Hˆro0)GN(cid:105)psatirltl
pqrs
remains in the perturbation part (cf. eq. (30)). In con-
To obtain the working equation for the PT amplitudes, trasttothePT2MDd-typemethodsdiscussedabove,the
wesubstituteeqs.(27)and(28)intoeq.(16). Notethat, perturbation Vˆ(cid:48) contains also modified off-diagonal ele-
N
incontrasttothePTmethodswithadiagonalzero-order mentsintheone-electronFockoperatorthatdonotvan-
Hamiltonian, the PT amplitudes are now obtained from ish in the PT equations. Therefore, we have to evaluate
a set of coupled equations. the wavefunction overlap (cid:104)AP1roG|AP1roG(cid:105) before we
a. Restricting the dual (cid:104)Ψ˜| to (cid:104)0|. If the dual state can determine the PT amplitudes. To obtain a compu-
is restricted to the reference determinant of the AP1roG tationally feasible model, we will approximate the over-
wavefunction, we benefit from the intermediate normal- lap (cid:104)AP1roG|AP1roG(cid:105)≈1+(cid:80) |ca|2, keeping only the
ia i
ization when evaluating the overlap (cid:104)0|AP1roG(cid:105). Analo- quadratic terms in the AP1roG amplitudes. This is usu-
gous to PT2SDd-type methods, the second order energy ally a good approximation as the AP1roG wavefunction
6
amplitudes are typically much smaller than 1 (ca (cid:28) 1). models scales as O(o2v2), where o is the number of oc-
i
Thus, the scaled off-diagonal Fock matrix elements can cupiedandv thenumberofvirtualorbitals,respectively.
be approximated by f¯ ≈ fpq . As for PT2MDd, Choosing an off-diagonal one-body zero-order Hamilto-
pq 1+(cid:80)ia|cai|2 nian (Hˆ ) increases the computational cost to O(o2v3)
the second-order energy E(2) can be determined from 0 N
eq.(26)withVˆ(cid:48) definedineq.(30). Notethatincontrast (see also Table I). Since we now have to solve a coupled
N set of linear equations iteratively, we have to consider an
to PT2SDo, the single excitations directly contribute to
additional prefactor. However, this prefactor is typically
the energy correction through both the one- and two-
muchsmallerthanv. Thus,PT2SDdandPT2MDdscale
electronoperatorsintheperturbationHamiltonian. The
PT models with an off-diagonal (o) Hˆ Hamiltonian and similar to AP1roG (or conventional electronic structure
0 methods like MP2), while in PT2SDo, PT2MDo, and
a multi-determinant (MD) wavefunction ((cid:104)AP1roG|) as
PT2b the computational cost increases by a factor of v.
dual will be abbreviated as PT2MDo, while the projec-
All PT approaches presented in this work are summa-
tionmanifoldwillbeagainindicatedinparentheses(dfor
rized in Table I for comparison.
doubles, sd for singles and doubles, respectively). Note
thatpair-excitationsareexcludedintheprojectionman-
ifold due to the orthogonality constraint.
B. Combining AP1roG with MP2
c. Relation to PTb theory. The PT2MDo model is
similar, but not equivalent, to the recently presented
PTb approach. [46] In the PTb method, the wavefunc- In our last PT model, we will combine conventional
tionoverlapintheperturbationHamiltonianisneglected, MP2 theory with AP1roG. Specifically, all wavefunction
that is, we assume 1 →0 so that f¯ →0. amplitudes corresponding to singly and doubly excited
(cid:104)AP1roG|AP1roG(cid:105) pq
determinantswithrespectto|0(cid:105)aredeterminedfromthe
By neglecting the scaled components of the Fock oper-
MP2 amplitude equations. To prevent double counting
ator, the perturbation Hamiltonian reduces to the full
quantum chemical Hamiltonian Vˆ(cid:48) = Hˆ(cid:48) . The second- of the correlation contribution associated with electron
N N pairs,weomittheamplitudesofallpair-exciteddetermi-
order energy correction is then given as
nantsintheMP2equations. InMP2,thePTamplitudes
E(2) =(cid:88)t¯ (cid:104)AP1roG|Hˆ(cid:48) |Φ (cid:105). (31) are determined from a set of uncoupled equations,
p N p
p (cid:88)(cid:48)t (cid:104)Φ |(Hˆ ) |Φ (cid:105)+(cid:104)Φ |Vˆ |0(cid:105)=0, (32)
p q 0 N p q N
Furthermore, in PTb theory, pair excitations are not p
excluded in the projection manifold and the first-order
correction of the wavefunction contains all double ex- where |0(cid:105) is the AP1roG reference determinant
citations. These pair excitations do not contribute to (cf. eq. (16)). The second-order energy corrections then
the energy corrections eq. (31). Their contribution van- accounts for correlation contributions beyond electron-
ishes as the corresponding terms in eq. (31) equal the pairs and can be determined from eq. (22), which con-
AP1roG amplitude equations. However, pair excitations tains both single and double excitations with respect to
indirectly enter the energy correction by coupling to the |0(cid:105). We will label this methods as AP1roG-MP2.
remaining PT amplitudes in the PT equations as well
as to the pair excitations of the AP1roG model. In this
work, we have extended the original PTb model by in- C. Linearized Coupled Cluster Corrections
cludingalsosingleexcitationsintheprojectionmanifold.
Furthermore, we also investigate how the pair excita- Besides PT-type methods, dynamic correlation effects
tions in the projection manifold influence the PTb en- can be built in the AP1roG wavefunction a posteriori
ergy correction. For that purpose we have excluded the using an exponential coupled cluster ansatz, [51]
pairexcitationsintheprojectionmanifoldwhenoptimiz-
ing the PTb amplitudes. To emphasize the order of the |Ψ(cid:105)=exp(Tˆ)|AP1roG(cid:105), (33)
energy correction in PTb, these PT models are abbrevi-
ated as PT2b, while the projection manifold is indicated where Tˆ = (cid:80) t τˆ is a general cluster operator. The
ν ν ν
in parentheses (d for doubles, sd for singles and doubles, corresponding time-independent Schr¨odinger equation
d\p for doubles without pairs, sd\p for singles and dou- then reads
bles excluding pair excitations).
d. Comment on computational scaling of the PT Hˆ exp(Tˆ)|AP1roG(cid:105)=Eexp(Tˆ)|AP1roG(cid:105). (34)
models. The computational scaling of all PT correc-
tions discussed above is determined by the first term of Insteadofconsideringthefullexpansionoftheclusterop-
eq. (16). Note that we can introduce suitable intermedi- erator in the exponential function, we can include terms
ates for the second term in eq. (16) so that summations only linear in Tˆ. If we truncate the Baker–Campbell–
are performed only once in the beginning of the calcula- Hausdorff expansion after the second term,
tion. If(Hˆ ) isrestrictedtobeadiagonalone-bodyop-
0 N
erator, the computational cost of the corresponding PT exp(−Tˆ)Hˆ exp(Tˆ)≈Hˆ +[Hˆ,Tˆ], (35)
7
TABLE I. Summary of PT models with zero-order Hamiltonian Hˆ , perturbation Vˆ, dual (cid:104)Ψ˜|, and excitation operator Tˆ. All
0
operators are defined in the text. The computational scaling is given in the last column.
Model Hˆ Vˆ (cid:104)Ψ˜| Tˆ scaling
0
PT2SDd Fˆd Fo +Wˆ(cid:48) (cid:104)0| Tˆ(cid:48), Tˆ +Tˆ(cid:48) O(o2v2)
N N N 2 1 2
PT2MDd Fˆd Fo +Wˆ(cid:48) (cid:104)AP1roG| Tˆ(cid:48), Tˆ +Tˆ(cid:48) O(o2v2)
N N N 2 1 2
PT2SDo Fˆd +Fˆo Wˆ(cid:48) (cid:104)0| Tˆ(cid:48), Tˆ +Tˆ(cid:48) O(o2v3)
N N N 2 1 2
PT2MDo Fˆd +Fˆo Fˆ −F¯ +Wˆ(cid:48) (cid:104)AP1roG| Tˆ(cid:48), Tˆ +Tˆ(cid:48) O(o2v3)
N N N N N 2 1 2
PT2b Fˆd +Fˆo Hˆ(cid:48) (cid:104)AP1roG| Tˆ , Tˆ(cid:48), Tˆ +Tˆ , Tˆ +Tˆ(cid:48) O(o2v3)
N N N 2 2 1 2 1 2
we obtain the linearized coupled cluster Schr¨odinger corrections of the dynamic correlation energy were per-
equation formed in our in-house quantum chemistry code. The
single-point energies of all investigated reactions were
(Hˆ +[Hˆ,Tˆ])|AP1roG(cid:105)=E|AP1roG(cid:105). (36) evaluatedusingthecc-pVDZ,cc-pVTZ,andcc-pVQZba-
sis sets of Dunning [68] and extrapolated to the basis set
The cluster amplitudes t are determined by solving a
ν limit. Two sets of calculations were performed, one with
linear set of coupled equations obtained from multiply-
thecanonicalHartree–Fockorbitalsandanotherwiththe
ing the above equation from left by determinants of the
AP1roG optimized orbital basis.
projection manifold (cid:104)ν|,
The dissociation curve of the C (ground state)
2
moleculewasdeterminedforDunning’saug-cc-pVTZba-
(cid:104)ν|(Hˆ +[Hˆ,Tˆ])|AP1roG(cid:105)=0. (37)
sis set, while for N (ground state) and BN (first ex-
2
cited state) the cc-pVTZ basis set was employed. All
The energy can be calculated by projecting against the
the dynamic energy corrections on top of AP1roG were
reference determinant of |AP1roG(cid:105),
calculated in the AP1roG optimized orbital basis. The
(cid:104)0|(Hˆ +[Hˆ,Tˆ])|AP1roG(cid:105)=E. (38) points of the potential energy curves of all investigated
diatomics were used for a subsequent generalized Morse
Specifically, in the AP1roG-LCCD approach, the clus- function [69] fit to obtain the equilibrium bond lengths
ter operator is restricted to double excitations with re- (re) and potential energy depths (De). The harmonic
spect to |0(cid:105) as defined in eq. (10), while the AP1roG- vibrational frequencies (ωe) were calculated numerically
LCCSD method also single excitations are included and using the five-point finite difference stencil [70].
hence Tˆ = Tˆ + Tˆ(cid:48) (cf. eqs.(10) and (11)). For more
1 2
information concerning the choice of the cluster opera-
B. Extrapolation to the basis set limit
tor and the projection manifold, we refer the reader to
ref. 51. Note that the computational scaling of AP1roG-
LCCD/LCCSDissimilartothecomputationalscalingof ThebasissetlimitoftheHartree–Fockenergywasob-
conventional LCCD/LCCSD methods, that is, O(o2v4). tainedbyfittinganexponentialfunctionoftheform[71]
ESCF(X)=ESCF+aexp(−bX) (39)
∞
III. COMPUTATIONAL DETAILS
to the Hartree–Fock energies obtained in the cc-pVDZ,
cc-pVTZ,andcc-pVQZbasissets. Intheaboveequation
A. Geometry optimization and single point
X indicates the cardinal number of the basis set (2 for
calculations
D, 3 for T, etc.). For all correlation calculations, the
basis set limit of the correlation energy was obtained by
All molecular structures used in the calculations of re-
a two-point fit using the fit function
action energies were fully optimized in the TURBO-
MOLE7.0 software package [57], employing the B3LYP Ecorr(X)=Ecorr+aX−3, (40)
∞
exchange–correlation functional [58, 59] and a polarized
valence triple-ζ basis set [60, 61]. The resulting geome- as suggested in refs. [71, 72]. In the above equation,
tries in xyz format are included in the Supporting Infor- Ecorr(X) indicates the correlation energy of a given
mation. method defined as Ecorr(X) = Etot(X)−ESCF(X). In
The single-point MP2, BCCD, and CCSD(T) calcula- the case of dynamic correction on top of AP1roG, the
tions have been carried out in the Molpro2012.1 soft- correlationenergyisthusasumofthecorrelationenergy
warepackage[62,63],whiletheNWChemsoftwaresuite of AP1roG and of the dynamic correlation model. Note,
(in combination with the tensor contraction engine [64– however, that extrapolating the sum of correlation ener-
67])wasusedforLCCD,LCCSD,andCR-CCSD(T)cal- gies or each contribution separately yields similar basis
culations. AllAP1roGcalculationsincludingaposteriori set limits of the (total) correlation energy Ecorr. For all
∞
8
correlation calculations, only the correlation energies of
TABLEII.Spectroscopicconstantsforthedissociationofthe
the cc-pVTZ (X = 3) and cc-pVQZ (X = 4) basis sets
C , N , and BN molecule for different quantum chemistry
were employed in the fitting procedure. 2 2
methods. The differences are with respect to MRCI-SD ref-
erence data [73] for C and N , and with respect to MRCI-
2 2
SD+Qreferencedata[74]forBN.E : groundstateenergyat
e
IV. RESULTS AND DISCUSSION r .
e
The previously discussed dynamic correlation mod- Method Ee [Eh] re [˚A] De [kmcoall] ωe [cm−1]
els are benchmarked against spectroscopic constants AP1roG −75.58569 1.227(−0.025) 132.9(−8.2) 1780(−56)
AP1roG-PT2SDd(d) −75.80222 1.251(+0.001) 160.4(+19.3) 1889(+53)
of multiply bonded diatomics and thermochemical
AP1roG-PT2SDd(sd) −75.81041 1.249(−0.003) 154.6(+13.5) 1915(+79)
data for compounds containing main-group elements. AP1roG-PT2MDd(d) −75.77832 1.239(−0.013) 121.6(−19.5)∗ 1940(+104)
AP1roG-PT2MDd(sd) −75.78630 1.238(−0.014) 115.6(−25.5)∗ 1963(+127)
Note that the former test systems are dominated by
AP1roG-PT2SDo(d) −75.81778 1.242(−0.010) 156.5(+15.4) 1919(+83)
static/nondynamic electron correlation effects, while dy- AP1roG-PT2SDo(sd) −75.81783 1.242(−0.010) 156.4(+15.3) 1919(+83)
namicelectroncorrelationbecomesimportantforthelat- AP1roG-PT2MDo(d) −75.79032 1.231(−0.021) 125.9(−15.2)∗ 2010(+174)
C2 AP1roG-PT2MDo(sd) −75.79108 1.230(−0.022) 119.9(−21.2)∗ 2019(+183)
ter. Both test sets allow us to assess the accuracy of the AP1roG-PT2b(d) −75.78350 1.235(−0.017) 127.2(−13.9)∗ 1938(+102)
proposed dynamic corrections for strongly- and weakly AP1roG-PT2b(sd) −75.79400 1.228(−0.024) 113.0(−28.1)∗ 2049(+213)
AP1roG-PT2b(d\p) −75.78381 1.231(−0.021) 123.9(−17.2)∗ 2016(+180)
correlated molecular systems. AP1roG-PT2b(sd\p) −75.79370 1.228(−0.024) 112.9(−28.2)∗ 2048(+212)
AP1roG-LCCD[51] −75.81125 1.240(−0.012) 139.3(−1.8) 1916(+80)
AP1roG-LCCSD[51] −75.81257 1.240(−0.012) 143.0(+1.9) 1926(+90)
NEVPT2[32] −75.78829 1.244(−0.008) 148.0(+6.9) 1886(+50)
A. Molecular systems dominated by CR-CCSD(T) −75.80484 1.242(−0.010) 152.1(+9.0) 1989(+153)
MRCI-SD[73] −75.78079 1.252 141.1 1836
static/nondynamic correlation
AP1roG −109.12686 1.087(+0.083) 255.9(−38.0) 2435(+94)
AP1roG-PT2SDd(d) −109.37326 1.094(−0.010) 249.8(−31.9) 2301(−40)
Wewillfirstbrieflydiscusstheperformanceandtrends AP1roG-PT2SDd(sd) −109.37383 1.093(−0.011) 233.4(+15.5) 2299(−42)
AP1roG-PT2MDd(d) −109.35918 1.104(+0.000) 210.3(−7.6)∗ 2744(+403)
of the proposed PT corrections in modeling dissocia- AP1roG-PT2MDd(sd) −109.35988 1.105(+0.001) 202.0(−15.1)∗ 2779(+438)
tion processes of three diatomic molecules, C , N , and AP1roG-PT2SDo(d) −109.39838 1.110(+0.006) 246.9(+29.0) 2278(−63)
2 2
AP1roG-PT2SDo(sd) −109.39838 1.108(+0.004) 248.8(+30.9) 2254(−87)
BN.Theirelectronicstructureshavebeenextensivelydis- AP1roG-PT2MDo(d) −109.38106 1.089(−0.015) 220.8(+2.9)∗ 2200(−141)
cussedintheliteratureandwerefertheinterestedreader AP1roG-PT2MDo(sd) −109.38132 1.090(−0.014) 216.7(−1.0)∗ 2203(−138)
to, for instances, refs. [73–77] for more details. Table II N2 AP1roG-PT2b(d) −109.37490 1.092(−0.012) 219.2(+1.3)∗ 2246(−95)
AP1roG-PT2b(sd) −109.37582 1.095(−0.009) 207.1(−10.8)∗ 2272(−69)
summarizes the spectroscopic constants for the C2, N2, AP1roG-PT2b(d\p) −109.37431 1.092(−0.012) 218.9(+1.0)∗ 2240(−101)
and BN molecules. The corresponding potential energy AP1roG-PT2b(sd\p) −109.37525 1.094(−0.010) 206.8(−11.1)∗ 2263(−78)
AP1roG-fLCCD −109.39671 1.102(−0.002) 206.7(−11.2) 2334(−7)
surfaces are shown in the Supplementary Information. AP1roG-fLCCSD −109.39871 1.103(−0.001) 211.6(−6.3) 2337(−4)
For C and N , MRCI-SD data by Peterson et al. [73] is CASSCF[73] −109.13190 1.106(+0.002) 211.6(−6.3) 2340(+1)
2 2 CR-CCSD(T) −109.39632 1.101(−0.003) 229.7(+11.8)∗ 2397(+56)
taken as reference, while MRCI-SD+Q results are used
MRCI-SD[73] −109.36162 1.104 217.9 2341
as reference for the BN molecule [74]. AP1roG −79.07582 1.304(+0.019) 114.3(−40.1) 1630(−52)
In general, PT2SDd and PT2SDo provide equilibrium AP1roG-PT2SDd(d) −79.28333 1.283(−0.002) 181.3(+26.9) 1751(+68)
AP1roG-PT2SDd(sd) −79.29473 1.271(−0.014) 166.1(+11.7)∗ 1784(+102)
bond lengths and vibrational frequencies that agree well
AP1roG-PT2MDd(d) −79.26001 1.273(−0.012) 160.7(+6.3) 1778(+96)
with the corresponding reference values for all investi- AP1roG-PT2MDd(sd) −79.27571 1.264(−0.021) 144.5(−9.9) 1760(+78)
AP1roG-PT2SDo(d) −79.29432 1.282(−0.003) 169.9(+15.5) 1718(+36)
gated diatomics. Furthermore, addition of single excita-
AP1roG-PT2SDo(sd) −79.29427 1.281(−0.004) 169.7(+15.3)∗ 1719(+37)
tions in the excitation operator Tˆ does not increase the AP1roG-PT2MDo(d) −79.26814 1.274(−0.011) 158.5(+4.1)∗ 1743(+61)
accuracy of PT2SDd(d) or PT2SDo(d). In contrast to AP1roG-PT2MDo(sd) −79.27584 1.269(−0.016) 139.0(−15.4)∗ 1725(+43)
BN AP1roG-PT2b(d) −79.25960 1.273(−0.012) 152.6(−1.8)∗ 1739(+57)
PT2SDdandPT2SDo,theperformanceofPTmodelsus- AP1roG-PT2b(sd) −79.27961 1.260(+0.025) 141.9(−12.5)∗ 1746(+64)
ingamulti-determinantdualstate(PT2MDd,PT2MDo, AP1roG-PT2b(d\p) −79.25938 1.273(+0.012) 156.6(+2.2)∗ 1741(+59)
AP1roG-PT2b(sd\p) −79.27934 1.261(−0.024) 142.3(−12.1)∗ 1746(+64)
and PT2b) is less satisfying. Specifically, these meth-
AP1roG-LCCD[51] −79.28205 1.277(−0.008) 174.1(−19.8) 1734(+52)
ods provide equilibrium bond distances and vibrational AP1roG-LCCSD[51] −79.28975 1.275(−0.010) 178.6(+24.2) 1749(+67)
CCSD(T)[78] — 1.244(−0.008) — 1886(+50)
frequencies that deviate most from MRCI-SD/MRCI-
CR-CCSD(T) −79.28080 1.275(−0.010) 164.9(+10.5) 1726(+44)
SD+Q referencedata. Similar to PT2SDdand PT2SDo, RMRCCSD(T)[78] — 1.284(−0.001) — 1681(+1)
addition of single excitations worsens spectroscopic con- MRCI-SD+Q[74] −79.26794 1.285 154.4 1682
stants and hence the excitation operator should be re-
stricted to double excitations only. Furthermore, in- ∗ Estimateddissociationenergy.
cluding pair excitations in Tˆ (PT2b-type methods) does
not significantly affect spectroscopic constants and both
models (with and without pair excitations) yield similar kcal/mol from MRCI-SD/MRCI-SD+Q reference data.
values for r and ω . Note that most PT corrections (PT2MD and PT2b) di-
e e
In contrast to equilibrium bond distances and vibra- verge in the vicinity of dissociation and hence only an
tionalfrequencies,theaccuratepredictionofdissociation estimated dissociation energy is given in Table II (indi-
energies D is more challenging. In general, none of the cated by the ∗ in the Table).
e
proposed PT models can reliably predict potential en- The performance of the LCCD/LCCSD correction on
ergy well depths, which differ by approximately 15 to 30 top of AP1roG is more robust in predicting spectro-
9
of all error measures). Compared to Hartree–Fock data,
TABLE III. List of reactions containing closed-shell
AP1roG predicts reaction energies that differ more from
molecules.
the CR-CCSD(T) reference. In general, AP1roG under-
F +H −→2HF {1}
2 2 estimates reaction energies and yields an RMSE of ap-
F O+H −→F +H O {2}
2 2 2 2 proximately15kcal/mol(seealsoTableIVforcc-pVQZ).
H O +H −→2H O {3}
2 2 2 2 Thisdiscrepancycanbeattributedtoelectroncorrelation
N +3H −→2NH {4}
2 2 3 effects beyond electron pairs that are not accounted for
N O+H −→N +H O {5}
2 2 2 2
C H +H −→C H {6} in the AP1roG method. Inclusion of open-shell determi-
2 2 2 2 4
BH +3HF−→BF +3H {7} nants in the electronic wavefunction improves reaction
3 3 2
CO+H O−→CO +H {8} energies. However, the performance of the PT models
2 2 2
CO+3H −→CH +H O {9} strongly depends on the choice of the zero-order Hamil-
2 4 2
2BH −→B H {10} ton and the perturbation. Specifically, all PT methods
3 2 6
2H2O−→(H2O)2 {11} with a diagonal one-electron Hamiltonian (PT2SDd and
2HF−→(HF) {12}
2 PT2MDd) provide slightly improved reaction energies
HCOOH−→CO +H {13}
2 2 with an RME of around 13 kcal/mol (see also Table IV
CO+CH −→CH CHO {14}
4 3 for cc-pVQZ). Yet, these methods are strongly basis-set
2NH −→(NH ) {15}
3 3 2 dependent if optimized AP1roG orbitals are used as or-
bitalbasis: theRMSEgraduallyincreaseswiththebasis
set size (see Figure 1). An a posteriori MP2 correction
in the AP1roG basis predicts the largest RMSE of ap-
scopicconstantsforallinvestigateddiatomics. Although
proximately 16 kcal/mol and therefore does not allow us
AP1roG-LCCSDyieldsr andω thatdeviatemorefrom
e e
to accurately model reaction energies. Choosing an off-
MRCI-SD/MRCI-SD+Q reference values compared to
diagonal one-electron zero-order Hamiltonian (PT2SDo
PT2SDd and PT2SDo results, an LCC correction al-
andPT2MDo)inthePTmodelfurtherreducestheerror
lows us to reliably model dissociation energies for C
2
withrespecttotheCCSD(T)referencewithanRMSEof
and N , while the errors in D increase for BN (dif-
2 e
approximately 4-5 kcal/mol. The error in reaction ener-
ferences are 2 kcal/mol for C , 6 kcal/mol for N , and
2 2
gies predicted by PT corrections can be minimized if the
24 kcal/mol for BN). Note that for the N molecule the
2
perturbation Hamiltonian of PT2MDo (see Table I) is
coupledclusterdoublesamplitudesforexcitationswithin
the valence shell (vs) have been kept frozen in the vicin- replaced by the full quantum chemical Hamiltonian HˆN
ity of dissociation (r > 2.0 ˚A with tvs = tvs ) to as in all PT2b-type methods (RMSE ≈ 3.0 kcal/mol),
r>2.0 r=2.0
prevent divergencies in the dissociation limit and to ob- where the accuracy is similar to MP2 theory (see Ta-
tain a smooth potential energy surface. This is indi- ble IV). Note that inclusion of single excitations in the
cate by fLCCD/fLCCSD in Table II. The corresponding excitation operator (Tˆ = Tˆ1 + Tˆ2(cid:48)) does not affect the
AP1roG-fLCCSDspectroscopicconstantsagreewellwith accuracyofthePTcorrectionsinpredictingreactionen-
CASSCFresultsbyPetersonetal.[73]andareingeneral ergies(seealsoTableIV).OnlytheperformanceofPT2b
closest to MRCI-SD reference values. is slightly improved.
For all investigated basis sets, the best agreement
withCR-CCSD(T)referencedataisobtainedforthelin-
B. Molecular systems dominated by dynamic earized coupled cluster correction including singles and
correlation doubles (see Figure 1) with an RMSE of 1.5 kcal/mol
and a mean absolute error of 1.2 kcal/mol. Note that
Our second test case includes 15 reactions contain- the errors of AP1roG-LCCSD is similar to those of con-
ing main-group elements summarized in Table III. For ventional CC methods like BCC, CCSD, CCSD(T), and
all studied systems, the reaction energies obtained from LCCD. Only LCCSD performs worse with an RMSE of
CR-CCSD(T) calculations are taken as reference values. 2.5 kcal/mol and a mean absolute error of 1.8 kcal/mol.
Furthermore,wewillfocusontwodifferentmolecularor- We should emphasize that, unlike the LCC method with
bitalbasissetsusedinourAP1roGcalculationsandinall an Hartree–Fock reference function, single excitations in
a posteriori corrections for dynamic correlation with an the LCC corrections on top of AP1roG are particularly
AP1roG reference function; these are optimized natural important and allow us to approach chemical accuracy
AP1roG orbitals and canonical Hartree–Fock orbitals. (approximately1kcal/mol)whenmodelingthermochem-
istry of main-group compounds.
1. Optimized natural AP1roG orbitals
2. Canonical Hartree–Fock orbitals
Figure 1 shows the reaction energies and root mean
squareerrors(RMSEs)withrespecttoCR-CCSD(T)ref- If canonical Hartree–Fock orbitals are chosen as or-
erencedata(greybars)forallinvestigatedmethodsinthe bitalbasis,wecanusethediagonalstructureoftheFock
optimized AP1roG basis (see Table IV for the definition matrixandsimplifythePTamplitudeequations. Specif-
10
30
cc−pVTZ
cc−pVQZ
CBS
20
ol]
m
al/ 10
c
k
T) [
D(
S
C 0
C
−
R
C
wrt 8
or −10 6
Err 4
2
0
−20
−2
−4
−6
−30
HF AP1roAGP1roPGT−2MSPDP2dT(2dS)DPdT(2sMd)DPdT(2dM)DPdT(2sSdD)PoT(2dS)DPoT(2sMd)DPoT(2dM)DPoT(2sbd()dP)T2b(sPdT)2b(dP\Tp2)b(sAdP\p1)roAGP−1LrCoCMG−PL2CCBCCCCSDCCSDL(TC)CDLCCSD
D S
D
FIG. 1. Errors in the optimized AP1roG orbital basis.
icallyforallPTmodels, thePTamplitudeequationsare dual state affects the performance of the PT models.
transformed into a set of uncoupled linear equations as
IncontrasttotheoptimizednaturalAP1roGbasis, all
Fˆo = 0 (cf. Table I). Furthermore, single excitations do
N PT models yield similar reaction energies with a mean
not contribute (either directly in the energy expression
absolute error of approximately 3 kcal/mol and a RMSE
or indirectly through coupling to the doubles manifold)
ofabout3.5kcal/mol,irrespectiveoftheexcitationoper-
if the AP1roG reference determinant is taken as dual atorTˆ(seeTableVandFigure2). Onlythecombination
(cid:104)Ψ˜|. Thus, the PT2SDd and PT2SDo methods become
of MP2 (for open-shell configurations) and AP1roG (for
equivalent(withandwithoutsingleexcitations)andonly
closed-shell configurations) performs worse resulting in
the PT2SDd(d) results are shown in Table V and Fig-
an RMSE of 4.4 kcal/mol. As observed above, the LCC
ure 2. Similarly, the PT2MDd and PT2MDo models are
correction outperforms all PT models if both single and
equivalent to each other. In contrast to PT2SD-type
double excitations are included in the cluster operator
approaches, single excitations do contribute directly to
(cf. Table V and Figure 2), where the RMSE reduces
the second-order energy correction as their contribution
to 1.0 kcal/mol. Furthermore, the accuracy of the PT
is determined from eq. (26), which also contains one-
methods does not deteriorate for increasing sizes of the
electron terms. Finally, in PT2b-type methods, the cou-
basis set (see Figure 2).
pling to pair-excited determinants vanishes in the PT
amplitude equations as they are coupled through off- Finally,weshouldnotethattheperformanceofalldy-
diagonal Fock matrix elements. Furthermore, the per- namic correlation corrections depends on the molecular
turbation Hamiltonian Hˆ(cid:48) of PT2b reduces to Vˆ(cid:48) and orbital basis. The errors with respect to CR-CCSD(T)
N N
both the amplitude equations and second-order energy referencedatasignificantlydecreaseforPTmethodswith
expressionareequivalenttothoseinPT2MD-typemeth- asingledeterminantasdual(cid:104)0|whencanonicalHartree–
ods. If (cid:104)AP1roG| is chosen as dual, all PT methods co- Fockorbitalsareusedasmolecularorbitalbasis. Incases
incide in the canonical Hartree–Fock basis. Therefore, of (cid:104)AP1roG|, the dependence on the molecular orbital
Table V and Figure 2 show only the PT2MDd(d) and basis is less severe. While reaction energies obtained by
PT2MDd(sd) results. The canonical Hartree–Fock basis PT2MDd/PT2MDo in the canonical Hartree–Fock ba-
thus allows us to directly assess how the choice of the sis slightly improve reducing the error by a factor of 2,
the mean absolute error and RMSE of PT2b remain al-
