Applications to Ionization Potential and Electron Affinity Calcula

Subscriber access provided by TEXAS A&M INTL UNIV

Article

The integration approach at the second-order perturbation theory: 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.5b00591 • Publication Date (Web): 09 Sep 2015 Downloaded from http://pubs.acs.org on September 10, 2015

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

The integration approach at the second-order perturbation theory: 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

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

Page 2 of 42

An integration approach is developed to calculate ionization potentials (IPs), and electron affinities (EAs), which is an extension of the D-∆MBPT(2) method [A. Beste, A. Vazquez-Mayagoitia, and J. V. Ortiz, J. Chem. Phys. 138, 074101(2013)]. The latter is an extension of the single-point method of Cohen et al. [A. J. Cohen, P. Mori-Sánchez, and W. Yang, J. Chem. Theory Comput. 5, 786 (2009)] from the perspective of fractional charges. While relaxation effects were included only at the Hartree-Fock (HF) level in the previous methods, such effects are fully taken into account in the present method up to the second-order Møller-Plesset (MP2) level. This is made possible by deriving the full MP2 energy gradient with respect to the orbital occupation numbers, which is solved through the coupled-perturbed HF (CP-HF) equations.

*Corresponding author: [email protected]

2


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 Ionization potential (IP) and electron affinity (EA) are essential molecular properties of fundamental importance.1-4 Many useful concepts, such as electronegativity,5 chemical potential,6 hardness and softness,7 as well as electrophilicity and nucleophilicity,8 etc., are defined on top of IP and EA. Vertical IP and EA, for the fixed atomic positions at the corresponding N0-electron system, are defined as:

IP = E ( N 0 − 1) − E ( N 0 )

(1)

EA = E ( N 0 ) − E ( N 0 + 1)

(2)

where E ( N 0 ) is the ground state energy of the N0-electron system with N0 being an integer, and

E ( N 0 m 1) correspond to the ground state energies of the

( N0 m 1) -electron

systems after an

electron is removed from or added to the frontier orbitals, i.e., the highest occupied molecular orbital (HOMO) or the lowest unoccupied molecular orbital (LUMO), of the N0-electron system. These equations involve two separate calculations, which can be carried out either at the Hartree-Fock (HF) level or at the post-HF level through e.g. Møller-Plesset (MP) perturbation theory, etc. and are generally referred to as the ∆ methods such as ∆HF9,10 or ∆MP11,12, respectively. One of the disadvantages of the ∆ methods is that the energy difference (e.g. IP or EA) is typically much smaller than the total energies of the initial and final states, i.e. a small number computed from the subtraction of two large numbers. Hence, there is potentially a precision problem, in particular, when the system size is large. The precision problem can be avoided in the direct (D) methods (see, e.g. 13-28), where IPs and EAs are computed directly without explicit involving total energy differences of the initial and final states. The simplest and once frequently used D method is the HF method under Koopmans’ 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

Page 4 of 42

theorem,13 according to which IPs and EAs are given by the negatives of the respective orbital energies for the occupied or unoccupied MOs, respectively, in the N0-electron system.29 This approximation corresponds to the frozen orbital approximation, where the effects of electron relaxation during the electron removal and addition are neglected. At the HF level, the effects of electron correlation are also neglected. There are a number of well-known cases, e.g. N2 and F2,30,31 where Koopmans’ theorem fails to predict the correct ordering of states, and the inclusion of correlation effects were found crucial.32-33 Over decades, many methods have been developed to compute directly IPs and EAs at the correlated level.14-28 For example, there are the equation-of-motion coupled cluster (EOM-CC) method,18-20 the coupled cluster Green’s function (CCGF) method,21 the coupled cluster linear response theory (CCLRT),22 the Fock space coupled cluster (FSCC) method,23,24 the SAC (symmetry-adapted cluster)/SAC-CI (configuration interaction) method,25 the electron propagator method, etc.26-28 These CC-based methods are closely related to each other and are generally computationally nontrivial. While higher level methods tend to be expensive, the HF method is known to be deficient without considering correlation effects. Thus MP2 is often applied as the lowest correlation level in the family of wave function methods.34-46 Recently, Cohen, Mori-Sánchez, and Yang developed an analytical derivative of MP2 by considering the explicit dependence of correlation energies on the orbital occupation numbers,45 leading to an expression that may be related to the MP2 single-particle energy. Their expression can also be obtained from the second-order self-energy in propagator theory,28 which corresponds to the lowest level Green’s function approach.29 More recently, Beste, Vázquez-Mayagoitia, and Ortiz proposed the D-∆MBPT(2) method for IP and EA calculations.46 4


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


They extended the work of Cohen et al,45 by further considering the explicit dependence of orbital energies on the occupation number,47 yielding some higher-order terms that are not present in Cohen’s expression. However, the expression of Beste et al. is still incomplete,46 as the orbital relaxation effects beyond HF are neglected. In analogous to the D-∆HF method,48 they employed a numeric integration technique of the energy derivative in the D-∆MBPT(2) method. Here, we will further extend the work of Beste et al.,46 by considering the full MP2 energy gradient with respect to the orbital occupation numbers. This is realized by solving the coupled-perturbed HF (CP-HF) equations. The reminder of this paper is organized as follows. In Sec. II, we briefly discuss the integration approach to obtain energy differences such as IPs, and EAs, where one-point and two-point equations are employed. We then derive the full MP2 energy derivatives with respect to the orbital occupation numbers as compared to the formalisms of Cohen et al.45 and Beste et al.46 Computation details are summarized in Sec. III, while test calculations are presented in Sec. IV for some selective atoms and molecules. Finally, we give some concluding remarks in Sec. V.

II. THEORY A. The Integration Approach The integration approach, its significance, as well as its connection to other methods, have been discussed in depth by Beste et al.46,48 Here we only recapitulate it briefly. The integration approach introduces a coupling parameter λ that defines a pathway to connect two different states of the system. The energy difference between the final state (λ=1) and the initial state (λ=0) is then calculated as an integral of the derivative of energy with respect to parameter λ: 5



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

1

dE ( λ )

0

dλ

∆E = E (1) − E ( 0 ) = ∫

As the energy derivative,

dE ( λ )

dλ

Page 6 of 42

(3)

, is on the same order of magnitude as the energy difference, ∆E ,

dλ

itself, the precision problem is, therefore, avoided. In terms of standard perturbation theory,29 one can always expand the energy as a polynomial function of λ: E ( λ ) = ∑ E ( k )λ k = E ( 0 ) + E (1) λ + E ( 2 )λ 2 + E (3) λ 3 + L

(4)

such that the derivative is given by dE ( λ ) dλ

= E (1) + 2 E (2) λ + 3E (3) λ 2 +L

(5)

If E ( λ ) changes linearly with parameter λ , ∆E is then defined by the initial slope of E ( λ ) , i.e.

E (1) =

dE ( λ ) . Normally, the E ( λ ) curve is not strictly linear, the initial slope only provides an d λ λ =0

approximation up to the first order

∫

1

0

dE ( λ ) dE ( λ ) d λ ≈ ∆ E [ →1] = dλ d λ λ =0

(6)

Usually, the quadratic term is the leading term that accounts for the nonlinearity of E ( λ ) . Hence we arrive at a two-point equation for the integral that is exact up to the second order in terms of coupling parameter λ 1

dE ( λ )

0

dλ

∫

d λ ≈ ∆E [ → 2] =

dE ( λ )  1  dE ( λ ) +   2  d λ λ = 0 d λ λ =1 

(7)

Generally, the cubic term makes relatively small contribution. Note that ∆ E [ → 2 ] leads to an error of 1 (3) E + L , where contributions from the higher order terms are exaggerated. 2

For an accurate and efficient calculation of the integral, the energy derivative with respect to 6


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


parameter λ shall be readily computed accurately and efficiently. For IP or EA calculations that we are interested in here, an obvious choice to define the integration pathway is to relate parameter λ to the orbital occupation number, ntλ , where the electron detachment or attachment occurs at the t-th orbital. Here and in the following, we use the labels p, q, r, s, and t to stand for any canonical HF spin-orbitals. While the orbital occupation number of p, q, r, and s can be either 0 or 1, but not in between, the orbital occupation number of the t-th orbital fulfills ntλ ∈ [ 0,1] . We also adopt the common labels where i, j, and k stand for the occupied canonical HF spin-orbitals, and a, b, and c for the virtual canonical HF spin-orbitals. In particular, t = k refers to an IP calculation, while t = c refers to an EA calculation, as will be discussed below. Let the N 0 -electron system correspond to the initial state, and the

( N0 m 1) -electron systems

correspond to the final state. For IP calculations from Eq. (3), we assume that the electron in the k-th orbital is ionized, and a series of hypothetical fractional charge systems are generated to connect the initial ( nk0 = 1 ) and final ( n1k = 0 ) states. Along this path, only the k-th orbital can be fractionally occupied, while the remaining orbitals keep integer occupations, such that we have

λ = 1 − nkλ i ≤ N 0 & i ≠ k : niλ = 1

(8)

λ

a > N 0 : na = 0

Analogously, a series of hypothetical fractional charge systems are generated in order to carry out an EA calculation, which connect the initial ( nc0 = 0 ) and final ( n1c = 1 ) states where only the c-th orbital can be fractionally occupied:

