Second-Order Perturbation Theory for Fractional Occupation

Subscriber access provided by GAZI UNIV

Article

Second-order perturbation theory for fractional occupation systems: Applications to ionization potential and electron affinity calculations Neil Qiang Su, and Xin Xu J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00197 • Publication Date (Web): 24 Mar 2016 Downloaded from http://pubs.acs.org on March 26, 2016

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 free 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 accessible to all readers and 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.

Journal of Chemical Theory and Computation 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 42

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

Second-order perturbation theory for fractional occupation systems: Applications to ionization potential and electron affinity calculations

Neil Qiang Su, Xin Xu* Collaborative Innovation Center of Chemistry for Energy Materials, Shanghai Key Laboratory of Molecular Catalysis and Innovative Materials, MOE Laboratory for Computational Physical Science, Department of Chemistry, Fudan University, Shanghai, 200433, China

ACS Paragon Plus Environment


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

ABSTRACT: Recently, we have developed an integration approach for the calculations of ionization potentials (IPs), and electron affinities (EAs) of molecular systems at the level of second-order Møller-Plesset (MP2)

[N. Q. Su, and X. Xu, J. Chem. Theory Comput. 11, 4677

(2015)], where the full MP2 energy gradient with respect to the orbital occupation numbers was derived, but only at integer occupations. The theory is completed here to cover the fractional occupation systems, such that Slater’s transition state concept can be used to have accurate predictions of IPs and EAs. Antisymmetrized Goldstone (ASG) diagrams have been employed for interpretations and better understanding of the derived equations, where two additional rules were introduced in the present work specifically for hole or particle lines with fractional occupation numbers.

*Corresponding author: [email protected]


Page 2 of 42

Page 3 of 42

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


I. INTRODUCTION In our previous work,1 we have derived the full second-order Møller-Plesset (MP2)2 energy gradients with respect to the orbital occupation numbers, ntλ , for integer occupation systems, where orbitals can only be integrally occupied with ntλ being either 1 or 0. Based on the derived gradients at different approximations1,3,4 and the integration approach,4,5 we have proposed one- and two-point equations for ionization potential (IP) and electronic affinity (EA) calculations.1 Here we derive the full theory to cover the fractional occupation systems. For IP calculations, the integration path connecting the initial ( nk0 = 1 ) and final ( n1k = 0 ) states may be chosen as nkλ = 1 − λ niλ = 1 for i ≠ k

(1)

λ

na = 0 for a ≠ k For EA calculations, the integration path connecting the initial ( nc0 = 0 ) and final ( n1c = 1 ) states may be chosen as ncλ = λ niλ = 1 for i ≠ c

(2)

λ

na = 0 for a ≠ c By convention, we use indices i, j, and l to stand for the occupied canonical Hartree-Fock (HF) spin-orbitals, a, b, and d for the virtual canonical HF spin-orbitals. We also use indices p, q, r, s, and t to stand for any canonical HF spin-orbitals which can be either occupied or virtual. In particular, k or c is used to represent a partially occupied orbital, which, as prescribed in Eqs. 1 and 2, is recognized as both an occupied and a virtual orbital indexed by both i (or j and l) and a (or b and d). In general, we use t to indicate the orbital of interest with a varying occupation. While t = k is a hole index used for IP calculations, labelling the orbital from which the electron is to be ionized, t = c is a



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 4 of 42

particle index used for EA calculations, labelling the orbital to which the electron is to be added. These notations will be used in the present work for fractional occupation systems as described below, which are similar to those in our previous work for integer occupation systems.1 As shown in the integration approach,1,4,5 IP and EA can be computed as 1

dE ( λ )

0

dnkλ

IP =  E (1) − E ( 0 )  = − ∫

dλ

(3)

and 1

dE ( λ )

0

dncλ

EA = −  E (1) − E ( 0 )  = − ∫

dλ

(4)

The first equal sign in Eq. 3 or 4 actually indicates the ∆ method as the total energy differences of the final and initial states, while the second equal sign suggests a direct method. Previously,1 the energy gradients,

dE ( λ ) dntλ

, were evaluated at either λ = 0 or λ = 1, leading to the one-point equations which

are exact up to the first order in terms of coupling parameter λ.1 By taking the average of both one-point results at λ = 0 and λ = 1, a two-point equation was obtained, which improved the accuracy for IP and EA calculations up to the second order.1 It has long been recognized and widely-used in density functional theory (DFT) methods, as introduced by Slater in his Xα approximation,6,7 the difference of two state energies can be computed by doing a single calculation in the transition state. That is, the state with an occupation number halfway between the initial and final states. It was shown that the one-point equation at λ = 1/2 has already led to results for energy differences that are exact up to the second order.8,9 Extension to a two-point equation at λ = (0, 2/3) has also been made,8 which formally improved the accuracy for energy difference calculations up to the third order.8,9 However, for such methods to be applicable to wavefunction-based correlated methods, such as MP2,2 as well as the newly developed doubly


Page 5 of 42

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


hybrid density functionals, such as XYG3,10-13 one needs to develop method where the energy gradients,

dE ( λ ) dntλ

, can be evaluated explicitly in fractional occupation systems. This is the subject

of the present work. The reminder of this paper is organized as follows. In Sec. II, we briefly discuss the MP2 method with fractional occupations at different levels of approximation.1,3,4 We then derive the full MP2 energy derivatives with respect to the orbital occupation numbers explicitly for fractional occupation systems. Antisymmetrized Goldstone (ASG) diagrams14 are employed to interpret the derived equations for better understanding. Computation details are summarized in Sec. III, and test calculations are presented in Sec. IV for some selective atoms and molecules. Some concluding remarks are given in Sec. V.

II. THEORY A. The MP2 method for fractional occupation systems at different levels of approximation The total energy of MP2 is given by EMP 2 ( λ ) = E HF ( λ ) + Ec , MP 2 ( λ )

(5)

where the HF total energy may be expressed as 1 EHF ( λ ) = ∑ nλp p λ hˆ p λ + ∑ nλp nqλ p λ q λ p λ q λ 2 pq p

(6)

and the corresponding MP2 correlation energy as

Ec , MP 2 ( λ ) = 1 ∑ 4 pqrs

nλp nqλ (1 − nrλ )(1 − nsλ ) pλ q λ r λ s λ

ε pλ + ε qλ − ε rλ − ε sλ

2

(7)

Here hˆ is the one-electron operator made of the kinetic operator −1 2 ⋅ ∇ˆ 2 and the external potential υˆext usually coming from the nuclei, while

{ε } λ p

refer to the HF orbital energies from the



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 42

self-consistent-field (SCF) calculation at a given λ with a chosen set of occupation numbers of either 1 or 0. The corresponding canonical HF spin orbitals {φ pλ } are denoted by

{n } λ p

{ p } , which λ

define the two-electron repulsion integrals as:

p λ qλ r λ s λ = ∫∫ φ pλ ( x1 )φrλ ( x1 ) r12−1φqλ ( x2 ) φsλ ( x2 ) dx1dx2

(8)

with x1 and x2 being combined spatial and spin coordinates. Eq. 8, in turn, defines its antisymmetrized combination:

p λ qλ r λ s λ = p λ q λ r λ s λ − p λ qλ s λ r λ

(9)

Eqs. 5-9 may be extended to partial occupation systems,1,3,4 such that one may differentiate the HF total energy with respect to the occupation number ntλ , leading to

dEHF ( λ ) λ

dnt

= ε tλ = t λ hˆ t λ + ∑ nλp p λ t λ p λ t λ

(10)

p

Eq. 10, when combined with Eqs. (3) and (4), gives IP and EA calculations at the HF level.1 One may also differentiate the MP2 correlation energy with respect to the occupation number

ntλ , arriving at1

dEc ,MP 2 ( λ ) ∂Ec ,MP 2 ( λ ) = dntλ ∂ntλ ∂Ec , MP 2 ( λ ) ∂ε λp +∑ ∂ε pλ ∂ntλ p

(I) (II)

(11)

 ∂Ec , MP 2 ( λ ) ∂ε pλ ∂Ec, MP 2 ( λ )  dφqλ + ∑ ∑ +  λ (III) ∂ε pλ ∂φqλ ∂φqλ q  p  dnt Part I accounts for the explicit dependence of Ec , MP 2 ( λ ) on the occupation numbers, giving1,3,4

(I) =

−

1 ∑ 2 prs 1 ∑ 2 pqs

n λp (1 − nrλ )(1 − nsλ ) t λ p λ r λ s λ

ε tλ + ε pλ − ε rλ − ε sλ n λp nqλ (1 − nsλ ) p λ q λ t λ s λ

ε pλ + ε qλ − ε tλ − ε sλ

Thus, the contribution from Part I to IP can be described as


2

2

(12)

Page 7 of 42

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


(I) IP ,k

λ λ λ λ λ λ λ 1 n j (1 − na )(1 − nb ) k j a b = ∑ 2 jab ε kλ + ε λj − ε aλ − ε bλ

1 ni n j (1 − nb ) i j k b ∑ ε iλ + ε λj − ε kλ − ε bλ 2 ijb λ

−

λ

λ

λ

λ

λ λ

2

(a) (13)

2

(b)

and that to EA as

(I) EA,c

λ λ λ λ λ λ λ 1 n j (1 − na )(1 − nb ) c j a b = ∑ 2 jab ε cλ + ε λj − ε aλ − ε bλ

1 ni n j (1 − nb ) i j c b ∑ ε iλ + ε λj − ε cλ − ε bλ 2 ijb λ

−

λ

λ

λ

λ

λ λ

2

(a) (14)

2

(b)

We would like to emphasize that the sums over occupied and virtual orbitals both include the fractionally occupied orbital k in Eq. 13 for the IP calculations or orbital c in Eq. 14 for the EA calculations as prescribed in Eqs. 1 and 2. Eq. 10 suggests that the orbital energy also depends explicitly on the occupation numbers: ∂ε λp = pλ t λ pλ t λ λ ∂nt

(15)

By taking this dependence into account, Part II of Eq. 11 can be derived as1,4

n λp nqλ (1 − nrλ )(1 − nsλ ) p λ q λ r λ s λ 1 (II) = − ∑ 2 4 pqrs (ε pλ + ε qλ − ε rλ − ε sλ )

2

(16)

⋅  t λ p λ t λ p λ + t λ q λ t λ q λ − t λ r λ t λ r λ − t λ s λ t λ s λ  . This result was first given in Ref. 4 and called the high-order terms (HOT) as they can be derived from the third-order energy diagrams (see below). Hence, the contribution from Part II to IP and EA can be formulated, respectively, as

(II) IP ,k

λ λ λ λ λ λ λ λ 1 ni n j (1 − na )(1 − nb ) i j a b =− ∑ 2 2 ijab (ε iλ + ε λj − ε aλ − ε bλ )

1 ni n j (1 − na )(1 − nb ) i j a b ∑ 2 2 ijab (ε iλ + ε λj − ε aλ − ε bλ ) λ

+

λ

λ

λ

λ

λ

λ λ


2

k λiλ k λiλ (17)

2

k λ aλ k λ aλ


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 42

and

(II) EA,c

λ λ λ λ λ λ λ λ 1 ni n j (1 − na )(1 − nb ) i j a b =− ∑ 2 2 ijab (ε iλ + ε λj − ε aλ − ε bλ )

1 ni n j (1 − na )(1 − nb ) i j a b ∑ 2 2 ijab (ε iλ + ε λj − ε aλ − ε bλ ) λ

+

λ

λ

λ

λ

λ

λ λ

2

cλ iλ cλ iλ (18)

2

cλ aλ cλ aλ

Approximations up to Part I and Part II have been coined as levels of MP2-I and MP2-II, respectively, in our previous publication.1 To realize the full theory for fractional occupation systems denoted as level of MP2-III,1 one has to derive the energy gradient from the consequence of orbital varying itself with occupation numbers, i.e. Part III in Eq. 11. This is accomplished as described below in the next subsection.

B. The full MP2 energy gradient for fractional occupation systems We expand the orbital variations with respect to changes of orbital occupation in the t-th orbital in terms of canonical molecular obitals {φqλ } :1 λ dφ pλ = ∑ φqλU qpnt λ dnt q

(19)

λ

where U nt is a rotation matrix that allows for orbital mixing as a result of orbital occupation change. The orthonormality of {φqλ } makes U nt an antisymmetric matrix, i.e. λ

λ S pq = p λ q λ = δ pq

λ λ λ dS pq λ = ∑ S rqλ U rpnt + S pr U rqnt λ dnt r

(

λ

(20)

)

(21)

λ

nt = U qpnt + U pq =0

which reduces nearly half of the number of matrix elements required for calculations. As frequently λ

shown, U nt should be obtained by solving the coupled-perturbed (CP) HF equations:1,15-23


Page 9 of 42

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


∑ ' A′λ

λ

λ

U bjnt = Baiλ ,t or A ′λ U nv-ot = B λv-o,t

(22)

ai ,bj

bj

The matrix elements are now extended to fractional occupation systems, defined by λ A′pqλ , rs = − ( ε rλ − ε sλ ) δ rpδ sq − Apq , rs

λ λ λ Apq , rs = ns (1 − nr )

(

(23)

pλ r λ qλ sλ + pλ sλ qλ r λ

)

(24)

and λ λ λ B pq qλ t λ ,t = p t

(25)

These equations are derived in the standard way,1,15-23 taking derivative of the Fock matrix Fpqλ with respect to the occupation number ntλ , while recognizing the fact that Fpqλ is diagonal,

Fpqλ = p λ hˆ q λ + ∑ nsλ p λ s λ q λ s λ = δ pqε pλ

(26)

s

leading to λ dFpqλ ntλ = ( ε pλ − ε qλ ) U pq + ∑ nλj (1 − nbλ ) p λ bλ q λ j λ + p λ j λ q λ b λ U bjnt + p λ t λ q λ t λ λ dnt jb

(

)

= 0 for p ≠ q

(27)

Like the previous ones for integer occupations,1 these equations are needed to calculate the λ

virtual-occupied (v-o) block U nv-ot . Hence, we set p = a and q = i in Eq. 27, which is just Eq. 22 with the matrix elements defined by Eq. 23‒25. Note that, however, as they are generalized to allow for fractional orbital occupations, there is a chance that both b and j stand for the same partially λ

occupied orbital, for which the matrix element U bjnt can be derived from Eq. 21 to be 0. Here, for convenience, the summation is marked with an apostrophe,

∑ ' , to exclude this element as in Eq. 22

and others shown below. There are other ways to extend the CPHF equations from integer to fractional systems. In particular, a different fractional expression has been proposed in the literature, which was based on the scaled orbitals φ%i → ni φi and φ%a →

(1 − na )φa .24-27 It can be shown that their rotation matrix



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 42

for the v-o block is related to ours as niλ (1 − naλ )U aint . λ

To obtain other blocks of the rotation matrix, more equations are needed. They are λ

U ljnt =

λ 1   ' Aljλ,aiU aint + Bljλ,t  λ ∑ ε j − ε l  ia 

(28)

λ

λ

for the occupied-occupied (o-o) block U no-ot , and λ

U bdnt =

λ 1   ' Abdλ ,aiU aint + Bbdλ ,t  λ ∑ ε d − ε b  ia 

λ

(29)

λ

for the virtual-virtual (v-v) block U nv-vt . In addition, the full derivatives of orbital energies can now be expressed as

d ε λp ntλ λ λ = ∑ ' App ,bjU bj + B pp ,t λ dnt bj λ

where Bpp ,t

(30)

∂ε λp . corresponds to Eq. 15 for ∂ntλ

By combining the above equations, Part III including orbital variation itself can be expressed as (III) = ∑ Pij( i≠ j

2 ),λ

MP 2, λ λ Bijλ,t + ∑ Pab( ) Bab U aint ,t + 2∑ ' Lai

λ

2 ,λ

a ≠b

(31)

ai

with 2,λ ,λ LMP = − 1 ∑ niλ n λj nlλ (1 − naλ )(1 − nbλ ) t ab i λ bλ j λ l λ ai jl 2 jlb

+ 1 ∑ niλ n λj (1 − naλ )(1 − nbλ )(1 − ndλ ) tijcb ,λ d λ b λ a λ j λ 2 jbd

(32)

+ ∑ niλ (1 − naλ ) Pjl(2),λ j λ a λ l λ i λ + ∑ niλ (1 − naλ ) Pbd(2),λ b λ a λ d λ i λ jl

bc

where tijab is the double substitution amplitude tijab ,λ = i λ j λ a λ b λ ( Here Pij

2 ), λ

(ε

λ i

+ ε λj − ε aλ − ε bλ )

(33)

and Pab( 2 ), λ are the elements of P ( 2 ),λ for the o-o block and the v-v block, respectively,

while P ( 2 ),λ is the MP2 correction matrix to the one-particle density matrix P λ , in which only the diagonal element can have a non-zero value given by the orbital occupation n λp :


Page 11 of 42

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


Ppqλ = δ pq n λp

(34)

The elements in the o-o and v-v blocks of P ( 2 ),λ are defined, respectively, as

Pij( 2),λ

 − 1 niλ nlλ (1 − naλ )(1 − nbλ ) t ab ,λ t ab ,λ for i = j il il  2∑ lab = ,λ − 1 ∑ niλ nλj nlλ (1 − naλ )(1 − nbλ ) tilab ,λ t ab for i ≠ j jl 2 lab 

(35)