λ = ncλ a > N 0 & a ≠ c : naλ = 0

(9)

λ

i ≤ N 0 : ni = 1

Hence 7



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

1

dE ( λ )

0

dλ

IP = ∫

1

dE ( λ )

0

dnkλ

d λ = −∫

Page 8 of 42

dλ

(10)

and 1

dE ( λ )

0

dλ

EA = − ∫

1

dE ( λ )

0

dncλ

d λ = −∫

dλ

(11)

We are now in a position to explore how energies of the systems depend on the occupation numbers.

B. The Hartree-Fock Method with Fractional Charges The HF total energy may be formulated as all

all

1 EHF ( λ ) = ∑ Ppqλ p λ hˆ q λ + ∑ Ppqλ Prsλ p λ r λ q λ s λ 2 pqrs pq

(12)

1 v where hˆ = − ∇ˆ 2 + υˆext ( r ) is the one-electron operator which contains the kinetic operator and an 2

external

potential

υˆext

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

usually

coming

from

the

nuclei,

while

pλ r λ qλ sλ

is an antisymmetrized combination of the regular two-electron

repulsion integrals: p λ r λ q λ s λ = ∫∫ φ pλ ( x1 )φqλ ( x1 ) r12−1φrλ ( x 2 ) φ sλ ( x 2 ) d x1d x 2 .

(13)

where x1 and x2 are combined spatial and spin coordinates. In these equations, p λ denotes the p-th canonical HF orbital ( φ pλ ) from the self-consistent-field (SCF) calculation which depends on λ at a chosen set of occupation number, and Ppqλ is an element of the one-particle density matrix P, where only the diagonal element can have a non-zero value that corresponds to orbital occupation n λp : Ppqλ = δ pq n λp

(14)

In the usual HF scheme, the occupation numbers assume values 1 or 0 for occupied and unoccupied orbitals, respectively. Here they are extended according to Eq. (8) or (9) to allow fractional occupations. Minimization of the HF total energy, Eq. (12), is achieved by varying one of the 8


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


spin-orbitals while preserving its normalization and orthogonality to the other spin-orbitals. This gives the HF equations, which yield the orbital energies as all

ε pλ = p λ hˆ p λ + ∑ Prsλ p λ r λ p λ s λ

(15)

rs ≠ p

On the other hand, one can differentiate the HF total energy, Eq. (12), with respect to the occupation number ntλ , leading to: all dEHF ( λ ) λ ˆ λ = t h t + ∑ Ppqλ p λ t λ q λ t λ dntλ pq ≠ t

(16)

In deriving Eq. (16), we have recognized the fact that the HF energy is stationary with respect to the orbital changes, as we use the HF canonical orbitals {φ pλ } . Comparison of Eq. (16) with Eq. (15) leads to dEHF ( λ ) = ε tλ dntλ

(17)

Hence, IPs and EAs can readily be calculated at the HF level as IPHF = ∫ ( −ε kλ ) d λ

(18)

EAHF = ∫ ( −ε cλ ) d λ

(19)

1

0

1

0

Approximations up to the first order and the second order for IPs are [ →1] IPHF (0) = −

dEHF ( λ ) dnkλ

= −ε k0 λ =0

(20)

occ

= − k 0 hˆ k 0 − ∑ k 0 i 0 k 0i 0 i≠k

and [ → 2] IPHF =−

dE HF ( λ )  1  dE HF ( λ ) 1 0 1 +   = − (ε k + ε k ) λ λ 2  dnk dn 2 k  λ =0 λ =1 

(21)

In Eq. (21), the SCF calculation needs to be performed also at the final state for λ = 1 , leading to [ →1] IPHF (1) = −

dE HF ( λ ) dnkλ

, where the original k-th occupied orbital for λ = 0 has just become λ =1

9



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

unoccupied. EAs can be calculated similarly as [ →1] EAHF (0) = −

dE HF ( λ ) dncλ

= − c hˆ c 0

= −ε c0

(22)

λ =0 occ

0

−∑ c i

0 0

0 0

ci

i

[ → 2] EAHF =−

dE HF ( λ )  1  dE HF ( λ ) 1 0 1 +   = − (ε c + ε c ) λ λ dnc 2  dnc 2  λ =0 λ =1 

(23)

while, in Eq. (23), the label c pertains to the neutral system for λ = 0 , which is originally [ →1] 1 unoccupied, has just become occupied for λ = 1 when calculating EAHF (1) = −ε c .

Note that Eq. (17) is related to Janak theorem,49 where HF may be considered as an approximate density functional, Eqs. (20) and (22) are related to Koopmans’ theorem,13 where orbitals are frozen at the N 0 -electron system for λ = 0 . The relaxation effects can be introduced via two-point equations of Eqs. (21) and (23).

C. The MP2 Method with Fractional Charges The total energy of MP2 is given by EMP 2 ( λ ) = E HF ( λ ) + Ec , MP 2 ( λ )

(24)

where the last expression for the MP2 correlation energy, generalized to include occupation numbers,45,46,50,51 is pλ qλ r λ sλ

all

2

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

(25)

We now need to find the derivative of the MP2 correlation energy Ec , MP 2 ( λ ) with respect to the occupation number ntλ : dEc , MP 2 ( λ ) λ

dnt

=

∂Ec , MP 2 ( λ ) λ

∂nt

all

+∑ p

∂Ec , MP 2 ( λ ) d ε pλ λ

∂ε p

λ

dnt

all

+∑ p

∂Ec , MP 2 ( λ ) dφ pλ ∂φ pλ

dntλ

(26) 10


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


Eq. (15) shows how ε λp depends explicitly on the occupation numbers of all orbitals other than n λp : ∂ε λp = pλ t λ pλt λ ∂ntλ

(27)

The full derivative of orbital energy with respect to occupation number is given by d ε pλ ∂ε pλ all ∂ε pλ d φqλ = +∑ λ λ dntλ ∂ntλ q ∂φq dnt

(28)

Hence we can group the full MP2 correlation energy derivative into three terms: dEc , MP 2 ( λ ) ∂Ec , MP 2 ( λ ) = dntλ ∂ntλ all ∂Ec , MP 2 ( λ ) ∂ε λp +∑ ∂ε pλ ∂ntλ p

(I) (II)

all  all ∂Ec , MP 2 ( λ ) ∂ε λp ∂Ec , MP 2 ( λ )  dφqλ + ∑∑ +  λ ∂φqλ ∂ε pλ ∂φqλ q  p  dnt

(29)

(III)

Term I accounts for the explicit dependence of Ec , MP 2 ( λ ) on the occupation numbers, leading to45,46 ∂Ec , MP 2 ( λ ) ∂ntλ

λ λ λ λ λ λ λ 1 all n p (1 − nr )(1 − ns ) p t r s = ∑ 2 prs ε pλ + ε tλ − ε rλ − ε sλ

−

all

1 ∑ 2 pqs

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

2

(30)

2

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

If contributions from Terms II and III are omitted, and the energy derivative is only calculated at the initial state for λ = 0 , this is practically the prescription given by Cohen, et al. as the MP2 single-particle energy:45 I ,[ →1] IPMP 2(0) = −

I dEMP 2 (λ ) dnkλ λ =0

2 2 0 0 0 0   occ vir i 0 j 0 k 0b0 1  occ vir k j a b 0  = −ε k − ∑ ∑ 0 −∑∑ 2  j ab ε k + ε 0j − ε a0 − ε b0 ij b ε i0 + ε 0j − ε k0 − ε b0   

(31)

11



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 ,[ →1] EAMP 2(0) = −

Page 12 of 42

I dE MP 2 (λ ) dncλ λ =0

2 2 0 0 0 0   occ vir i 0 j 0 c 0b 0 1  occ vir c j a b 0  = −ε c − ∑ ∑ 0 −∑∑ 2  j ab ε c + ε 0j − ε a0 − ε b0 ij b ε i0 + ε 0j − ε c0 − ε b0   

(32)

Note that we have put a superscript, I, in Eqs. (31) and (32) to show that only Term I has been considered in computing the MP2 correlation energy derivatives. These expressions can also be obtained from the second-order self-energy in propagator theory,28 which corresponds to the lowest level Green’s function approach,29 and has also been used to calculate the quasiparticle band gap of solids.52,53 Under the Term I approximation, IPs and EAs can also be evaluated up to the second order by using the two-point equation: I I  dEMP 1  dEMP 2 ( λ ) 2 (λ ) I ,[ → 2] IPMP =−  +  2  2  dnkλ dnkλ λ =0 λ =1 

(33)

I I  dEMP 1  dEMP 2 ( λ ) 2 (λ ) I ,[ → 2] =−  + EAMP  2  2  dncλ dncλ λ =0 λ =1 

(34)

where the energy derivatives need to be calculated also at the final state for λ = 1 . Beste et al. extended the work of Cohen et al,45 by further considering the explicit dependence of orbital energies on the occupation number (i.e. Term II), yielding some higher-order terms (HOT47) that are not present in Cohen’s expression: all

∑ p

∂ε pλ

2

pλ qλ r λ sλ 1 all λ λ λ λ = − n n 1 − n 1 − n ∑ p q ( r )( s ) λ λ λ λ 2 ∂ntλ 4 pqrs (ε p + ε q − ε r − ε s )