Pab( )

 1 niλ n λj (1 − naλ )(1 − ndλ ) tijad ,λ tijad ,λ for a = b  2 ∑ ijd =  1 ∑ niλ n λj (1 − naλ )(1 − nbλ )(1 − ndλ ) tijad ,λ tijbd ,λ for a ≠ b  2 ijd

(36)

and

2 ,λ

While those in the o-o and v-v blocks can be obtained directly by definitions, elements in the v-o block of P ( 2),λ can only be solved via numerical solution of a modified Z-vector method19

∑ ' P ( ) λ A′λ 2 ,

bj

bj , ai

( ) 2, λ = LMP or Pv-o A ′λ = LMP 2,λ , ai 2 ,λ

(37)

bj

(

)

( ) ( 2 ),λ Meanwhile, the o-v block of P ( 2),λ can easily be acquired by Po-v = transpose Pv-o . Through 2 ,λ

CP-HF equations described in Eq. 22, the last term in Eq. 31 can be reformulated as ( 2 ),λ ( 2 ),λ λ 2LMP 2,λ ⋅ U nv-ot = Pv-o ⋅ B λv-o,t + Po-v ⋅ B o-v, t λ

(38)

Finally, Part II and Part III can be combined, yielding ( ) ( ) ( ) ( ) λ λ λ λ (II)+(III) = Po-o ⋅ B o-o, t + Pv-v ⋅ B v-v,t + Pv-o ⋅ B v-o,t + Po-v ⋅ B o-v,t 2 ,λ

2 ,λ

2 ,λ

2 ,λ

( II+III )IP,k = ∑ Pij( 2),λ

k λ i λ k λ j λ + ∑ Pab( )

k λ a λ k λ bλ

+ ∑ ' Pai( )

k λ a λ k λ i λ + ∑ ' Pia(

( II+III ) EA,c = ∑ Pij( 2),λ

c λ i λ c λ j λ + ∑ Pab( )

(39)

Hence, their contributions to IP and EA are, respectively,

ij

2 ,λ

ai

2 ,λ

ab

2 ),λ

k λiλ k λ aλ

(40)

ai

and

ij

+ ∑ ' Pai

( 2 ), λ

ai

2 ,λ

c λ a λ c λ bλ

ab

λ

λ

λ λ

c a c i

+ ∑ ' Pia(

2 ),λ

ai


cλ i λ cλ a λ

(41)


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

Now the full MP2 energy gradient for fractional occupation systems is established.

C. Interpretation with antisymmetrized Goldstone diagrams We use ASG diagrams for the interpretation of the derived equations for better understanding. Some conventional ASG rules are summarized as follows,14 with two additional rules being introduced here specifically for hole or particle lines with fractional occupation number n λp : (1) Any unlabeled line in an ASG diagram is summed over all possible hole or particle indices, where upward arrows correspond to particles in virtual orbitals and downward arrows to holes in occupied orbitals. (2) Each two-body vertex is interpreted as an antisymmetrized integral with the usual arrangement of . (3) Denominators are derived from the generating energy diagrams, while the resolvent lines are not shown by convention. An energy denominator is included between each successive pair of vertices, calculated as the sum of the orbital energies of the hole lines minus those of the particle lines all crossed by the same resolvent line. (4) A phase factor

( − 1)

h +l

is assigned, where h and l are the numbers of the hole lines and the

closed loops, respectively. (5) A weight factor (1 / 2 )

m

is assigned, where m is the number of equivalent line pairs, while

two lines are said to be equivalent if they connect the same two vertices in the same direction. (6) A factor (1 − n λp ) is assigned for each

, where the particle line with occupation number

n λp connects two vertices; while a factor n λp is assigned for each

occupation number n λp connects two vertices;


, where the hole line with

Page 12 of 42

Page 13 of 42

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


(7) A factor 1 (1 − n λp ) is assigned for each bracket

, while a factor 1 n λp for each

. Here the right

indicates that the two lines connected have the same line indexes, while the former and

the latter correspond to particle and hole lines, respectively, with occupation number n λp . Applying Rules 1-6 to Eqs. 7, 13 and 14, we obtain the corresponding ASG diagrams as shown in Figures 1. Diagrams on the left of the horizontal arrows correspond to the MP2 correlation energy terms (Eq. 7). The differentiation with respect to the occupation number corresponds to a systematic opening of the energy diagrams, where the scissors indicate which line is to cut. Thus, diagrams on the right of the horizontal arrows correspond to the derivatives at the MP2-I level that considers only the explicit dependence of the MP2 correlation energy on orbital occupations (Eqs. 13 and 14). Note that,4 if a particle line is cut, the line direction in the IP term shall be changed to express the correspondence between Eq. 13a and Figure 1a, since the open line has to have the hole character to describe an ionization. Similarly, if a hole line is cut, the line direction in the EA term shall be changed, i.e., the correspondence between Eq. 14b and Figure 1b, such that the open line has the particle character. The diagrammatic expressions for the elements of the MP2 correction matrix P (

2 ), λ

are shown in

( ) Figure 2. While Figure 2a corresponds to the o-o block, Po-o , defined in Eq. 35, Figure 2b 2 ,λ

( ) corresponds to the v-v block, Pv-v , defined in Eq. 38. As compared to the off-diagonal elements in 2 ,λ

( ) ( ) Eqs. 35 and 38, the diagonal elements of Po-o and Pv-v 2 ,λ

2 ,λ

drop a factor of nλj and (1 − nbλ ) ,

respectively, which can be easily explained through their diagrammatic expressions in Figure 2 by additionally applying Rule 7. Now it becomes clear that Part II of the energy gradient, i.e. Eq. 17 and 18 as developed originally at the MP2-II level,4 contains only the diagonal contribution, while the off-diagonal contribution comes from Part III (Eq. 31). By combining Part II and the first two terms



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 14 of 42

in Part III, we arrive at the full description of the o-o and v-v contributions to IP and EA as illustrated by diagrams in Figure 3. Similar to those noted by Beste and co-workers as the high-order terms (HOT),4 these diagrams can be viewed as deriving from a subset of third-order energy terms by cutting the bubbles as illustrated in Figure 3. Although solving Eq. 37 provides numerical solution of the v-o block of P ( 2 ),λ needed for full MP2 energy gradients, it does not render us more clear understanding of the elements in this block. For better understanding, we again resort to diagrams. As we shall see below, higher up to infinity order of energy terms are included via the CP equations. With some adjustment, Eq. 37 is transformed to Pai(

2 ), λ

MP 2, λ Abjλ , ai 2 ,λ = Lλai λ + ∑ ' Pbj( ) λ ε i − ε a jb ε i − ε aλ

(42)

For its validity and without loss of generality, here and hereinafter, it is required that i and a represent the conventional occupied and virtual orbitals, respectively, and should not run into the same orbital index as a partially occupied orbital. Both sides of Eq. 42 contain the v-o block elements, such that this set of equations should be solved iteratively. The process starts by inserting the initial value of Pai( ,()0 ) = 0 2 ,λ

into the right hand side of Eq.

2, λ Pai( 2,()1,)λ = LMP ai

(ε

λ i

42 to

obtain the first iteration result

− ε aλ ) , with which the contributions to IP or EA from Pai( 2,()1,)λ can be formulated as


Page 15 of 42

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 ), λ λ λ Pv-o, (1) ⋅ B v-o,t + Po-v,(1) ⋅ B o-v,t

j λ l λ a λ bλ iλ bλ j λ l λ a λ t λ i λ t λ 1 λ λ λ λ λ = − ∑ ' ∑ ni n j nl (1 − na )(1 − nb ) 2 ai jlb (ε λj + ε lλ − ε aλ − ε bλ )(ε iλ − ε aλ )

(a )

i λ j λ d λ bλ d λ bλ a λ j λ a λ t λ i λ t λ 1 λ λ λ λ λ + ∑ ' ∑ ni n j (1 − na )(1 − nb )(1 − nd ) 2 ai jbd (ε iλ + ε λj − ε dλ − ε bλ )(ε iλ − ε aλ ) −

a λ bλ j λ l λ j λ l λ i λ bλ i λ t λ a λ t λ 1 λ λ λ λ λ ' 1 − 1 − n n n n n ∑ ∑ i j l ( a )( b ) (ε λj + ε lλ − ε aλ − ε bλ )(ε iλ − ε aλ ) 2 ai jlb

(c )

d λ bλ i λ j λ a λ j λ d λ bλ i λ t λ a λ t λ 1 λ λ λ λ λ + ∑ ' ∑ ni n j (1 − na )(1 − nb )(1 − nd ) 2 ai jbd (ε iλ + ε λj − ε dλ − ε bλ )(ε iλ − ε aλ ) (2), λ jl

+ ∑ ' ∑ ni (1 − na ) P

 j λ aλ l λ iλ   ε iλ − ε aλ 

(2),λ bd

+ ∑ ' ∑ ni (1 − na ) P

 bλ a λ d λ iλ   ε iλ − ε aλ 

λ

ai

jl

λ

ai

λ

λ

bd

λ λ

λ λ

+

i t a t

λ λ

λ λ

i t a t

+

j λiλ l λ aλ

ε iλ − ε aλ bλ iλ d λ a λ λ

λ

εi − εa

(b )

 aλ t λ iλ t λ     a t i t    λ λ

λ λ

(d )

(e ) (43)

(f )

Its diagrammatic expressions are illustrated in Figure 4, where diagrams (a)-(f) correspond to terms (a)-(f) of Eq. 43 respectively. The process continues by taking the result from the previous iteration as the new stating point: Pai( ,()n ) = ∑ ' Pbj( ,()n −1) 2 ,λ

2 ,λ

jb

Abjλ , ai ε iλ − ε aλ

To interpret equations as diagrams, it is recognized that being multiplied by

equation is equivalent to adding

and

(44)

Abjλ , ai in the ε iλ − ε aλ

to the original diagram. Hence, the contributions

( 2 ), λ ( 2 ), λ λ λ from the v-o and o-v blocks to IP or EA in the n-th iteration, i.e. Pv-o, ( n ) ⋅ B v-o,t + Po-v,( n ) ⋅ B o-v,t , can be

expressed diagrammatically as in Figure 5. Eventually, Pai(

2 ), λ

can be obtained by adding term by

term up to infinity, i.e. Pai( 2 ), λ = Pai( 2,()1,)λ + Pai( 2,()2,λ) + Pai( 2,()3, λ) + L .

(45)

By defining the dressed interaction as a wavy line shown in Figure 6, the contributions to IP and EA



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

from Parts II and III, due to the inexplicit dependence on the occupation variations, can now be compactly expressed with six diagrams as shown in Figure 7. They resemble the standard third-order energy diagrams but with dressed interactions up to the infinity order.

III. COMPUTATIONAL DETAILS All calculations were performed in the spin-unrestricted formalisms using a modified version of NorthWest computational Chemistry (NWChem) software package, version 6.1.1.28 For atomic calculations, no requirement for spherical symmetry was imposed on electron densities. For benchmarking purpose, the variable, tol2e in NWChem,28 was set very tightly to be 1E-16. This determined the integral screening threshold for the evaluation of the energy and related Fock-like matrices. The SCF energies were always converged to be less than 1E-11 au. Calculations have been performed on Set 1, a set of atoms (i.e., Li, Be, B, C, N, O and F), as well as some molecules (i.e., F2, OH, NH2, CH3, CN and O2), where the frontier orbital energies of HOMO (highest occupied molecular orbital) and LUMO (lowest unoccupied molecular orbital) were compared with the reference data of vertical IPs and EAs. The basis set used is cc-pVQZ,29,30 with the spherical angular functions. The reference data were generally the experimental values taken from Refs. 31, and 32, while for EAs of Be, N, and molecules, the CCSD(T) data from Ref. 33 were used as references, as no vertical data from experiments are available. This data set for frontier orbitals was compiled by Cohen, Mori-Sánchez and Yang,34 which has been used in discussing their fractional charge perspective on the band gap problem in density functional theory,34 as well as MP2.3 The same set has also been used to examine the fractional charge behavior of some doubly hybrid functionals.35 Here in this work, the reference data for IP of NH2 has been revised from 11.14 to 12.00 eV, as the former actually corresponds to the adiabatic IP,31,36 while the latter is a vertical


Page 16 of 42

Page 17 of 42

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


value.31,37 Also, the reference data for both IP and EA of CN have been revised due to lack of vertical values from experiments. The reference data now are from CCSD(T) extrapolated to the complete basis set limits in the same way as in Ref. 33. Calculations have also been performed on Set 2, a set of molecules (i.e., H2O, CH2O, CH4, NH3, and N2), where not only the valence orbital energies, but also the core orbital energies were compared with the experimental data of vertical IPs.4,38,39 This data set for vertical IPs of core & valence orbitals was compiled by Beste, Vázquez-Mayagoitia and Ortiz,4 which has been used in discussing their direct ∆MBPT(2) method that used the energy gradients at the MP2-II level.4 Following their work,4 we calculated the vertical IPs using a completely uncontracted cc-pVTZ basis set with the spherical angular functions. Molecular geometries for both Set 1 and Set 2 were consistently optimized at the MP2/cc-pVTZ level, where the results are provided as the supporting information. Previously,1 the experimental geometries were adopted for Set 1 and the Cartesian angular functions were in company with the cc-pVQZ basis set.1,3,35 This choice is found to generally change the statistics for the whole set at the second decimal point, although some calculated EAs, e.g. ∆MP2 for F2, are prone to geometry changes, where the F-F distance from MP2/cc-pVTZ is 1.3958 Å, as compared to the experimental one of 1.4119 Å.40

IV. RESULTS AND DISSCUSION A. The MP2 Energy Derivatives at Different Levels of Approximations

Deviations for the calculated MP2 energy derivatives with respect to the orbital occupation numbers ntλ are depicted in Figure 8 for (a) 1s, (b) 2s, (c) 2p, and (d) 3s of Ne atom. Here 2p corresponds to



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

HOMO and 3s to LUMO. MP2 numerical energy derivatives, ∆EMP 2 ∆nt , are taken as the references, where ∆nt was chosen as 0.001.1,3,35 Data points were actually recorded in intervals of 0.1 electrons (See supporting information for details), although smooth lines are drawn for eye guidance only. As shown in Figure 8, values from the MP2-I level, where only the explicit dependence of the MP2 energy on the occupation number is considered, deviate from the reference data most significantly for all orbitals in Ne atom. The deviations at the MP2-I level for 1s, 2s, and 2p orbitals are up to 1.67, 0.70, and 0.67 eV, respectively. This suggests that higher-order terms that take into account the inexplicit orbital dependence on occupation numbers are important. For the example of Ne atom shown here, the MP2-II level is always an improvement over MP2-I. The largest deviations at the MP2-II level for 1s, 2s, and 2p orbitals are reduced to 0.91, 0.11, and 0.34 eV, respectively. Note that, the MP2-II correction can overshoot, showing positive deviations. The MP2-III level is actually the full theory, whose gradients overlap with the numeric references for all orbitals under investigation. Previously, we have shown,1 for a set of atoms and molecules,3,27,28 that the MP2-III level using the analytical gradient formalisms at integer numbers of λ = 0 and 1 faithfully reproduce the corresponding numeric gradient results. Figure 8 demonstrates that the analytical gradients are correctly derived and implemented as the full MP2 theory for the fractional occupation systems (See supporting information for details). From Figure 8 it is seen that the deviations at the MP2-I and MP2-II levels are usually higher at the fractional points than at the integer points, highlighting the importance to explore the full theory for the fractional occupation systems. We will now apply this full MP2 theory for IP and EA calculations, discussing no more results


Page 18 of 42

Page 19 of 42

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


from approximations at the MP2-I and MP2-II levels.

B. Ionization Potentials and Electron Affinities for Set 1 [ →1] [ →1] Previously,1 we have used the symbols, IPHF ( 0 ) and IPHF (1) , to indicate that IPs are calculated at the

HF level, where the energy gradients with respect to orbital occupations are evaluated at λ = 0 and 1, respectively, whose results are exact up to the first order [→1] in terms of coupling constant λ. [ → 2] Accordingly, the two-point equation for IP calculations at λ = (0, 1) is symbolized as IPHF ( 0,1) and

calculated according to4 dE HF ( λ )  1  dE HF ( λ ) [ → 2] IPHF +   (0,1) = − λ dnkλ 2  dnk  λ =0 λ =1 

(46)

The result is exact up to the second order as indicated by the superscript [→2]. Here in this work, Slater’s transition state method will be used,6,7 which corresponds to the [ → 2] one-point equation with energy gradients evaluated at λ = 1/2, i.e., IPHF (1/ 2 ) . An extended Slater’s

transition state method will also be used here,8 which is a two-point equation at λ = (0, 2/3), being exact up to the third order:8,9  dE HF ( λ ) 1  dE HF ( λ ) [ → 3] IPHF + 3⋅   (0,2/3) = − λ λ 4  dnk dnk  λ =0 λ = 2/3 

(47)

The HF method can be updated to the MP2 method using formulae for analytical gradients developed in the previous sections. The calculated IP results for Set 1 are summarized in Table 1 with the HF and MP2 methods using various approximations. Note that a short hand notation is 2 actually used in the table, where, for example, IP( MP 0,2 /3) is used to stand for a two-point IP calculation