∂Ec , MP 2 ( λ ) ∂ε λp

(35)

⋅  t λ p λ t λ p λ + t λ q λ t λ q λ − t λ r λ t λ r λ − t λ s λ t λ s λ 

If the t-th orbital happens to be an occupied orbital, we have, at λ = 0 ,

12


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


∂Ec , MP 2 ( λ ) ∂ε pλ ∑p ∂ε pλ ∂nkλ

2

i 0 j 0 a 0b 0 1 occ vir = − ∑∑ 4 ij ab ( ε i0 + ε 0j − ε a0 − ε b0 )2

all

λ =0

(36)

⋅  k 0i 0 k 0i 0 + k 0 j 0 k 0 j 0 − k 0 a 0 k 0 a 0 − k 0b 0 k 0b 0 

If the t-th orbital is an unoccupied orbital, we have, at λ = 0 , all

∑ p

2

∂Ec , MP 2 ( λ ) ∂ε λp ∂ε pλ

∂ncλ

i 0 j 0 a 0b 0 1 occ vir = − ∑∑ 4 ij ab ( ε i0 + ε 0j − ε a0 − ε b0 )2

λ =0

(37)

⋅  c 0i 0 c 0i 0 + c 0 j 0 c 0 j 0 − c 0 a 0 c 0 a 0 − c 0b 0 c 0b 0 

Similar expressions can be obtained at λ = 1 . Adding Eqs. (30) and (35) together gives an approximate MP2 correlation energy derivative as developed by Beste et al.47 Here we label it as the Term II approximation with a superscript, II: dEcII, MP 2 ( λ ) dntλ

=

∂Ec , MP 2 ( λ ) ∂ntλ

all

+∑ p

∂Ec , MP 2 ( λ ) ∂ε λp ∂ε pλ

∂ntλ

(38)

II ,[ →2] II ,[ → 2] II ,[ →1] II ,[ →1] , and EAMP can be calculated, up to the first order at λ = 0 Hence, IPMP 2 2( λ ) , EAMP 2( λ ) , IPMP 2

or 1, and the second order, respectively, according to Eqs. (31) ‒ (34) after updating

I dEMP 2 (λ ) to dntλ

II dE MP 2 (λ ) . dntλ

We now proceed to develop the full MP2 derivative to take into account the contributions from all three terms shown in Eq. (29) or (27).

D. The Full MP2 Derivative with Respect to Orbital Occupations The full MP2 correlation energy derivative with respect to orbital occupations is realized here by solving the coupled-perturbed (CP) HF equations, where the formalism can be related to the well-established algorithm for MP2 geometry optimization or calculations of other response properties.54-58 Thus Term I can be analogously referred to as the direct part, depending explicitly on 13



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

orbital occupations, while Terms II and III as the indirect part coming from the orbitals’ responses, which depend inexplicitly on orbital occupations. We assume that the orbital derivatives with respect to the t-th orbital occupation can be expanded by the canonical HF orbitals {φqλ } : d φ pλ all λ λ ,nt = ∑ φq U pq dntλ q

(39)

λ ,n where U pq is an element in U λ ,nt , which is a rotation matrix that allows for orbital mixing due to t

the changing of orbital occupations. As we are using one-point and two-point equations to approximate the integration approach, only the derivatives at λ = 0 or λ = 1 are desired. We will then drop the λ-dependence for partial charges and only consider the cases for integer occupations for the time being. t The virtual-occupied (v-o) block U nv-o have to be determined by solving the CP-HF

equations:54-58 A′U nv-ot = B v-o,t

(40)

For a given spin σ (α or β), the matrix elements are defined by A′pσσ qσ , rs = − ( ε r − ε s ) δ rpσ δ sqσ − Aσpσ qσ , rs

(41)

Aσpσ qσ , rs = pσ r qσ s + pσ s qσ r

(42)

B σpσ qσ ,t = pσ t qσ t

(43)

and

t To settle the full U nt matrix, we need additional blocks, i.e. U no-o from occupied-occupied (o-o) and

t U nv-v virtual-virtual (v-v) orbitals: 54-58

U bcnt =

1   Abc ,aiU aint + Bbc ,t  ∑  ε c − ε b  ia 

(44)

14


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


U kjnt =

1   Akj ,aiU aint + Bkj ,t  ∑  ε j − ε k  ia 

(45)

The U nt matrix has the property nt U pq + U qpnt = 0

(46)

In analogous, the full derivatives of orbital energies can be expressed as dε p = ∑ App , aiU aint + B pp ,t dnt ia

(47)

while B pp ,t , as suggested by Eq. (47), corresponds to Eq. (27) for

∂ε p . ∂nt

The contributions of the indirect part from Terms II and III can be formulated as  ∂E c , MP 2 ( λ ) d ε pλ ∂E c , MP 2 ( λ ) d φ pλ ∑p  ∂ε pλ dntλ + ∂φ pλ dntλ  all

Here Pij(

2)

 (2) (2) MP 2 nt  = ∑ Pij Bij ,t + ∑ Pab Bab ,t + 2 ∑ Lai U ai ij ab ia 

(48)

and Pab( ) are the elements of the standard MP2 correction matrix P(2) to the one-particle 2

density matrix, which correspond to the o-o block and the v-o block, respectively. They are defined as54-58 Pij( ) = − 1 ∑ tikab t ab jk 2 kab

(49)

Pab( 2 ) = 1 ∑ tijac tijbc 2 ijc

(50)

2

and

where tijab is the double substitution amplitude54-58 tijab = ij ab

(ε i + ε j − ε a − ε b )

(51)

MP 2 59-61

2 refers to the element of the well-known MP2 Lagrangian L LMP ai

,

2 LMP = − 1 ∑ t ab jk ib + 1 ∑ tijcb aj cb + ∑ Pjk( 2) ja ki + ∑ Pbc( 2 ) ba ci ai jk 2 jkb 2 jcb jk bc

(52)

which is related to the v-o block of P(2), well-known as the Z-vector method59 ( 2) Pv-o A′ = LMP 2

(53)

From Eq. (40), one has 15



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

LMP 2 ⋅ Unv-ot = LMP2 A′−1Bv-o,t .

Page 16 of 42

(54)

( )

( ) ( 2) = transpose Pv-o , and Eqs. (53) and (54), the last term of Eq. (48) can be reformulated Noting Po-v 2

as ( 2) ( 2) 2LMP 2 ⋅ Unv-ot = Pv-o ⋅ Bv-o,t + Po-v ⋅ Bo-v,t

(55)

Finally, terms on the RHS of Eq. (48) can be joined together to have  ∂E c , MP 2 ( λ ) d ε pλ ∂E c , MP 2 ( λ ) d φ pλ ∑p  ∂ε pλ dntλ + ∂φ pλ dntλ  all

all  all ( 2 ) (2)  = ∑ Ppq B pq ,t = ∑ Ppq tp tq pq pq 

(56)

Thus, instead of directly solving the CP-HF equations for U nt , we use the Z-vector method to find the MP2 correction matrix P(2), which fixes the contributions from Terms II and III for λ = 0 and λ = 1 . Combining Eqs. (56) and (30), we arrive at the full MP2 correlation energy derivative with

respect to the orbital occupation number

dEc , MP 2 ( λ ) , which can be used to calculate IPs and EAs dntλ

III ,[ → 2] III ,[ → 2] III ,[ →1] III ,[ →1] with one-point (i.e., IPMP , and EAMP ) 2 2( λ ) , EAMP 2 ( λ ) at λ = 0 or λ = 1) and two-point (i.e. IPMP 2

equations. The superscript, III, remains, which is just a reminder of the fact that the MP2 derivative with respect to the orbital occupation number has been developed successively to include the explicit dependence of the orbital occupation numbers (I),45 orbital eigenvalues (II)46 and now orbitals themselves (III).

III. COMPUTATIONAL DETAILS The afore-described formalisms at different levels of approximations have been implemented in the NorthWest computational Chemistry (NWChem) software package, version 6.1.1.62 All calculations were performed in the spin-unrestricted formalisms. No spherical symmetry requirements were imposed on electron densities in atoms. In the numerical calculations of the 16


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


derivatives

∆E , ∆nt was chosen as 0.001.45 For benchmarking purpose, the variable, tol2e in ∆nt

NWChem,62 which determined the integral screening threshold for the evaluation of the energy and related Fock-like matrices was set very tightly to be 1E-16. The SCF convergence was satisfied, if the total energy difference of two successive cycles was less than 1E-11 au. We discussed here only the frontier orbitals, i.e., HOMO or LUMO, although our formalisms described above are general that can be also applied to core electrons. All orbitals during SCF were occupied according to the aufbau principle, where only HOMO/LUMO were allowed to be fractionally occupied. Note that the orbital label k or c pertains to the N0 system for λ = 0, where k = HOMO and c = LUMO. For the final state of λ = 1, the k-th orbital becomes LUMO of the (N0-1) system when its electron in is ionized, while the c-th orbital becomes HOMO of the (N0+1) system when an electron is added to this orbital. Numerical testing is performed for 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). This data set was compiled by Cohen, Mori-Sánchez and Yang in discussing their fractional charge perspective on the band gap problem in density functional theory, 63 as well as MP2.45 The same set has been used recently to examine the fractional charge behavior of some doubly hybrid functionals.64-66 The experimental IP and EA data were generally taken from Refs. 67, and 68, which were used as reference data for statistics analyses. There are no experimental data available for EAs of Be, N and F2. Hence the CCSD(T) data from Ref. 69 were used as references instead. As we are interested in VIPs, geometries for cations were fixed to be those in the respective neutral molecules,69-71 i.e., F2 (1.4119 Å), OH (0.9697 Å), NH2 (1.02402 Å, 103.40°), CH3 (1.079 Å, 120.0°), CN (1.1718 Å), and O2 (1.2075 Å). For EA calculations, the same experimental geometries of the neutral molecules were also adopted for the 17



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

anions. The basis set used in all calculations reported in the present work is cc-pVQZ.72,73

IV. RESULTS AND DISSCUSION A. The MP2 Correlation Energy Derivatives at Different Levels of Approximations The calculated MP2 correlation energy derivatives with respect to the orbital occupation numbers ntλ are collected in Tables 1 and 2. While data in Table 1 correspond to the case t = k where electron

is to be ionized from the k-th occupied orbital, those in Table 2 correspond to the case t = c where electron is to be added to the c-th unoccupied orbital. The derivatives,

dEcL, MP 2 ( λ ) , are evaluated at dntλ

both λ = 0 and 1 using different levels of approximations (L = I, II, and III). Here L = I refers to the approximation where only explicit dependence of the MP2 correlation energy on the occupation number is considered as in Ref. 45. L = II refers to additional consideration of the explicit dependence of the HF eigenvalues on the occupation number as in Ref. 46. L = III is the full MP2 correlation energy derivative developed in the present work. Tables 1 and 2 also summarize the numerical results from finite differences, labeled as L = Num. These numbers serve as references here, although they themselves may suffer from numerical unsteadiness. From Table 1, it is seen that the MP2 correlation derivatives for electron ionizations are all positive at λ = 0. They become all negative at λ = 1. Table 2 shows the opposite trend, where the MP2 correlation derivatives for electron attachment are all negative at λ = 0, and they become all positive at λ = 1. These observations show that the correlation effects to the HF orbital energies depend on at which point (e.g. initial state for λ = 0 or final state for λ = 1) the evaluations are carried out and which orbitals (i.e., occupied or unoccupied) are concerned. 18


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


For atomic systems, correlation effects generally become larger along with the increasing electronegatives. There are shell structures when fully/half occupied/unoccupied orbitals are involved. The negative value is generally higher than the positive one in magnitude for a given N0-species. For molecular systems, the correlation contributions also become more significant if more electronegative atoms present in the molecule (e.g. CH3 < NH2 < OH). As compared to L = Num, L = I presents a reasonably good approximation with mean absolute deviations (MADs) in the range of 0.16 – 0.30 eV. The majority of the data from L = I is smaller in magnitude than those from L = Num, regardless how the derivatives are evaluated. The maximum deviation (MAX, Ref. – Calc.) amounts to –1.34 eV occurring at CN when calculating the initial derivative of the electron attachment. Data in Tables 1 and 2 suggest that L = II is not necessary a better approximation than L = I. MADs for L = II range from 0.22 to 0.41 eV, which almost consistently underestimates the correlation contribution as compared to L = I when the derivatives are positive for an occupied orbital (i.e., nk0 or n1c ); whereas overestimates the correlation contribution as compared to L = I when the derivatives are negative for an unoccupied orbital (i.e., n1k or nc0 ).  ∆EcNum dEcL, MP 2 ( λ )  , MP 2 ( λ ) Figure 1 depicts the derivative differences,  −  , at different levels of λ λ ∆ n dn t t  

approximations using those from L = Num as references. While Figure 1 (a) and (b) show the cases of t = k for IP calculations using the initial (λ = 0) and final (λ = 1) derivatives, respectively, Figure 1 (c) and (d) show the cases of t = c for EA calculations using the initial (λ = 0) and final (λ = 1) derivatives, respectively. As is clear, L = II displays, in many cases, a higher error bar than does L = I. Ideally, one should find that L = III delivers identical results to those from L = Num. This is indeed the case. As compared to the data from L = Num, L= III consistently leads to MAD = 0.00 eV, 19



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

Page 20 of 42

demonstrating that the analytical gradient formalism developed in the present work faithfully takes into account the full MP2 correlation effects upon changing the orbital occupation numbers.

B. Fractional Charge Behavior for Total Energy vs. Electron Number N From the perspective of fractional charges, it is now well-known that the total energy of a system with a fractional number of electrons N = N 0 ± δ with 0 ≤ δ ≤ 1 is a straight line connecting the total energies at the adjacent integers. This was first proved in density functional theory for the grand canonical ensemble,74 and then confirmed for the pure state in consideration of finite systems in the dissociation limit.75 Recently, the extension of wave function theories to fractional charges have been explored by Steinmann and Yang.76 In recognition of the concave behavior of the HF total energy with respect to the fractional occupations, the straight line behavior is increasingly improved with increasingly sophistical treatment of correlation effects from MP2 to coupled-cluster singles and doubles (CCSD) to additional perturbative triples, CCSD(T). For any system with an approximate method, one can check the fractional charge behaviors against the linearity condition.45,63,64,76,77 Figure 2 shows, using the O atom as an example, the ground state energy E as a function of the electron number N. The exact straight lines are obtained using the experimental IP (13.61 eV65,66) and EA (1.46 eV65,66) of the O atom. As is clear from Figure 2, MP2 exhibits a much improved linearity as compared to HF, signifying the great importance of the correlation effects. The inset in Figure 2 shows the deviation from linearity for a given method in another way. It depicts the energy deviations from its own straight line interpolation between adjacent integer values for the O atom (i.e. O+ O O‒). Here the positive deviations suggest a concave error associated 20


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


with the HF method, where the largest deviations occur at the respective middle points with the k-th or c-th orbital being half occupied. The errors associated with MP2 for the O atom are obviously smaller than those of HF. While the curves from the HF method are almost quadratic, the MP2 errors change sign at the middle points, having minimum errors at the middle points.

C. Ionization Potential and Electron Affinity Calculated with the HF method Table 3 summarizes the IPs calculated with the HF method using various approximations. Data in the [ →1] are just the opposite of the orbital energies for the HOMOs of the column under the title of IPHF (0 ) [ →1] 0 N0-systems, i.e., IPHF . This corresponds to the well-known Koopmans’ theorem.13 As ( 0 ) = − ε HOMO [ →1] orbital relaxation effects are not taken into account, IPHF yields numbers which are consistently (0 )

too high as compared to the data from ∆HF. On the other hand, data in the column under the title of [ →1] are just the opposite of the orbital energies for the LUMO of the (N0‒1)-systems, i.e., IPHF (1) [ →1] 1 . If orbitals are fixed at the (N0‒1)-system, which underestimates the stability of the IPHF (1) = − ε LUMO

[ →1] N0-system, IPHF yields IPs which are consistently too low as compared to the data from ∆HF. The (1) [→2] [ →1] [ →1] combines the two opposite effects of IPHF and IPHF , showing two-point equation of IPHF (0 ) (1) [→2] good agreement with ∆HF. As shown in Table 3, MAD_1/MAX_1 from IPHF are 0.15/0.31 eV, [ →1] [ →1] and IPHF . which are one order of magnitude smaller than those of IPHF (0 ) (1)

Traditionally, Koopmans’s theorem is often used to predict VIPs. This benefits from the lack of both relaxation effect and correlation effect which tend to cancel each other.11,33 Thus, as shown in Table 3, ∆HF consistently underestimates IPs as compared to the experimental values, leading to MAD_2 = 0.91 eV

[ →1] with MAX_2 = 1.77 eV, whereas IPHF consistently overestimates IPs as (0 )

compared to the experimental values, leading to MAD_2 = 1.08 eV with MAX_2 = ‒2.88 eV. Note 21



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

Page 22 of 42

[ →1] that values from IPHF are generally too small, whose MAD_2 and MAX_2 are 1.97 and 4.01 eV, (1)

[→2] respectively. IPHF seems to provide the smallest MAD and MAX of 0.80 and 1.48 eV,

respectively. Figure 2 provides the other way, from the perspective of fractional charges, to look into the HF errors for IP predictions. From N0 = 8 to the left, electron flows out of the O atom, corresponds to the [ →1] ionization process. IPHF yields the initial derivative which is quite accurate (erroneous by 0.58 (0 ) [ →1] eV). It gradually bends downwards after the HOMO orbital is half occupied. IPHF yields the final (1)

derivative which is quite poor (erroneous by 3.24 eV). Ideally, the linearity condition demands 0 1 IP , HF [ →1] [ →1] . Hence, the difference ∆ straight defines, in a way, the divergence from ε HOMO = ε LUMO = IPHF ( 0 ) − IPHF (1) [→2] the linearity. Figure 2 indicates that E HF ( λ ) is more quadratic than linear. Hence IPHF provides

a better way to approximate IPs. Table 3 also summarizes the EAs calculated with the HF method using various approximations. [ →1] are just the opposite of the orbital energies for the Data in the column under the title of EAHF (0)