at λ = (0, 2/3), which is similar to Eq. 47 but at the MP2-III level of the full MP2 theory, being exact MP 2 III ,[ → 3] up to the third order. That is, IPMP 2 (0,2 /3) is simplified as IP( 0,2 /3) .



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

Data in Table 1, in the column under the title of IP( 0HF) , correspond to the well-known Koopmans’ 0 theorem,41 which are just the opposite of the orbital energies for the HOMOs, i.e., IP( 0HF) = − ε HOMO .

Due to the lack of orbital relaxation effects, IP( 0HF) yields numbers which are consistently too high as compared to the ∆HF results. This can also be related to the concave errors associated with the HF method from the perspective of fractional charge.3,34,35,42,43 The orbital relaxation effect can be taken into account via the two-point equation, e.g., at λ = (0, 1). Hence IP(0HF,1) leads to a significant improvement. We calculate the mean absolute deviation (MAD_1) with respect to the corresponding ∆ method results for Set 1. As shown in Table 1, MAD_1 from IP( 0HF) is as high as 1.59 eV, which is drastically reduced to 0.15 eV by using IP(0HF,1) . Analogously, the maximum deviation (MAX_1, Ref. ‒ Calc.) with respect to the ∆HF results is ‒2.82 eV for IP( 0HF) , which drops to only ‒0.31 eV for IP(0HF,1) . Due to the lack of correlation effects, ∆HF leads to MAD_2 of 0.98 eV with respect to the

experimental reference data, while IP( 0HF) and IP(0HF,1) give MAD_2 of 1.00 and 0.85 eV, respectively. Table 1 shows that both IP(1/HF2 ) and IP( HF 0,2 /3) are able to take into account the relaxation effects properly. They faithfully reproduce the ∆HF results, yielding MAD_1 of 0.08 and 0.01 eV for IP(1/HF2 ) HF and IP(0,2 / 3) , respectively.

It is very interesting to compare the MP2 results with the HF ones when correlation effects are accounted for. IP( 0MP) 2 gives IPs which are consistently smaller than those from ∆MP2. This is an indication that the MP2 method has a convex error as was found before for conventional DFT methods.13,34,35,42-44 MAD_1 from IP( 0MP) 2 is 0.65 eV, which is less than half of that from IP( 0HF) , indicating that the magnitude of the delocalization error for the former is much smaller than the localization error for the latter. Although both being exact up to the second order, the values


Page 20 of 42

Page 21 of 42

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 predicted at λ = (1/2) is always better than that at λ = (0, 1) for either HF or MP2 for Set 1. IP( MP 0,2 /3)

reproduces even more satisfactorily the ∆MP2 results, giving MAD_1 of 0.10 eV. With respect to the experimental reference data, MP2 is a significant improvement over HF, showing the importance of the correlation effect. On average, ∆MP2 gives MAD_2 of only 0.15 eV, which significantly improves over MAD_2 of HF (0.98 eV). However, IP( 0MP) 2 leads to MAD_2 of MP 2 0.77 eV, which is unacceptably high, suffering from the convex error of MP2. IP(0,1) , IP(1/MP2)2 , and

2 MP 2 , yielding MAD_2 of 0.40, 0.19 and 0.20 eV, IP( MP 0,2 /3) overcome this difficulty encountered by IP( 0 )

respectively. The latter two are close to the ∆MP2 results. Correlation effects may change the interpretation. For example, HF leads to an ionization from the π orbital of CN, predicting a 3Π as the ground state of CN+, whereas MP2 predicts an ionization from the σ orbital of CN, suggesting a 1Σ+ as the ground state of CN+. Although the MP2 assignment is supported by the extensive configuration interaction (CI) calculations,45 one has to note that the unrestricted MP2 method suffers from the spin contamination and have difficulty in treating species of multi-reference characters. The EA results calculated with the HF and MP2 methods using various approximations are presented in Table 2. As is well known, EA(HF that uses LUMO orbital energy to approximate EA is 0) in serious error. MAD_2 from EA(HF is 2.92 eV, which is too high to be useful in predicting an 0) experimental EA, while MAD_1 from EA(HF is also as high as 1.39 eV, being consistently too large 0) HF HF in magnitude as a result of localization errors. EA(0,1) , EA(1/HF2) , and EA(0,2 give continuously /3)

improved approximation to ∆HF with MAD_1 of 0.24, 0.11 and 0.02 eV, respectively. Nevertheless, ∆HF itself is not a useful theory for EA prediction, giving MAD_2 of 1.53 eV, and erroneously

predicting many anion states (e.g. Li, B, O, F2, OH, NH2, CH3 in Table 2) unstable.



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

Correlation effects also play an important role in EA calculations. ∆MP2 corrects the qualitative errors associated with ∆HF, reducing MAD_2 from 1.53 for ∆HF to 0.45 eV. As shown by data in 2 Table 2, EA(MP seems to provide a good approximation to ∆MP2, where the associated error, 0)

2 MAD_1, is only 0.28 eV. Errors become more pronounced for EA(MP 0 ,1) , where the associated error,

MAD_1, is increased to 0.47 eV. Results at λ = (1/2) are always better than those from λ = (0, 1), while those from λ = (0, 2/3) are always more satisfactory than those from λ = (1/2). Hence, MAD_1 MP 2 from EA(1/MP2 2) , and EA(0,2 /3) are 0.21 and 0.07 eV, respectively, increasingly approaching the results

from ∆MP2. We conclude that higher order contributions are treated more accurately using the transition state schemes.

C. Vertical Ionization Potentials for Core and Valence Orbitals in Set 2

Table 3 summarizes the vertical IP results for Set 2 calculated with the HF and MP2 methods using various approximations. As shown in Table 3, the relaxation effects are most significant when ionizations occur at the core orbitals, while the relaxation effects for the valence orbitals are similar to those for the HOMOs discussed in the previous subsection. For example, IP( 0HF) overestimates the core IP by 20.43 eV in H2O, while the errors for its other valence orbitals including HOMO are all around 2.5 eV. MAD_1c, i.e., the error statistics for 7 core orbitals, from IP( 0HF) for Set 2 is 14.86 eV, while MAD_1v, i.e., the error statistics for 19 valence orbitals, from IP( 0HF) for Set 2 is 2.08 eV. HF HF Relaxation effects can be properly taken into account by IP(0,1) , IP(1/HF2 ) , and IP(0,2 / 3) , leading to results

that approach to those from ∆HF with similar MADs for both core and valence electrons. MAD_1c HF HF or MAD_1v for IP( HF 0,1) , IP(1/ 2 ) , and IP( 0,2 /3) are around 0.17, 0.08 and 0.02 eV, respectively.

With respect to the experimental reference data,4,38,39 ∆HF gives MAD_2c of 3.26 eV and


Page 22 of 42

Page 23 of 42

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


MAD_2v of 1.36 for core and valence electrons, respectively. IP( 0HF) deteriorates the situation for core electrons, yielding MAD_2c of 16.96 eV, due to lack of the relaxation effects. It is, therefore, always faulty to use HF orbital energies, i.e. IP( 0HF) , to approximate the experimental vertical IPs for core ionizations, emphasizing the importance of taking into account of the relaxation effects. For valence electrons, IP( 0HF) even gives a smaller MAD_2v of 1.09 eV for this set, due to error HF cancellation, which is, nonetheless, not guaranteed for every case. IP(0,1) , IP(1/HF2 ) , and IP( HF 0,2 /3) can better

reproduce the ∆HF results, giving MAD_2c/MAD_2v of 3.42/1.22, 3.19/1.42 and 3.27/1.34 eV, respectively. Clearly, correlation effects are important in reproducing the experimental data. MP2 is the simplest wavefunction method that includes the correlation effects. ∆MP2 leads to important improvement over ∆HF, reducing MAD_2c and MAD_2v to 1.68 and 0.31 eV, respectively. The MP2 results are much more satisfactory for the valence IPs than for the core IPs. It is widely recognized that,46 due to the lack of correlation effects, HF leads to a wrong order of the valence orbitals of N2, predicting that 1eu is higher lying than 3a1g. Inclusion of the correlation effects allows for the correct description of the ionization spectrum of N2, showing that 1eu is actually more than 1 eV lower lying than 3a1g. As shown in Table 3, IP( 0MP) 2 is not a good approximate to ∆MP2 even for valence IPs, where MP 2 2 MAD_1v is 1.69 eV. IP(0,1) , IP(1/MP2 )2 , and IP( MP 0,2 /3) provide increasingly better approximation to ∆MP2,

yielding MAD_1v of 0.56, 0.24 and 0.13 eV, respectively. However, the core IPs can deteriorate the statistics. Taking the innermost core IP of N2 as an example, IP( 0MP) 2 deviates from the ∆MP2 MP 2 reference by ‒3.62 eV, while the deviation is reduced to ‒0.81 eV with IP(0,1) . Nonetheless, further