13,29 [ →1] 0 LUMOs of the N0-systems, i.e., EAHF ( 0 ) = − ε LUMO , corresponding to the Koopmans’ theorem, [ →1] while data in the column under the title of EAHF are just the opposite of the orbital energies for (1) [ →1] 1 the HOMO of the (N0+1)-systems, i.e., EAHF . Ideally, the linearity condition demands (1) = − ε HOMO 0 1 EA , HF [ →1] [ →1] . Hence the difference ∆ straight may also be used as a way to ε LUMO = ε HOMO = EAHF ( 0 ) − EAHF (1)

quantify the divergence from the linearity. [ →1] Table 3 shows the statistics data. EAHF ( 0 ) gives data which are too small or even negative,

because orbitals are fixed in the N0-system while evaluating the total energy of the (N0+1)-system. As [ →1] compared to ∆HF, EAHF leads to MAD_1 and MAX_1 of 1.37 and 2.44 eV, respectively. On the (0) [ →1] other hand, EAHF overestimates EAs, because orbitals are fixed in the (N0+1)-system while (1)

22


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


[ →1] evaluating the total energy of the N0-system. As compared to ∆HF, EAHF leads to MAD_1 and (1) [ →1] [ →1] MAX_1 of 1.90 and ‒3.37 eV, respectively. Because of the opposite errors of EAHF and EAHF , (2) (1)

the two-point equation is a good approximation to ∆HF, yielding MAD_1 and MAX_1 of 0.27 and 0.47 eV, respectively. Note that ∆HF is known to be a poor approximation due to the lack of correlation effects, which is particularly important in stabilizing the anion systems. As compared to [ →1] the experimental data, ∆HF underestimates EAs (MAD = 1.46 eV). EAHF diverges even more (0)

from the experimental data (MAD = 2.75 eV) due to the accumulating effects of both relaxation [ →1] effects and correlation effects. EAHF , on the other hand, is the most satisfactory method at the (1)

level of HF for EA predictions, having MAD of 0.85 eV. Recall that it is well-established in density functional theory that the HOMO energy of the (N0+1)-system presents a good approximation of EA of the corresponding N0-system. Figure 2 shows that the HF energy curve against fractional charges. From N0 = 8 to the right, electron flows into the O atom, corresponds to the electron attachment process. The curve is clearly [ →1] in a concave manner, which gives too small initial derivatives. In fact, EA predicted by EAHF ( 0 ) is

negative (‒2.64 eV), as oppose to the experimental value of 1.46 eV.65,66 The final derivative is [ →1] usually satisfactory. EA predicted by EAHF is 1.73 eV, in fair agreement with the experimental (1)

value.

D. Ionization Potential and Electron Affinity Calculated with the MP2 method Table 4 summarizes the IPs calculated with the MP2 method using various approximations. Immediately, one finds that the quality of IP calculations are much improved at all levels of MP2 approximations as compared to those for HF. ∆MP2 yields MAD of only 0.25 eV, as compared to the 23



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

Page 24 of 42

I ,[ →1] III ,[ →1] experimental values. If the one-point equation is employed, IPMP and IPMP lead to 2(1) 2(1)

MAD_2/MAX_2 of 0.35/1.79 and 0.38/1.55 eV, respectively. These performances are ~0.2/~1.5 eV 28,29,45,52,53 I ,[ →1] III ,[ →1] better than those from IPMP the initial derivative, 2(0) and IPMP 2(0) . In the literature, I ,[ →1] II ,[ →1] II ,[ →1] lead to comparable MAD_2 around 0.5 eV, IPMP 2(0) , is most commonly used. IPMP 2(0) and IPMP 2(1)

whereas MAX_2 is 1.13 eV smaller for the former than for the latter. If the two-point equation is employed, the quality of the calculations is, in general, similar to that from the one-point equation. III ,[ → 2] yields MAD_2 of 0.42 eV, whose MAX_2 is reduced to 0.81 eV. Note that if more data IPMP 2

points are used in the integration approach, only L = III can faithfully reproduce the ∆MP2 results due to its nature of accurate analytical MP2 energy gradient with respect to the orbital occupation numbers. In fact, ∆MP2 sets up a reference with which all levels of other MP2 approximations can be [ →1] [ →1] compared. With the HF method, we see from Table 3 that IPs evaluated at IPHF and IPHF are (0 ) (1)

too large and small, respectively, as compared to ∆HF. Such trends have been compensated by the opposite contributions of the MP2 correlation energy derivatives for λ = 0 and λ = 1 (c.f. Table 1), [ →1] [ →1] leading to similar performances for IPMP and IPMP as compared to ∆MP2 shown in Table 4. 2( 0 ) 2 (1)

The two-point equation usually leads to a smaller error, as compared to the one-point equation. IP , MP 2 [ →1] [ →1] From the perspective of fractional charges, we compare the calculated ∆ straight = IPMP 2 ( 0 ) − IPMP 2 (1) IP IP , HF [ →1] [ →1] with ∆ straight . For a better theory, a smaller ∆ straight is expected according to the = IPHF ( 0 ) − IPHF (1)

IP linearity condition. We depict ∆ straight at different levels of theory in Figure 3(a). As demonstrated in

IP , MP 2 the figure, ∆ straight is nearly zero for atoms, in agreement with the observation that the MP2

fractional charge curve for total energy is nearly linear (c.f. Figure 2 for the O atom). The tendency IP for increasing ∆ straight from Li to F at the HF level is largely eliminated at the MP2 level, although

24


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


sizable errors still exist for the F atom at MP2. MP2 is less satisfactory in molecular species in IP contrast to its performance in the atomic species. This is particularly the case for F2. ∆ straight

amounts to ‒3.08 eV for MP2 at L = III, which is ‒1.91 and ‒2.58 eV, respectively, at L = I and II. As L = III is rigorous at the MP2 level, its values reflect the intrinsic deviations of MP2 from the IP , benefiting from error cancellations. linearity. It is possible that L = I, and II lead to smaller ∆ straight

The EA results calculated with the MP2 method using various approximations are presented in Table 5, which can be compared with the HF results listed in Table 3. While ∆HF is generally faulty for EA predictions, e.g. erroneously predicting an unstable O anion state (see also Figure 2), ∆MP2 [ →1] corrects such qualitative errors, reducing MAD_2 from 1.46 eV for ∆HF to 0.42 eV. EAHF is (0)

practically useless (MAD_2 = 2.75 eV), whereas adding correlation contributions removes the major part of errors, reducing MAD_2 to 0.56, 0.46 and 0.60 eV, respectively, for MP2 at the levels of L = I, [ →1] II, and III. EAHF seems to work reasonably well (MAD_2 = 0.85 eV), whereas adding correlation (1)

contributions destroys the error cancellation, increasing MAD_2 to 1.20, 0.94 and 1.36 eV, EA respectively. Figure 3 (b) depicts ∆ straight at various levels of approximations. The improvement is

EA most evident for atoms, where ∆ straight at the MP2 level approaches zero. Deviations from linearity

can still be sizable for electronegative atom like F and the molecular species, where higher order EA could be smaller for L = I, II than for L = III. correlations are desired.76 As noted in Figure 3, ∆ straight

Nonetheless, the derivatives are only exact for L = III, which reflect the intrinsic errors associated with the MP2 method.

IV. CONCLUDING REMARKS In the present work, we have derived the formalism to calculate the full MP2 energy derivatives with 25



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

Page 26 of 42

respect to the orbital occupation numbers, which is implemented and solved through the CP-HF equations. An integration approach is employed to calculate IPs, and EAs using one-point and two-point equations. The present method is an extension of the D-∆MBPT(2) method by Beste, et al.,46 which, in turn, is an extension of the single-point method of Cohen et al.45 from the perspective of fractional charges. While relaxation effects were included only at the HF level in the previous methods, such effects are fully taken into account in the present method up to the MP2 level. Numerical tests have been carried out for a selective set of atoms and molecules, which verifies the accuracy and usefulness of the present methodology. Comparison with the HF results signifies the importance to include the correlation effects. The formalism derived in the present work for calculating the full MP2 derivatives applies only to the integer points for λ = 0 and λ = 1. It is important to extend it also to consider the fractional occupations.80 Density functional theory has become the leading method in studying the IPs and EAs in materials sciences. Extension of the present work to the conventional functionals (i.e. local density approximations, generalized gradient approximations, hybrid functionals) and the recently developed doubly hybrid functionals is an important future direction.

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

REFERENCES

26


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


1. Parr, R.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. 2. Chattaraj, P. K., Ed. Chemical Reactivity Theory. A Density Functional View; CRC Press, Taylor & Francis Group: New York, 2009. 3. Pearson, R. G. Chemical Hardness; Wiley-VCH, Weinheim: Germany, 1997. 4. Geerlings, P.; De Proft, F.; Langenaeker, W. Chem. Rev. 2003, 103, 1793–1874. 5. Mulliken, R. S. J. Chem. Phys. 1934, 2, 782–293. 6. Ayers, P. W. J. Math. Chem. 2008, 43, 285–303. 7. Parr, R. G.; Pearson, R. G. J. Am. Chem. Soc. 1983, 105, 7512–7516. 8. Parr, R. G.; Szentpaly, L. V.; Liu, S. J. Am. Chem. Soc. 1999, 121, 1922–1924. 9. Bagus, P. S. Phys. Rev. 1965, 139, A619–A634. 10. Guest, M. F.; Saunders, V. R. Mol. Phys. 1975, 29, 873–884. 11. Chong, D. P.; Herring, F. G.; McWilliams, D. J. Chem. Phys. 1974, 61, 78–84. 12. Bacskay, G. B. Chem. Phys. 1977, 26, 47–57. 13. Koopmans, T. Physica 1933, 1, 104–113. 14. Pickup, B. T.; Goscinski, O. Mol. Phys. 1973, 26, 1013–1035. 15. Jørgensen, P. Annu. Rev. Phys. Chem. 1975, 26, 359–380. 16. Cederbaum, L. S.; Domcke, W. Adv. Chem. Phys. 1977, 36, 205–344. 17. Niessen von, W.; Schirmer, J.; Cederbaum, L. S. Compt. Phys. Rep., 1984, 1, 57–125. 18. Bartlett, R. J. Modern Electronic Structure Theory; World Scientific: Singapore, 1995. 19. Stanton, J. F.; Bartlett, R. J. J. Chem. Phys. 1993, 98, 7029–7039. 20. Stanton, J. F.; Gauss, J. J. Chem. Phys. 1994, 101, 8938–8944. 27



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

21. Nooijen, M.; Snijders, J. G. Int. J. Quantum Chem. 1993, 48, 15–48. 22. Koch, H.; Jensen, H. A.; Jørgensen, P.; Helgaker, T. J. Chem. Phys. 1990, 93, 3345–3350. 23. Lindgren, I. Int. J. Quantum Chem. 1978, 14(S12), 33–58. 24. Haque, M.; Mukherjee, D. J. Chem. Phys. 1984, 80, 5058–5069. 25. Nakatsuji, H. Computational Chemistry-Reviews of Current Trends, vol. 2; Leszczynski, J. Ed.; World Scientific: Singapore, 1997. 26. Ortiz, J. V. “The electron propagator picture of molecular electronic structure,” in Computational Chemistry: Reviews of Current Trends; World Scientific: Singapore, 1997; pp 1–61. 27. Flores-Moreno, R.; Zakrzewski, V. G.; Ortiz, J. V. J. Chem. Phys. 2007, 127, 134106. 28. Linderberg, J.; Öhrn, Y. Propagators in Quantum Chemistry, 2nd ed.; Wiley, 2004. 29. Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; Dover, 1996. 30. Cade, P. E.; Sales, K. D.; Wahl, A. C. J. Chem. Phys. 1966, 44, 1973–2003. 31. Cederbaum, L. S.; Hohlneicher, G.; von Niessen, W. Chem. Phys. Lett. 1973, 18, 503–508. 32. Purvis, G. D.; Öhrn, Y. J. Chem. Phys. 1974, 60, 4063–4069. 33. Chong, D. P.; Herring, F. G.; McWilliams, D. J. Chem. Phys. 1974, 61, 3567–3570. 34. Nooijen, M.; Snijders, J. G. J. Chem. Phys. 1995, 102, 1681–1688. 35. Stanton, J. F.; Gauss, J. J. Chem. Phys. 1995, 103, 1064–1076. 36. Christiansen, O.; Koch, H.; Jørgensen, P. Chem. Phys. Lett. 1995, 243, 409–418. 37. Nakajima, T.; Nakatsuji, H. Chem. Phys. Lett. 1999, 300, 1–8. 38. Gwaltney, S. R.; Nooijen, M.; Bartlett, R. J. Chem. Phys. Lett. 1996, 248, 189–198. 39. Nooijen, M.; Bartlett, R. J. J. Chem. Phys. 1997, 106, 6441–6448. 28


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


40. Foresman, J. B.; Head-Gordon, M.; Pople, J. A.; Frisch, M. J. J. Phys. Chem. 1992, 96, 135–149. 41. Head-Gordon, M.; Rico, R. J.; Oumi, M.; Lee, T. J. Chem. Phys. Lett. 1994, 219, 21–29. 42. Head-Gordon, M.; Maurice, D.; Oumi, M. Chem. Phys. Lett. 1995, 246, 114–121. 43. Cioslowski, J.; Piskorz, P.; Liu, G. J. Chem. Phys. 1997, 107, 6804–6811. 44. Flores-Moreno, R.; Zakrzewski, V. G.; Ortiz, J. V. J. Chem. Phys. 2007, 127, 134106. 45. Cohen, A. J.; Mori-Sánchez, P.; Yang, W. J. Chem. Theory Comput. 2009, 5, 786–792. 46. Beste, A.; Vázquez-Mayagoitia, Á.; Ortiz, J. V. J. Chem. Phys. 2013, 138, 074101. 47. Giesbertz, K. J. H.; Baerends, E. J. J. Chem. Phys. 2010, 132, 194108. 48. Beste, A.; Harrison, R. J.; Yanai, T. J. Chem. Phys. 2006, 125, 074101. 49. Janak, J. F. Phys. Rev. B 1978, 18, 7165–7168. 50. Blaizot, J.-P.; Ripka, G. Quantum Theory of Finite Systems; The MIT Press: Cambridge, 1986. 51. Casida, M. Phys. Rev. B, 1999, 59, 4694–4698. 52. Pino, R.; Scuseria, G. E. J. Chem. Phys. 2004, 121, 2553–2557. 53. Ayala, P. Y.; Kudin, K. N.; Scuseria, G. E. J. Chem. Phys. 2001, 115, 9698–9707. 54. Pople, J. A.; Krishnan, R.; Schlegel, H. B.; Binkley, J. S. Int. J. Quantum Chem. 1979, 16(S13), 225–241. 55. Frisch, M. J.; Head-Gordon, M.; Pople, J. A. Chem. Phys. Lett. 1990, 166, 275–280. 56. Salter, E. A.; Trucks, G. W.; Fitzgerald, G.; Bartlett, R. J. Chem. Phys. Lett. 1987, 141, 61–70. 57. Cioslowski, J.; Ortiz, J. V. J. Chem. Phys. 1992, 96, 8379–8389. 58. Su, N. Q.; Zhang, I. Y.; Xu, X. J. Comput. Chem. 2013, 34, 1759–1774. 59. Handy, N. C.; Schaefer, H. F., III J. Chem. Phys. 1984, 81, 5031–5033. 60. Handy, N. C.; Amos, R. D.; Gaw, J. F.; Rice, J. E.; Sirnandiras, E. D. Chem. Phys. Lett. 1985, 29



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

120, 151–158. 61. Head-Gordon, M. Mol. Phys. 1999, 96, 673–679. 62. 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. 63. Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Phys. Rev. B 2008, 77, 115123. 64. Su, N. Q.; Yang, W.; Mori-Sánchez, P.; Xu, X. J. Phys. Chem. A 2014, 118, 9201–9211. 65. Zhang, I. Y.; Xu, X.; Goddard, W. A., III Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 4963–4968. 66. Zhang, I. Y.; Xu, X.; Jung, Y.; Goddard, W. A., III Proc. Natl. Acad. Sci. U.S.A 2011, 108, 19896–19900. 67. Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Pople, J. A. J. Chem. Phys. 1997, 106, 1063– 1079. 68. 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). 69. Lin, Y. S.; Tsai, C. W.; Li, G. D.; Chai, J. D. J. Chem. Phys. 2012, 136, 154109. 70. Huber, K. P.; Herzberg, G. Constants of Diatomic Molecules, Molecular Spectra and Molecular Structure, Vol. 4; Van Nostrand Reinhold: Princeton, 1979. 71. Coriani, S.; Marchesan, D.; Gauss, J.; Hättig, C.; Helgaker, T.; Jørgensen, P. J. Chem. Phys. 2005, 123, 184107. 72. Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007–1023. 73. Woon, D. E.; Dunning, T. H., Jr. J. Chem. Phys. 1993, 98, 1358–1371. 74. Perdew, J. P.; Parr, R. G.; Levy, M.; Balduz, J. L., Jr. Phys. Rev. Lett. 1982, 49, 1691–1694. 30


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


75. Yang, W.; Zhang, Y. K.; Ayers, P. W. Phys. Rev. Lett. 2000, 84, 5172–5175. 76. Steinmann, S. N.; Yang, W. J. Chem. Phys. 2013, 139, 074107. 77. Vydrov, O. A.; Scuseria, G. E.; Perdew, J. P. J. Chem. Phys. 2007, 126, 154109. 78. Sham, L. J.; Schlüter, M. Phys. Rev. Lett. 1983, 51, 1888–1891. 79. Chai, J.-D.; Chen, P.-T. Phys. Rev. Lett. 2013, 110, 033002. 80. Peng D.; Yang, W. J. Chem. Phys. 2013, 138, 184108.

31



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 caption

Table 1. Comparison of the MP2 correlation energy derivatives with respect to the orbital occupation number nkλ for IP calculations at different levels of approximations.a,b

Table 2. Comparison of the MP2 correlation energy derivatives with respect to the orbital occupation number ncλ for EA calculations at different levels of approximations.a,b

Table 3. Comparison with the experimentala and ∆HF resultsb for ionization potential (IP) and electron affinity (EA) calculations using the one-pointc,d and two-pointe equations.

Table 4. Comparison with the experimentala and ∆MP2 resultsb for ionization potential (IP) calculations using the one-pointc and two-pointd equations.