error reduction is not seen at IP(1/MP2 )2 . Instead, IP(1/MP2 )2 leads to an increased error up to 1.65 eV, 2 although the deviation is eventually reduced to ‒0.02 eV with IP( MP 0,2 /3) . The abnormal behavior at



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

IP(1/MP2 )2 is worthy of note. Even though the total energy seems to change smoothly as a function of the

occupation numbers for the innermost core orbital as shown in Figure 9(a), the deviations from linearity are clearly shown in Figure 9(b). In particular, Figure 9(c) depicts the energy derivatives as a function of the occupation number, where the abnormality is clearly seen in between occupation numbers of (0.4, 0.7). The reason for such erratic behaviors at certain values of ntλ at the MP2 level has been discussed by Beste, et al.4 If the negative orbital energy for an inner orbital k is to be subtracted when it takes the role of an unoccupied orbital, the denominator in energy derivatives can become very small or even tend to zero. Such an effect is particularly large in Part II in calculating the diagonal elements of the MP2 correction matrix P ( 2),λ , as the sum of the orbital energies in the denominator is squared.4

IV. CONCLUDING REMARKS In the present work, we have derived the formalism to calculate the full MP2 energy derivatives with respect to the orbital occupation numbers. The present formalism is applicable not only to the integer points for λ = 0 and λ = 1 as described in our previous work,1 but also to the systems of fractional occupations. Antisymmetrized Goldstone (ASG) diagrams have been employed for interpretations and better understanding of the derived equations, where two additional rules have been introduced in the present work specifically for hole or particle lines with fractional occupation numbers. In addition to the explicit dependence of energy on the occupation variations, the inexplicit dependence has been compactly expressed with six diagrams, which resemble the standard third-order energy diagrams but with dressed interactions up to the infinity order. The formalism has been implemented and solved through the CP-HF equations, which has been combined with the Slater’s transition state


Page 24 of 42

Page 25 of 42

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


concept to have accurate predictions of IPs and EAs at the MP2 level. The present study can also be extended to compute the excitation energies,4 Fukui function and response function for nonlocal and fractional systems.25 Extension to the recently developed doubly hybrid functionals can also be envisioned.

ACKNOWLEDGEMENTS

This work was supported by National Natural Science Foundation of China (21133004, 91427301), and the Ministry of Science and Technology (2013CB834606).

ASSOCIATED CONTENT Supporting Information

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.xxxxxxx.

AUTHOR INFORMATION Corresponding Author

*Email: [email protected]. Tel.: +86 21 65643529 Notes

The authors declare no competing financial interest.

REFERENCES

1. Su, N. Q.; Xu, X. J. Chem. Theory Comput. 2015, 11, 4677–4688. 2. Møller, Chr., Plesset, M. S. Phys. Rev. 1934, 46, 618–622. 3. Cohen, A. J.; Mori-Sánchez, P.; Yang, W. J. Chem. Theory Comput. 2009, 5, 786–792. 4. Beste, A.; Vázquez-Mayagoitia, Á.; Ortiz, J. V. J. Chem. Phys. 2013, 138, 074101.



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

5. Beste, A.; Harrison, R. J.; Yanai, T. J. Chem. Phys. 2006, 125, 074101. 6. Slater, J. C.; Mann, J. B.; Wilson, T. H.; Wood, J. H. Phys. Rev. 1969, 184, 672–694. 7. Slater, J. C. Advan. Quantum. Chem.1972, 6, 1–92. 8. Williams, A. R.; DeGroot, R. A.; Sommers, C. B. J. Chem. Phys. 1975, 63, 628–631. 9. Chong, D. P. Chem. Phys. Lett. 1995, 232, 486–490. 10. Zhang, Y.; Xu, X.; Goddard, W. A., III Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 4963–4968. 11. Zhang, I. Y.; Xu, X.; Jung, Y.; Goddard, W. A., III Proc. Natl. Acad. Sci. U.S.A 2011, 108, 19896–19900. 12. Zhang, I. Y.; Xu, X. Int. Rev. Phys. Chem. 2011, 30, 115–160. 13. Su, N. Q.; Xu, X. Int. J Quantum Chem. 2015, 115, 589–595. 14. Shavitt, I.; Bartlett, R. J. Many-Body Methods in Chemistry and Physics, MBPT and Coupled-Cluster Theory; Cambridge University Press: New York, 2009. 15. Pople, J. A.; Krishnan, R.; Schlegel, H. B.; Binkley, J. S. Int. J. Quantum Chem. 1979, 16(S13), 225–241. 16. Frisch, M. J.; Head-Gordon, M.; Pople, J. A. Chem. Phys. Lett. 1990, 166, 275–280. 17. Salter, E. A.; Trucks, G. W.; Fitzgerald, G.; Bartlett, R. J. Chem. Phys. Lett. 1987, 141, 61–70. 18. Cioslowski, J.; Ortiz, J. V. J. Chem. Phys. 1992, 96, 8379–8389. 19. Handy, N. C.; Schaefer, H. F., III J. Chem. Phys. 1984, 81, 5031–5033. 20. Yamaguchi, Y.; Osamura, Y.; Goddard, J. D.; Schaefer, III, H. F. A New Dimension To Quantum Chemistry: Analytic Derivative Methods in Ab Initio Molecular Electronic Structure Theory; Oxford Univ. Press: Oxford, 1994. 21. Su, N. Q.; Zhang, I. Y.; Xu, X. J. Comput. Chem. 2013, 34, 1759–1774.


Page 26 of 42

Page 27 of 42

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