Table 5. Comparison with the experimentala and ∆MP2 resultsb for electron affinity (EA) calculations using the one-pointc and two-pointd equations.

32


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


Figure caption

Figure 1. Comparison of the MP2 correlation energy derivatives with respect to the orbital occupation number ntλ at different levels of approximations for L = I, II, and III. Here the differences of numerical and analytical derivatives are displayed, where the numerical derivatives are the references. The energy unit is in eV, where t = k for IP calculations, and t = c for EA calculations. (a) The initial derivatives corresponding to HOMO of the N0-electron system, (b) the final derivatives corresponding to LUMO of the (N0‒1)-electron system, (c) the initial derivatives corresponding to LUMO of the N0-electron system; (d) The final derivatives corresponding to HOMO of the (N0+1)-electron system.

Figure 2. Fractional charge behaviors for the oxygen atom with HF and MP2. The inset depicts the deviations from the corresponding linear interpolations of the methods. IP EA Figure 3. Comparison of the derivative differences, (a) ∆ straight and (b) ∆ straight , evaluated at the

initial state for λ = 0 and the final state for λ = 1 with methods of different levels of approximations.

33



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

Table 1. Comparison of the MP2 correlation energy derivatives with respect to the orbital occupation number nkλ for IP calculations at different levels of approximations.a,b λ = 0b Numc

Ic

λ = 1b IIc

IIIc

Numc

Ic

IIc

IIIc

Li

-0.03

-0.03

-0.04

-0.03

-0.04

-0.03

-0.04

-0.04

Be

-0.27

-0.56

-0.47

-0.27

-1.27

-1.28

-1.29

-1.27

B

0.50

0.22

0.18

0.50

-0.85

-0.89

-1.02

-0.85

C

0.84

0.61

0.55

0.84

-1.40

-1.39

-1.55

-1.40

N

1.26

1.09

1.01

1.26

-2.04

-1.98

-2.18

-2.04

O

1.25

1.17

1.07

1.26

-2.76

-2.72

-2.90

-2.76

F

2.11

2.13

1.92

2.11

-3.81

-3.68

-3.94

-3.81

F2

4.62

4.43

4.02

4.62

-2.24

-1.25

-2.34

-2.24

OH

1.95

1.85

1.54

1.95

-3.56

-3.44

-3.78

-3.56

NH2

1.42

1.18

0.91

1.42

-3.20

-3.11

-3.44

-3.20

CH3

1.26

0.83

0.64

1.26

-2.16

-2.23

-2.53

-2.16

CN

1.78

2.69

2.14

1.79

-0.90

-1.13

-2.05

-0.88

O2

5.05

4.12

3.72

5.05

-1.08

-0.83

-2.09

-1.08

MADd

-

0.30

0.41

0.00

-

0.16

0.30

0.00

d

-

0.93

1.33

0.00

-

-0.98

1.15

-0.02

MAX a

The data set is taken from Refs. 45, 63, 64. The unit is eV.

b

The MP2 correlation energy derivative with respect to the orbital occupation number is defined by

dEcL,MP 2 ( λ ) , where the orbital index k refers to the HOMO dnkλ

of the N0-electron system, which can be partially occupied. The initial state (λ=0) corresponds to the N0-electron system where nk0 = 1 , while the final state (λ=1)

corresponds to the (N0‒1)-electron system where n1k = 0 , indicating that the k-th orbital has become the LUMO. c

L indicates the level of approximation. L = Num stands for a finite difference calculation with

∆EcNum , MP 2 ( λ ) using ∆ n kλ = 0.001. L = I refers to the approximation ∆nkλ

where only explicit dependence of the MP2 correlation energy on the occupation number is considered as in Ref. 45. L = II refers to additional consideration of the explicit dependence of the HF eigenvalues on the occupation number as in Ref. 46. L = III is the full MP2 correlation energy derivative developed in the present work. Ideally, one has Num = III. d

MAD and MAX refer to the mean absolute deviation and the maximum deviation (Ref. – Calc.), respectively, from the finite difference results.

34


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


Table 2. Comparison of the MP2 correlation energy derivatives with respect to the orbital occupation number ncλ for EA calculations at different levels of approximations.a,b λ = 0b Numc

Ic

λ = 1b IIc

IIIc

Numc

Ic

IIc

IIIc

Li

-0.52

-0.52

-0.52

-0.52

-0.01

-0.24

-0.15

-0.01

Be

-0.43

-0.48

-0.55

-0.43

0.22

0.03

0.05

0.22

B

-0.97

-0.96

-1.08

-0.96

0.60

0.38

0.36

0.60

C

-1.66

-1.60

-1.78

-1.66

1.14

0.96

0.89

1.14

N

-2.03

-2.00

-2.15

-2.03

0.77

0.82

0.76

0.78

O

-3.19

-3.05

-3.29

-3.19

1.95

2.03

1.81

1.95

F

-4.53

-4.24

-4.60

-4.53

3.18

3.37

2.99

3.19

F2

-3.49

-2.68

-3.68

-3.49

4.98

4.83

4.28

4.99

OH

-3.99

-3.78

-4.22

-3.99

2.86

2.78

2.31

2.86

NH2

-3.13

-3.04

-3.43

-3.13

1.94

1.67

1.29

1.94

CH3

-2.08

-2.05

-2.28

-2.08

0.72

0.40

0.23

0.72

CN

-4.23

-2.89

-3.27

-4.23

2.50

1.80

1.39

2.50

O2

-2.77

-2.29

-3.06

-2.77

2.58

2.47

2.03

2.58

MADd

-

0.28

0.22

0.00

-

0.21

0.40

0.00

d

-

-1.34

-0.96

-0.01

-

0.70

1.11

-0.01

MAX a

The data set is taken from Refs. 45, 63, 64. The unit is eV.

b

The MP2 correlation energy derivative with respect to the orbital occupation number is defined by

dEcL,MP 2 ( λ ) , where the orbital index c refers to the LUMO dncλ

of the N0-electron system, which can be partially occupied. The initial state (λ=0) corresponds to the N0-electron system where nc0 = 0 , while the final state (λ=1) corresponds to the (N0+1)-electron system where n1c = 1 , indicating that the c-th orbital has become the HOMO . c

L indicates the level of approximation. L = Num stands for a finite difference calculation with

∆EcNum , MP 2 ( λ ) using ∆ ncλ = 0.001. L = I refers to the approximation ∆ncλ

where only explicit dependence of the MP2 correlation energy on the occupation number is considered as in Ref. 45. L = II refers to additional consideration of the explicit dependence of the HF eigenvalues on the occupation number as in Ref. 46. L = III is the full MP2 correlation energy derivative developed in the present work. Ideally, one has Num = III. d

MAD and MAX refer to the mean absolute deviation and the maximum deviation (Ref. – Calc.), respectively, from the finite difference results.

35



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

Table 3. Comparison with the experimentala and ∆HF resultsb for ionization potential (IP) and electron affinity (EA) calculations using the one-pointc,d and two-pointe equations. Expt-IPa

∆HF-IPb

[ →1] IPHF ( 0)

c

[ →1] IPHF (1)

d

[→2] IPHF

e

Expt-EAa

∆HF-EAb

[ →1] EAHF (0 )

c

[ →1] EAHF (1)

d

[ → 2] EAHF

Li

5.39

5.34

5.34

5.34

5.34

0.62

-0.17

-0.29

0.26

-0.01

Be

9.32

8.04

8.42

7.78

8.10

-0.36

-0.92

-1.19

-0.57

-0.88

B

8.30

8.04

8.67

7.51

8.09

0.28

-0.43

-1.09

0.53

-0.28

C

11.26

10.80

11.94

9.85

10.90

1.26

0.33

-0.78

1.84

0.53

N

14.54

13.89

15.52

12.52

14.02

-0.22

-2.26

-3.37

-0.46

-1.91

O

13.61

12.02

14.19

10.37

12.28

1.46

-0.88

-2.64

1.73

-0.46

F

17.42

15.65

18.47

13.41

15.94

3.40

0.90

-1.54

4.27

1.36

F2

15.70

16.13

18.13

14.36

16.24

0.42

0.03

-2.37

2.78

0.21

OH

13.01

11.36

13.95

9.39

11.67

1.83

-0.59

-2.61

2.36

-0.12

NH2

11.14

10.45

12.59

8.79

10.69

0.77

-1.39

-2.89

0.91

-0.99

CH3

9.84

8.98

10.47

7.75

9.11

0.08

-1.96

-2.83

-0.43

-1.63

CN

13.60

12.47

14.65

10.96

12.55

3.86

2.83

0.85

5.11

2.98

O2

12.30

13.38

15.18

11.77

13.48

-0.08

-1.21

-2.80

0.73

-1.04

-

-

1.58

1.29

0.15

-

-

1.37

1.90

0.27

MAD_1f f

-

-

-2.82

2.24

0.31

-

-

2.44

-3.37

0.47

MAD_2g

-

0.91

1.08

1.97

0.80

-

1.46

2.75

0.85

1.26

g

-

1.77

-2.88

4.01

1.48

-

2.50

4.94

-3.53

2.04

MAX_1

MAX_2 a

e

The data set is taken from Refs. 45, 63, 64. The unit is eV. The experimental data taken from Refs. 67, 68, except EAs for Be, N, F2 and O2, which are CCSD(T)

data from Ref. 69. b

In the ∆HF method, ionization energy is calculated as IP = EHF(N0–1) – EHF(N0), while electron affinity is calculated as EA = EHF(N0) – EHF(N0+1).

c

IP or EA is calculated up to the first order using −

dE HF ( λ ) = −ε t0 , which corresponds to the HOMO energy of the N0-electron system for IP or the LUMO dntλ λ =0

energy of the N0-electron system for EA. d

IP or EA is calculated up to the first order using −

dE HF ( λ ) = −ε t1 , which corresponds to the LUMO energy of the (N0-1)-electron system for IP or HOMO dntλ λ =1

energy of the (N0+1)-electron system for EA. 1 0 ( ε t + ε t1 ) . 2

e

IP or EA is calculated up to the second order using −

f

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

g

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

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


Table 4. Comparison with the experimentala and ∆MP2 resultsb for ionization potential (IP) calculations using the one-pointc and two-pointd equations. Li

Expt.a

∆MP2b

I ,[ →1] IPMP 2(0)

5.39

5.38

5.37

c

I ,[ →1] IPMP 2(1)

5.38

c

I ,[ → 2] IPMP 2

5.38

d

II ,[ →1] IPMP 2(0)

5.38

c

II ,[ →1] IPMP 2(1)

c

5.38

II ,[ → 2] IPMP 2

5.38

d

III ,[ →1] IPMP 2(0)

c

5.37

III ,[ →1] IPMP 2(1)

5.38

c

III ,[ → 2 ] IPMP 2

5.38

Be

9.32

8.88

8.98

9.06

9.02

8.89

9.06

8.98

8.69

9.05

8.87

B

8.30

8.31

8.45

8.40

8.42

8.49

8.52

8.51

8.17

8.35

8.26

C

11.26

11.30

11.34

11.24

11.29

11.39

11.40

11.39

11.10

11.25

11.18

N

14.54

14.63

14.44

14.50

14.47

14.52

14.70

14.61

14.27

14.57

14.42

O

13.61

13.42

13.02

13.09

13.05

13.11

13.27

13.19

12.93

13.13

13.03

F

17.42

17.37

16.34

17.09

16.72

16.55

17.35

16.95

16.36

17.22

16.79

F2

15.70

15.44

13.70

15.61

14.65

14.11

16.70

15.40

13.51

16.59

15.05

OH

13.01

13.06

12.10

12.84

12.47

12.42

13.18

12.80

12.00

12.95

12.48

NH2

11.14

11.99

11.41

11.90

11.66

11.68

12.23

11.96

11.17

11.98

11.58

CH3

9.84

9.76

9.64

9.98

9.81

9.83

10.28

10.05

9.21

9.91

9.56

CN

13.60

12.93

14.54

15.39

14.97

14.83

16.32

15.57

13.67

15.15

14.41

O2

12.30

11.79

11.06

12.60

11.83

11.46

13.86

12.66

10.13

12.85

11.49

MAD_1e

-

-

0.59

0.39

0.41

0.47

0.64

0.40

0.71

0.42

0.36

MAX_1

e

-

-

1.74

2.46

2.04

1.90

3.39

2.64

1.93

2.22

1.48

MAD_2

f

-

0.25

0.61

0.35

0.44

0.54

0.63

0.42

0.69

0.38

0.42

MAX_2f

-

0.85

2.00

1.79

1.37

1.59

2.72

1.97

2.19

1.55

0.81

a

The data set is taken from Refs. 45, 63, 64. The unit is eV. The experimental data are taken from Refs. 67, 68.

b

In the ∆MP2 method, ionization energy is calculated as IP = EMP2(N0–1) – EMP2(N0).

c

IP is calculated up to the first order using −

d

L dE MP 2 (λ ) at λ = 0 or λ = 1 under different levels of approximations (L = I as in Ref. 45, II as in Ref. 46, and III dnkλ

developed in the present work). d

L L  dEMP 1  dE MP 2 ( λ ) 2 (λ ) IP is calculated up to the second order using −  . + λ   2  dnk dnkλ λ =0 λ =1 

e

MAD_1 and MAX_1 refer to the mean absolute deviation and the maximum absolute deviation, respectively, from the ∆MP2 results.

f

MAD_2 and MAX_2 refer to the mean absolute deviation and the maximum absolute deviation, respectively, from the experimental results.

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

Page 38 of 42

Table 5. Comparison with the experimentala and ∆MP2 resultsb for electron affinity (EA) calculations using the one-pointc and two-pointd equations. Expt.a

∆MP2b

I ,[ →1] EAMP 2(0)

Li

0.62

0.32

0.23

0.51

0.37

0.23

0.42

0.32

0.22

0.27

0.25

Be

-0.36

-0.75

-0.71

-0.61

-0.66

-0.64

-0.62

-0.63

-0.76

-0.79

-0.78

B

0.28

0.04

-0.13

0.15

0.01

0.00

0.17

0.08

-0.12

-0.07

-0.10

C

1.26

1.09

0.82

0.88

0.85

1.00

0.95

0.98

0.88

0.70

0.79

N

-0.22

-0.87

-1.37

-1.29

-1.33

-1.22

-1.22

-1.22

-1.34

-1.24

-1.29

O

1.46

0.95

0.40

-0.29

0.05

0.65

-0.08

0.28

0.55

-0.22

0.16

F

3.40

3.14

2.70

0.90

1.80

3.06

1.28

2.17

2.99

1.08

2.03

c

I ,[ →1] EAMP 2(1)

c

I ,[ → 2 ] EAMP 2

d

II ,[ →1] EAMP 2(0)

c

II ,[ →1] EAMP 2(1)

c

EAMIIP,[2→ 2 ]

d

III ,[ →1] EAMP 2( 0)

c

III ,[ →1] EAMP 2(1)

c

EAMIIIP,[2→ 2 ]

F2

0.42

0.04

0.31

-2.05

-0.87

1.31

-1.50

-0.10

1.13

-2.20

-0.54

OH

1.83

1.51

1.17

-0.42

0.37

1.61

0.06

0.83

1.39

-0.50

0.44

NH2

0.77

0.40

0.14

-0.76

-0.31

0.54

-0.37

0.08

0.24

-1.03

-0.40

CH3

0.08

-0.47

-0.78

-0.83

-0.81

-0.55

-0.66

-0.60

-0.75

-1.15

-0.95

CN

3.86

4.80

3.74

3.31

3.52

4.12

3.72

3.92

5.08

2.60

3.84

O2

-0.08

-0.49

-0.51

-1.74

-1.13

0.26

-1.30

-0.52

-0.03

-1.85

-0.94

MAD_1e

-

-

0.33

0.99

0.62

0.32

0.74

0.34

0.30

1.08

0.55

MAX_1

e

-

-

1.06

2.24

1.34

1.27

1.86

0.97

1.09

2.24

1.11

MAD_2

f

-

0.42

0.56

1.20

0.88

0.46

0.96

0.61

0.60

1.36

0.83

MAX_2f

-

0.94

1.15

2.50

1.60

1.00

2.12

1.23

1.22

2.62

1.39

a

d

The data set is taken from Refs. 45, 63, 64. The unit is eV. The experimental data taken from Refs. 67, 68, except EAs for Be, N, F2 and O2, which are CCSD(T)

data from Ref. 69. b

In the ∆MP2 method, electron affinity is calculated as EA = EMP2(N0) – EMP2(N0+1).

c

EA is calculated up to the first order using −

L dEMP 2 (λ ) at λ = 0 or λ = 1 under different levels of approximations (L = I as in Ref. 45, II as in Ref. 46, and III dncλ

developed in the present work). d

L L  dE MP 1  dEMP 2 ( λ ) 2 (λ ) EA is calculated up to the second order using −  . + λ  2  dncλ dn c λ =0 λ =1 

e

MAD_1 and MAX_1 refer to the mean absolute deviation and the maximum absolute deviation, respectively, from the ∆MP2 results.

f

MAD_2 and MAX_2 refer to the mean absolute deviation and the maximum absolute deviation, respectively, from the experimental results.

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 1. Comparison of the MP2 correlation energy derivatives with respect to the orbital occupation number ntλ at different levels of approximations for L = I, II, and III. Here the differences of numerical and analytical derivatives are displayed, where the numerical derivatives are the references. The energy unit is in eV, where t = k for IP calculations, and t = c for EA calculations. (a) The initial derivatives corresponding to HOMO of the N0-electron system, (b) the final derivatives corresponding to LUMO of the (N0‒1)-electron system, (c) the initial derivatives corresponding to LUMO of the N0-electron system; (d) The final derivatives corresponding to HOMO of the (N0+1)-electron system.

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 2. Fractional charge behaviors for the oxygen atom with HF and MP2. The inset depicts the deviations from the corresponding linear interpolations of the methods.

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


IP EA Figure 3. Comparison of the derivative differences, (a) ∆ straight and (b) ∆ straight , evaluated at the

initial state for λ = 0 and the final state for λ = 1 with methods of different levels of approximations.

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

TOC

The integration approach at the second-order perturbation theory: 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

42


Applications to Ionization Potential and Electron Affinity Calcula

Recommend Documents