22. Su, N. Q.; Xu, X. Sci. Sin. Chim. 2013, 12, 1761–1779. 23. Su, N. Q.; Adamo, C.; Xu, X. J. Chem. Phys. 2013, 139, 174106. 24. Mori-Sánchez, P.; Cohen, A. J.; Yang, W. Phys.Rev.A 2012, 85, 042507. 25. Peng, D.; Yang, W. J. Chem. Phys. 2013, 138, 184108. 26. Steinmann, S. N.; Yang, W. J. Chem. Phys. 2013, 139, 074107. 27. Yang, W.; Mori-Sánchez, P.; Cohen, A. J. J. Chem. Phys. 2013, 139, 104114. 28. Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; van Dam, H. J. J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T. L.; de Jong, W. A. Comput. Phys. Commun. 2010, 181, 1477–1489. 29. Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007–1023. 30. Woon, D. E.; Dunning, T. H., Jr. J. Chem. Phys. 1993, 98, 1358–1371. 31. Linstrom P. J.; Mallard W. G., NIST chemistry WebBook, National Institute of Standards and Technology, Gaithersburg MD, 20899, (http://webbook.nist.gov). 32. Johnson R. D., III NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101, release 15b. http://cccbdb.nist.gov/ (August 2011). 33. Lin, Y. S.; Tsai, C. W.; Li, G. D.; Chai, J. D. J. Chem. Phys. 2012, 136, 154109. 34. Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Phys. Rev. B 2008, 77, 115123. 35. Su, N. Q.; Yang, W.; Mori-Sánchez, P.; Xu, X. J. Phys. Chem. A 2014, 118, 9201–9211. 36. Gibson, S.; Greene, J.; Berkowitz, J. J. Chem. Phys., 1985, 83, 4319–4328. 37. Locht, R.; Servais, C.; Ligot, M.; Derwa, F.; Momigny, J. Chem. Phys., 1988, 123, 443–454. 38. Jolly, W. L.; Bomben, K. D.; Eyermann, C. J. At. Data Nucl. Data Tables 1984, 31, 433–493. 39. Kimura, K.; Katsumata, S.; Achiba, Y.; Yamazaki, T.; Iwata, S. Handbook of HeI Photoelectron



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

Spectra of Fundamental Organic Molecules. Ionization Energies, Ab Initio Assignments, and Valence Electronic Structure for 200 Molecules; Halsted: New York, 1981. 40. Huber, K. P.; Herzberg, G. Constants of Diatomic Molecules, Molecular Spectra and Molecular Structure, Vol. 4; Van Nostrand Reinhold: Princeton, 1979. 41. Koopmans, T. Physica 1933, 1, 104–113. 42. Cohen, A. J.; Mori-Sánchez, P.; Yang, W. J. Chem. Phys. 2007, 126, 191109. 43. Vydrov, O. A.; Scuseria, G. E.; Perdew, J. P. J. Chem. Phys. 2007, 126, 154109. 44. Su, N. Q.; Xu, X. J. Phys. Chem. A 2015, 119, 1590–1599. 45. Hirst, D. M. Mol. Phys. 1994, 82, 359–368. 46. Cederbaum, L. S. Chem. Phys. Letters 1974, 25, 562–563.


Page 28 of 42

Page 29 of 42

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


TOC

Second-order perturbation theory for fractional occupation systems: Applications to ionization potential and electron affinity calculations

Neil Qiang Su, Xin Xu* Collaborative Innovation Center of Chemistry for Energy Materials, Shanghai Key Laboratory of Molecular Catalysis and Innovative Materials, MOE Laboratory for Computational Physical Science, Department of Chemistry, Fudan University, Shanghai, 200433, China



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 30 of 42

Table 1. Comparison with results from the experiments and the ∆ methods: IP from HOMO as calculated using the one-point and two-point equations.a ∆HF-IP

b

IP( HF 0)

c

IP( HF 0 ,1)

d

HF IP(1/2 )

e

IP( 0HF,2/3)

f

∆MP2-IP b

2 IP( MP 0)

c

2 IP( MP 0,1)

d

MP 2 IP(1/2 )

e

2 IP( MP 0 ,2/3)

f

Ref.g

Li

5.34

5.34

5.34

5.34

5.34

5.36

5.36

5.36

5.36

5.36

5.39

Be

8.04

8.42

8.10

8.02

8.05

8.87

8.67

8.85

8.87

8.85

9.32

B

8.04

8.67

8.09

8.02

8.04

8.31

8.18

8.26

8.33

8.30

8.30

C

10.80

11.94

10.90

10.75

10.80

11.30

11.11

11.18

11.35

11.28

11.26

N

13.89

15.52

14.02

13.83

13.90

14.63

14.27

14.42

14.73

14.60

14.54

O

12.02

14.19

12.28

11.89

12.04

13.41

12.93

13.03

13.59

13.35

13.61

F

15.65

18.47

15.94

15.51

15.67

17.36

16.36

16.79

17.63

17.29

17.42

F2

16.04

18.03

16.15

15.98

16.04

15.42

13.53

15.04

15.61

15.39

15.70

OH

11.37

13.95

11.67

11.22

11.39

13.06

12.02

12.49

13.33

12.98

13.01

NH2

10.47

12.62

10.71

10.36

10.49

12.00

11.19

11.60

12.19

11.95

12.00

CH3

8.99

10.46

9.12

8.93

8.99

9.76

9.23

9.57

9.86

9.75

9.84

CN

12.73

14.71

12.99

12.60

12.76

13.57

13.48

14.39

14.26

14.41

13.79

O2

13.52

15.31

13.61

13.47

13.52

11.82

10.14

11.53

11.96

11.80

12.30

MAD_1h

-

1.59

0.15

0.08

0.01

-

0.65

0.31

0.17

0.10

-

MAX_1h

-

‒2.82

‒0.31

0.15

‒0.03

-

1.89

-0.81

-0.69

-0.83

-

i

MAD_2

0.98

1.00

0.85

1.04

0.97

0.15

0.77

0.40

0.19

0.20

-

i

1.77

‒3.01

1.48

1.91

1.75

0.48

2.17

0.77

-0.47

-0.62

-

MAX_2 a

The data set is taken from Refs. 3, 34, 35. The unit is eV. The cc-pVQZ basis set is used for all calculations.

b

In the ∆ method, IP is calculated as ∆Emethod = Emethod(N0–1) – Emethod(N0), where method = HF or MP2.

c

IP is calculated up to the first order using one-point equation −

d

e

f

g

dEmethod ( λ ) at λ = 0. dnkλ

dEmethod ( λ )  1  dEmethod ( λ ) IP is calculated up to the second order using two-point equation −  + .  2  dnkλ dnkλ λ =0 λ =1  IP is calculated up to the second order using one-point equation −

dEmethod ( λ ) at λ = 1/2. dnkλ

 dEmethod ( λ ) 1  dEmethod ( λ ) IP is calculated up to the third order using two-point equation −  + 3⋅ . λ  4  dnkλ dn k λ =0 λ = 2/3  The experimental data are taken from Refs. 31, 32. Vertical IP for CN is CCSD(T) value calculated as in Ref. 33. For HF related methods, CN+ is in the triplet.

For MP2 related methods, CN+ is in the singlet. h

MAD_1 and MAX_1 refer to the mean absolute deviation and the maximum deviation (Ref. – Calc.), respectively, with respect to the ∆ method results.

i

MAD_2 and MAX_2 refer to the mean absolute deviation and the maximum deviation (Ref. – Calc.), respectively, with respect to the experimental results.


Page 31 of 42

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


Table 2. Comparison with results from the experiments and the ∆ methods: EAs from LUMO as calculated using the one-point and two-point equations.a ∆HF-EA b

EA(HF 0)

Li

-0.21

-0.38

-0.11

-0.26

-0.23

0.29

0.24

0.24

0.28

0.30

0.62

Be

-0.95

-1.22

-0.92

-0.97

-0.95

-0.78

-0.79

-0.81

-0.77

-0.78

-0.36

B

-0.47

-1.14

-0.34

-0.53

-0.48

0.00

-0.13

-0.12

0.05

0.02

0.28

C

0.28

-0.82

0.46

0.20

0.27

1.03

0.87

0.78

1.15

1.06

1.26

N

-2.34

-3.52

-2.03

-2.49

-2.37

-0.97

-1.36

-1.37

-0.82

-0.89

-0.22

O

-0.94

-2.73

-0.55

-1.12

-0.97

0.86

0.55

0.14

1.18

0.97

1.46

F

0.84

-1.61

1.26

0.63

0.81

3.04

3.00

2.04

3.52

3.17

3.40

F2

-0.35

-2.73

-0.18

-0.43

-0.35

-0.31

0.78

-0.86

-0.03

-0.28

0.42

c

EA(HF 0 ,1)

d

HF EA(1/ 2)

e

EA(HF 0 ,2 /3)

f

∆MP2-EA

b

2 EA(MP 0)

c

2 EA(MP 0 ,1)

d

MP 2 EA(1/ 2)

e

2 EA(MP 0 ,2 /3)

f

Ref.g

OH

-0.67

-2.70

-0.25

-0.87

-0.70

1.38

1.38

0.43

1.82

1.52

1.83

NH2

-1.49

-3.04

-1.14

-1.66

-1.52

0.26

0.22

-0.45

0.57

0.38

0.74

CH3

-2.03

-2.98

-1.75

-2.16

-2.06

-0.57

-0.77

-1.00

-0.41

-0.48

-0.07

CN

2.79

0.85

2.90

2.73

2.80

4.33

4.99

3.81

4.61

4.26

3.93

O2

-1.12

-2.70

-0.96

-1.20

-1.13

-0.46

0.06

-0.88

-0.26

-0.42

-0.08

MAD_1h

-

1.39

0.24

0.11

0.02

-

0.28

0.47

0.21

0.07

-

MAX_1h

-

2.44

-0.42

0.21

0.03

-

-1.09

1.01

-0.47

-0.14

-

MAD_2i

1.53

2.92

1.29

1.64

1.55

0.45

0.56

0.87

0.30

0.39

-

i

2.56

5.01

2.14

2.77

2.59

0.75

1.14

1.40

-0.68

0.67

-

MAX_2 a

The data set is taken from Refs. 3, 34,35. The unit is eV. The cc-pVQZ basis set is used for all calculations.

b

In the ∆ method, EA is calculated as ∆Emethod = Emethod(N0) – Emethod(N0+1).

c

EA is calculated up to the first order using one-point equation −

d

e

f

g

dEmethod ( λ ) at λ = 0. dncλ

dEmethod ( λ )  1  dEmethod ( λ ) EA is calculated up to the second order using two-point equation −  + .  2  dncλ dncλ λ =0 λ =1  EA is calculated up to the second order using one-point equation −

dEmethod ( λ ) at λ = 1/2. dncλ

 dEmethod ( λ ) 1  dEmethod ( λ ) EA is calculated up to the third order using two-point equation −  + 3⋅ . λ  4  dncλ dn c λ =0 λ = 2/3  The reference data for atoms other than Be and N are taken from Refs. 31,32. Others are taken from Ref. 33 except for CN, which is the CCSD(T) value

calculated in the way as in Ref. 33. h

MAD_1 and MAX_1 refer to the mean absolute deviation and the maximum deviation (Ref. – Calc.), respectively, with respect to the ∆ method results.

i

MAD_2 and MAX_2 refer to the mean absolute deviation and the maximum deviation (Ref. – Calc.), respectively, with respect to the experimental results.



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 32 of 42

Table 3. Comparison with results from the experiments and the ∆ methods: Vertical IPs of core and valence orbitals as calculated using the one-point and two-point equations.a Orb. H2O

CH2O

CH4

NH3

N2

∆HF-IP

b

IP( HF 0)

c

IP( HF 0 ,1)

d

HF IP(1/2 )

e

IP( HF 0 ,2/3)

f

∆MP2-IP b

2 IP( MP 0)

c

2 IP( MP 0,1)

MP 2 IP(1/2 )

e

2 IP( MP 0 ,2/3)

f

Ref.g

1a1

538.91

559.34

538.82

538.95

538.89

539.97

528.25

539.93

540.01

540.16

539.86

2a1

34.03

36.63

34.22

33.94

34.04

33.85

31.42

33.13

34.23

33.77

-

1b1

17.35

19.27

17.52

17.27

17.36

19.02

17.94

18.63

19.21

18.98

18.51

3b1

13.21

15.77

13.47

13.08

13.23

14.93

13.52

14.36

15.21

14.87

14.74

1b2

10.91

13.74

11.21

10.77

10.93

12.69

11.08

12.03

13.01

12.62

12.78

1a1

537.93

559.84

537.91

537.94

537.91

539.68

525.08

539.01

540.08

539.99

539.48

2a1

293.81

308.54

293.61

293.89

293.79

295.48

291.22

296.50

296.12

296.04

294.47

3a1

34.83

38.17

35.16

34.71

34.89

34.35

33.33

34.24

34.88

34.62

-

4a1

22.28

23.72

22.40

22.22

22.29

22.25

20.60

21.64

22.38

22.02

-

1b2

17.70

18.88

17.76

17.67

17.70

17.33

15.90

17.02

17.50

17.32

16.6

5a1

14.54

17.69

14.69

14.46

14.54

16.45

14.45

15.82

16.71

16.40

16.0

1b1

12.31

14.53

12.72

12.22

12.45

14.81

13.74

14.21

14.66

14.43

14.5

2b2

9.41

12.04

9.54

9.34

9.38

11.26

9.43

10.64

11.54

11.25

10.88

1a1

290.53

304.86

290.21

290.70

290.53

290.73

286.51

292.65

290.22

290.88

290.86

2a1

24.23

25.73

24.33

24.19

24.24

23.96

22.84

23.66

24.09

23.92

-

1t2

13.52

14.87

13.60

13.48

13.52

14.46

13.78

14.29

14.54

14.44

14.40

1a1

405.03

422.66

404.81

405.14

405.02

405.81

397.30

411.37

408.27

409.19

405.6

2a1

28.94

31.01

29.09

28.86

28.95

28.67

26.79

28.12

28.93

28.59

-

1e

15.28

17.02

15.40

15.21

15.28

16.64

15.71

16.36

16.78

16.61

16.5

3a1

9.35

11.61

9.59

9.23

9.36

10.87

9.73

10.40

11.09

10.81

10.85

1σg

419.34

426.83

419.46

419.28

419.34

404.91

401.29

404.10

406.56

404.89

409.90

1σu

419.23

426.74

419.35

419.17

419.23

404.82

401.23

404.41

405.91

404.88

409.90

2σg

37.20

39.75

37.33

37.13

37.20

35.71

28.48

32.70

36.39

34.65

-

2σu

20.13

21.27

20.19

20.10

20.14

18.26

16.91

18.00

18.38

18.23

18.78

1πu

15.09

16.51

15.17

15.05

15.09

17.25

16.45

17.11

17.32

17.24

16.98

3σg

15.59

17.18

15.65

15.56

15.59

15.28

13.90

15.09

15.36

15.25

15.60

MAD_1ch

-

14.86

0.16

0.08

0.01

-

7.22

1.49

0.97

0.67

-

h

-

-21.91

0.32

-0.17

0.02

-

14.60

-5.56

-2.46

-3.38

-

h

MAD_1v

-

2.08

0.17

0.07

0.02

-

1.69

0.56

0.24

0.13

-

MAX_1vh

-

-3.34

-0.41

0.14

-0.14

-

7.23

3.01

-0.68

1.06

-

MAD_2ci

3.26

16.96

3.42

3.19

3.27

1.68

8.46

3.06

1.86

2.29

-

i

MAX_2c

-9.44

-20.36

-9.56

-9.38

-9.44

5.08

14.4

5.8

3.99

5.02

-

MAD_2vi

1.36

1.09

1.22

1.42

1.34

0.31

1.12

0.35

0.42

0.28

-

MAX_2vi

2.19

-2.49

1.81

2.28

2.05

-0.73

1.87

0.78

-0.9

-0.72

-

MAX_1c

a

The data set is taken from Ref. 4. The unit is eV. A completely uncontracted cc-pVTZ basis set is used for all calculations.

b

In the ∆ method, IP is calculated as ∆Emethod = Emethod(N0–1) – Emethod(N0).

c

IP is calculated up to the first order using one-point equation −

d

d

dEmethod ( λ ) at λ = 0. dnkλ

dEmethod ( λ )  1  dEmethod ( λ ) IP is calculated up to the second order using two-point equation −  + .  2  dnkλ dnkλ λ =0 λ =1 


Page 33 of 42

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


e

f

IP is calculated up to the second order using one-point equation −

dEmethod ( λ ) at λ = 1/2. dnkλ

 dEmethod ( λ ) 1  dEmethod ( λ ) IP is calculated up to the third order using two-point equation −  + 3⋅ . λ  dnkλ dn 4  k λ =0 λ = 2/3 

g

Experimental ionization potentials are taken from Refs. 4,38 for core and Refs. 4,39 for valence

h

MAD_1 and MAX_1 refer to the mean absolute deviation and the maximum deviation (Ref. – Calc.), respectively, with respect to the ∆ method results, where c

and v refer to core and valence ionization, respectively i

MAD_2 and MAX_2 refer to the mean absolute deviation and the maximum deviation (Ref. – Calc.), respectively, with respect to the experimental results,

where c and v refer to core and valence ionization, respectively



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

Figure 1

Diagrams on the left of the horizontal arrows correspond to the MP2 correlation energy

terms (Eq. 7), while those on the right correspond to the derivatives at the MP2-I level that considers only the explicit dependence of the MP2 correlation energy on orbital occupations (Eqs. 13 and 14). The scissors indicate which line is cut to go from the left to the right.


Page 34 of 42

Page 35 of 42

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


Figure 2

Diagrammatic expressions for the elements of the MP2 correction matrix (a) in the o-o

block (defined in Eq. 35) and (b) in the v-v block (defined in Eq. 36). Here the right bracket indicates that the two lines connected, i.e., (i,j) or (a,b), have the same indices, while the right bracket with a cross

indicates that the line indices for the two lines connected are different.



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

Figure 3

Page 36 of 42

Diagrammatic expressions corresponding to contributions from (a) the o-o and (b) the v-v

( 2 ), λ ( 2 ), λ blocks of the MP2 correction matrix, Po-o and Pv-v , to IP (when t represents a hole line index),

or EA (when t represents a particle line index). Diagrams on the left of the horizontal arrows correspond to third-order energy terms, where the scissors indicate which line is cut to obtain the diagrams on the right of the horizontal arrows.

36


Page 37 of 42

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


Figure 4

Diagrammatic expressions corresponding to contributions from the v-o and o-v blocks of

( ) ( ) and Po-v,(1) , to IP (when t represents a hole the MP2 correction matrix in the first iteration, Pv-o,(1) 2 ,λ

2 ,λ

line index), or EA (when t represents a particle line index), where diagrams (a)-(f) correspond to terms (a)-(f) of Eq. 43.

37



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

Figure 5

Page 38 of 42

Diagrammatic expressions corresponding to contributions from the v-o and o-v blocks of

( ) ( ) the MP2 correction matrix in the n-th iteration, Pv-o,(n) and Po-v,(n) , to IP (when t represents a hole 2 ,λ

2 ,λ

line index), or EA (when t represents a particle line index), where diagrams (a)-(f) are generated from (a)-(f) of FIG. 4 by adding

and

order by order.

38


Page 39 of 42

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


Figure 6. The dressed interaction expressed by a wavy line.

39



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 40 of 42

Figure 7. Diagrammatic expressions corresponding to contributions from Parts II and III to IP (when t represents a hole line index), or EA (when t represents a particle line index), where diagrams (a)-(f) are obtained by summing up all terms order by order as shown in Figure 5, while diagrams (e) and (f) include also the o-o and v-v contributions shown in Figure 3.

40


Page 41 of 42

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


Figure 8.

Deviations (Ref. – Calc.) of energy derivatives as a function of the occupation number of

partially occupied orbitals: (a) 1s, (b) 2s, (c) 2p, and (d) 3s of Ne atom. MP2 numerical energy derivatives are taken as the references, while MP2-I, II, and III refer to the three levels of approximation defined in Eq. 11. The energy unit is in eV, where data in (d) has been amplified by 103. Data points were recorded in intervals of 0.1 electron.

41



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 42 of 42

Figure 9. Behaviors of (a) relative MP2 total energies with respect to the neutral system, (b) deviations from the corresponding linearity interpolations, and (c) full MP2 energy derivatives as functions of occupation number of the innermost orbital in N2. The energy unit is in eV. Data points were recorded in intervals of 0.1 electron.

42


Second-Order Perturbation Theory for Fractional Occupation

Recommend Documents