Application of self-consistent $α$ method to improve the performance of model exchange potentials

0.13 MB·English

by  Valentin V. Karasiev

Application of self-consistent α method to improve the performance of model
exchange potentials
Valentin V. Karasiev,1,∗ Eduardo V. Luden˜a,1 and Art¨em E. Masunov2
1Centro de Qu´ımica, Instituto Venezolano de Investigaciones Cient´ıficas, Apartado 21827,
Caracas 1020-A, Venezuela
2NanoScience Technology Center, Department of Chemistry and Department of Physics,
University of Central Florida, 12424 Research Parkway, Orlando, Florida 32826
(Dated: 17 October2009)
2
Self interaction error remains an impotrant problem in density functional theory. A number of
1
approximations to exact exchange aimed to correct for this error while retainining computational
0
efficiencyhadbeensuggestedrecently. Wepresentacriticalcomparisonbetweenmodelexchangepo-
2
tentialsgeneratedthroughtheapplicationoftheasymptotically-adjustedself-consistentα,AASCα,
n method and BJ effective exchange potential advanced in [A.D. Becke and E.R. Johnson, J. Chem.
a
Phys. 124, 221101 (2006)] and [V.N. Staroverov, J. Chem. Phys. 129, 134103 (2008)]. In particu-
J
lar we discuss their compliance with coordinate-scaling, virial and functional derivativeconditions.
9 We discuss the application of the AASCα method to generate the AA-BJ potential. A numerical
2 comparison is carried out through the implementation of a fully-numerical diatomic molecule code

yielding molecular virial energies and ionization potentials approximated by the energies of the
]
h HOMO orbitals. It is shown that some of the shortcomings of these model potentials, such as the
p non-compliancewiththeLevy-Perdewvirialrelation,maybeeliminatedbymultiplyingtheresponse
- termbyanorbital-dependentfunctionalα,whichcanbesimplifiedtoaconstantdeterminedduring
m theself-consistent procedure (self-consistent α).
e
h PACSnumbers:
c
.
s I. INTRODUCTION tial related to an approximate (and arbitrary) exchange
c
i functional Ex0[ρ] = d3rρ(r)ǫ0x([ρ];r) yielding the local
ys Animportantproblemarisingintheapplicationofthe exchange potential vRx0(r) = δEx0[ρ]/δρ; (however, more
h Kohn-Sham equations is that of the construction of the generally,theresponsetermv˜r0esp alsomaybemodelled).
p local exact exchange potential. In principle, this poten- In Eq. (1), αx[{ψi}] is the functional
[ tial can be exactly calculated following the procedure
1 outlined in the OEP method. [1–9] However, in prac- α [{ψ }]= Ex[{ψi}]−ExLP vS([ρ]),ρ (2)
v tice, there are a number of difficult numerical imped- x i ExLP v˜r0esp([(cid:2)ρ]),ρ (cid:3)
8 iments that bar the way to the realization of this ex- (cid:2) (cid:3)
3 actapproach[10–13]inadditiontosomedeeperproblems where Ex[{ψi}] is the exact orbital expression for
0 stemmingfromgeneralinstabilitiesofKohn-Shampoten- the exchange energy, and ExLP v([ρ];r),ρ(r) =
6 tials in finite dimensional subspaces.[13] d3rv(r)[3ρ(r)+r·∇ρ(r)])istheLevy(cid:2)-Perdewexpre(cid:3)ssion
.
1 For this reason, particularly in recent years, much at- Rfor the exchange energy functional corresponding to the
0 tention has been devoted to the development of alterna- potential v(r) [25, 26]. The local potential given by Eq.
2 tive local-exchange potentials which are simple to apply (1)iscalled“asymptoticallyadjusted” becausewhilethe
1 but which at the same time yield sufficiently accurate first term guarantees the correct asymptotic behavior of
:
v results.[3, 14–23] −1/r for large r the second term does not contribute in
Xi A few years ago, we introduced such an approach, the asymptotic region.
which we called the “asymptotically-adjusted self- The second characteristic follows from Eq. (2):
r
a consistentα” (AA-SCα)method. Althoughthedetailed E [{ψ }] = E v ([ρ]),ρ +α [{ψ }]E v˜0 ([ρ]),ρ
theoretical justification for this method is presented in x i xLP S x i xLP resp
Ref. [24], we comment here on two of its characteristics. ≡ EAAS(cid:2)Cα[{ψ }] (cid:3) (cid:2) (3(cid:3))
x i
The first is that in the AA-SCα method the following
The variationalderivative of this functionalwith respect
model potential is postulated:
to the Kohn-Sham orbital yields:
vAASCα(r)=v (r)+α [{ψ }]v˜0 (r) (1)
x S x i resp δEAASCα[{ψ }]
x i = vAASCα(r)+∆v (r) ψ (r) (4)
where vS(r) is the local Slater potential and where δψj (cid:18) x x (cid:19) j
v˜r0esp(r) = vx0(r) −2ǫ0x([ρ];r) is a local response poten- where ∆v (r)ψ (r) = [v (r)−v (r)−b α v˜0 (r)]ψ (r)
x j xj S x resp j
is a non-local correction to the vAASCα(r) potential. It
x
has beenbshown [24] tbhat the contribution of ∆v (r)
x
∗Electronicaddress: [email protected]fl.edu to the energy is ∆Ex[{ψi}] = 0 where ∆Ex[{ψi}] =
b
2
E [{ψ }] − E [v ] − α [{ψ }]E [v˜0 ]. Hence, the E v ([ρ ]),ρ +α E v˜0 ([ρ ]),ρ andthe fact
x i xLP S x i xLP resp xLP S λ λ x xLP resp λ λ
omission of ∆vx(r) in the Kohn-Sham equation only af- that (cid:2)δExAASCα[ρλ(cid:3)] = λδExA(cid:2)ASCα[ρ], and (cid:3)following the
fectsthequalityoftheconvergedorbitalsbutitdoesnot same arguments as in Ref. [28], it can be readily shown
change the expbression for the energy. that vAASCα(r) also satisfies the scaling property.
x
The AASCα method has been applied previously to In the case of the potentials introduced by Becke and
improve the potentials and energies of several DFT ex- Johnson [20] and Staroverov [23], it is clear that they
change functionals. In particular, we have analyzed the satisfy the scaling requirement. However, the fact that
improvements brought about by this method with re- thesepotentialsdonotsatisfytheLevy-Perdewcondition
spect to the LDA and PW91 functionals and proposed implies that the virialrelationin the Born-Oppenheimer
two models for v˜r0esp term [24]. approximation[29]isnotsatisfiedeither(seeRef. [30]for
However, quite recently a very simple potential de- details). Moreover,as a consequenceofthis, “thereis no
noted as the BJ effective exchange potential has been unique choice of functional for the evaluation of the ex-
proposed in a heuristic way by Becke and Johnson[20]. changeenergy”forstructureoptimizationortotalenergy
ThispotentialcontainsaSlatertermplusalocalresponse comparison (see Ref. [31], where the Becke-Johnson po-
one. ThisworkhasbeenextendedbyStaroverov[23]who tentialwasemployedforbandgapcalculationsinsolids).
hasadvancedafamilyofmodelpotentialswhichhavethe Also in the case of the BJ and other model potentials
form of a Slater potential plus a response term which is discussed in[20, 23], the third requirement from [27, 28],
modeled. In all these cases the energy is evaluated using namely,thatthemodelpotentialmustcorrespondtothe
the exactexpressionEx[{ψi}] for the exchange function- functionalderivativeofthemodelfunctionalwithrespect
als constructedfromthe N occupiedorbitalswhich have to the density is not satisfied.
self-consistentlyconvergedforthegivenmodelpotential. We denote as in Ref. [23] by E the total energy
conv
In the present work, due to the similarities between computed with the exact exchange expression E [{ψ }]
x i
this family ofpotentials andthe AASCα one, wemake a using the converged orbitals and by E the corre-
vir
criticalcomparisonbetweenthesepotentials. Wealsoap- spondingonewhichincludestheLevy-Perdewexpression
ply the AASCα method using the model exchange func- E [v ,ρ] for exchange. Since for the BJ and related
xLP x
tional associated with the BJ potential to generate the potentials E [{ψ }] 6= E [v ,ρ] the total energy val-
x i xLP x
AA-BJone. Thesystemschosenforthe presentcompar- ues E and E are also different. The claims that
conv vir
isonaresomeselecteddiatomicmolecules. Weshowthat “(E −E ) gives an indication of how close v is to
vir conv xσ
applicationofthe AASCα methoddoesindeedbringim- the exact functional derivative”[23] or that it serves as
provements, albeit slight, on both the BJ energies and anindicationoftheaccuracyofthecalculationitself(see
ionization potentials. Also, we show that it yields in- Ref. [32]) do not seemto holdgenerally,as there are ap-
ternucleardistancesthatareinexcellentagreementwith proximate KS exchange potentials (such as the AASCα
the exact Kohn-Sham x-only results. one, for example), which while differing from the exact
one, yield, nonetheless, E =E .
vir conv
In the present article, in order to compare the results
II. COMPARISON WITH THE BJ AND
obtained using the AASCα model with those coming
RELATED MODELS
from the BJ and related potentials, we provide numeri-
calvalues,inparticular,fordiatomicmolecules. Inorder
There are three formal conditions that the exact op- to carry out this numerical comparison, the BJ effective
timized effective potential (OEP) for exchange must model potentials proposed in [20] and [23] were imple-
satisfy[27,28]. These are: the variationalderivative con- mented in a fully-numerical diatomic molecule code [34]
dition, the virial relation, and the scaling requirement. (duetoitsnumericalinstabilitystemmingfrom“trouble-
Quite clearly,the AASCα localpotential vAASCα(r) is some ”, terms, a fact that was corroborated in our test
x
an approximate one and, hence, it does not satisfy all calculations, the gradient-corrected model [23] was not
three of these conditions. As it is shown in Eq. (4), implemented).
the variational derivative of EAASCα[{ψ }] with respect For completeness, we include some results obtained in
x i
to ψ yields the potential vAASCα(r)+∆v (r). Let us the context of the GLLB model proposed in Ref. [14],
i x x
emphasize, however, that ∆v (r) does not contribute the localized Hartree-Fock (LHF) model [15, 33], and
x
to the exchange energy, omission of thisbterm in the thecommonenergydenominatorapproximation(CEDA)
Kohn-Sham equation makes tbhe local potential vAASCα [35]. We compare these results with those of the BJ and
x
the approximate one. With respect to the virial rela- related potentials as well as with the AA-BJ ones ob-
tion, it follows from the definition of vAASCα(r) given tained by applying the AASCα model to the BJ func-
x
by Eq. (1) that vAASCα(r) satisfies by construction tional. It is shown that the AA-BJ results are only
x
the Levy-Perdew virial condition. Also, it is easy to slightly improved due to this application. We also in-
show using the fact that both E [ρ ] = λE [ρ] and clude some previous AASCα model results obtained for
x λ x
E v([ρ ]),ρ = λE v([ρ]),ρ that α [{ψ }], as the PW91 andGLLB modelfunctionals (in particularof
xLP λ λ xLP x i
define(cid:2)d by Eq. (cid:3)(2), is invar(cid:2)iant und(cid:3)er coordinate scal- the AA-PW91 and AA-m2 types, see Ref. [24]). These
ing (does not depend on λ). Then using EAASCα[ρ ] = resultsshowthatallthesemodelsareclosertoEXXthan
x λ
3
the BJ or AA-BJ ones in the case ofdiatomic molecules. Kohn-Sham potential. The asymptotically-adjusted po-
tentials based on the DFT approximation for the re-
sponseterm(the AA-PW91)areinvariantw.r.t. unitary
III. RESULTS AND DISCUSSION transformation of orbitals (as is also the BJ potential);
however, the AA-m2 and GLLB ones are not. The ad-
vantages of the previously proposed AA-PW91, AA-m2
Table I shows the Hartree-Fock total energies and
andGLLBmodels(andoftheCEDAandLHFmethods)
the energy differences for KS-x-only methods calcu-
are the following: (i) For these models E =E (no-
lated when the virial relation is used (except for the vir conv
tations from [23]); these models constitute examples of
“BJ(conv)” column). The energy differences of the ap-
situations where although (E −E )=0, the poten-
proximate methods (eight last columns) should be com- vir conv
tialisstillonlyapproximate;(ii)theenergiesobtainedby
pared to the exact exchange (EXX) values shown in the
theAA,GLLBmodelsandbytheCEDAandLHFmeth-
second column or to the KS(x-only) values obtained by
ods satisfy the variational principle: EHF ≤ EOEP ≤
the iterative procedure described in Ref. [36] (shown in
Eapprox, which is not the case for the BJ model; (iii) the
the third column), which are very close to the EXX val- vir
virial energy for the AA and for the GLLB models is an
ues.
excellent approximation to the OEP/EXX energy. This
All approximate methods shown in Table I, except for
isnotthecaseforthemodelpotentialofRefs. [20,23],in
the BJ potential, are exactfor the singlet state of a two-
spite of the fact that the BJ model is in excellent agree-
electron system (H molecule). The orbital-dependent
2
ment with the OEP exchange potential; (iv) the ioniza-
methodsfortheresponsepotentialterm(GLLB,AA-m2,
tion potential approximated by the HOMO energy is a
CEDA and LHF) provide very good approximations to
property entirely defined by the effective potential. The
theEXXenergies,theCEDAandLHFvaluesareonly 2
recentlyproposedBJmodel potentialfails to adequately
mHartrees higher than the corresponding EXX energies.
describe this property, in contrast with the AA-PW91,
The AA-PW91 potential, taken here as an example of
GLLB and AA-m2 models and with CEDA and LHF
an asymptotically-adjustedpotential where the response
methods, where the agreement with the EXX values is,
termismodeledbyaconventionalPW91DFTfunctional
in most cases, excellent.
was found to yield a good approximation to the EXX
Someoftheshortcomingsofthemodelpotentialsfrom
energy;the largestdifference is 39.4mHartees forthe F
2
Refs. [20,23]maybeeliminated,however,byscalingthe
molecule (as compared to the EXX difference which is
model response term by the AASCα method, as it was
8.6 mHartrees).
done in Ref. [24]. A rigorous variational justification for
The BJ differences still are significantly larger than
this type of scaling is given in Eqs. (12) through (15) of
those arising from all other approximate methods. The
Ref. [24]. By applying this procedure to the BJ model,
negative value in Table I shows that the corresponding
the new AA-BJ potential (for the spin-unpolarizedcase)
energy is lower than the HF value, i.e. the variational
reads
principle is not satisfied. Large errors in the total en-
ergy also affect significantly the calculated atomization vAA−BJ =v +α [{ψ }] τ[{ψ }]/ρ, (5)
energies and the predicted equilibrium geometries. x S x i i
p
The ionization potentials approximated by the high- whereτ is the kinetic energydensity. The self-consistent
est occupied molecular orbital (HOMO) energy are pre- constant α is defined by
x
sented in Table II. The situation is similar: the GLLB,
AA-m2, CEDA and LHF methods provide an excellent Ex[{ψi}]−ExLP[vS,ρ]
α [{ψ }]= . (6)
approximationtotheEXXvalues. TheAA-PW91values x i E [ τ[{ψ }]/ρ,ρ]
xLP i
arealsoveryclosetotheEXXforallmoleculespresented
p
inTableexceptforthe N2 andFH.TheBJpotentialun- The energies for the new AA-BJ model (we emphasize
derestimates the ionization potential by an amount of that E = E for AA-BJ) are presented in last col-
vir conv
∼30-50%. The BJ potential is calculated without shift umn of Table I. The AA-BJ energies are slightly closer
as suggested in [20]. By applying the shift, the HOMO to the EXX values than E values for the original BJ
conv
energies will be equal to the corresponding HF values model(BJ(conv)columninTableI).Moreover,modified
(but not to the OEP ones). AA-BJ model eliminates the ambiguity with regard to
thechoiceoffunctionalfortheevaluationoftheexchange
energy(conventionalorvirial). Theionizationpotentials
presented in Table II also are slightly improved in the
AA-BJ model as compared to the original BJ values.
In Table III we present some bond lengths values ob-
tained from severaldifferent x-only methods. We do not
include the diatomic molecule F as it does not bind at
2
The asymptotically-adjusted, GLLB, CEDA and LHF the level of an x-only approximation. The results show
exchange potentials, as well as the BJ model potential, that for H , FH and N , the AA-BJ bond lengths coin-
2 2
provideanexcellentapproximationtotheexactexchange cide up to three decimals with those obtained by means
4
TABLE I: Full-numerical Hartree-Fock (HF) total energies (in a.u.) and differences between KS-x-onlyand HF total energies
(in mHartrees) calculated at theexperimental geometries. aValuesare taken from Ref. [6]. bFrom Ref. [24]. cFrom Ref. [35] .
dFrom Ref. [33].
HF EXXa KS(x-only) AA-PW91b GLLBb AA-m2b CEDAc LHFd BJ(vir) BJ(conv) AA-BJ
H2 -1.1336 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -80.6 0.8 0.0
FH -100.0708 2.0 2.2 17.0 5.9 3.7 533.1 10.4 6.6
N2 -108.9931 5.2 5.7 30.4 11.4 8.3 7.7 7.3 239.4 9.9 9.5
CO -112.7909 5.1 5.6 31.1 12.2 9.1 7.6 7.2 323.0 12.6 11.5
F2 -198.7722 8.6 9.3 39.4 16.8 13.6 1124.9 23.1 18.1
TABLEII:IonizationpotentialapproximatedbythenegativeoftheHOMOenergies(ineV).aFromRef. [5]. bRef. [24]. cRef.
[6]. dN2→N+2(2Σg). eN2 →N+2(2Πu). fFrom Ref. [35] . gFrom Ref. [15].
HF EXXa AA-PW91b GLLBb AA-m2b CEDAf LHFg BJ AA-BJ
H2 16.2 16.2 16.2 16.2 16.2 16.2 16.2 10.1 16.1
FH 17.7 17.4 13.3 16.5 18.2 9.4 11.0
N2d 17.3 17.2c 12.7 15.5 17.0 17.1 9.9 10.4
N2e 16.7 18.1c 13.9 16.2 18.0 18.5 10.7 11.2
CO 15.1 14.1 11.1 13.7 15.0 15.0 15.0 8.4 9.0
F2 18.2 14.5 13.4 15.9 18.7 9.3 11.1
of the exact Kohn-Sham x-only method. For the case exact-exchange energy term has to be calculated. The
of CO, there is a difference of 0.002 angstroms between same term, however, must also be calculated for the BJ
the AA-BJ and the KS(x-only) result. The BJ(vir) ge- modelpotentialwhenthetotalenergyisobtainednotby
ometriesdifferfromtheKS(x-only)resultsforallfourdi- means of the virial relation, (where it is denoted as E
vir
atomics presented in Table. The BJ(conv) bond lengths in Ref. [23]), but through the exact exchange term (in
coincides with the exact KS(x-only) results for the case whichcase,itisdenotedasE ). The AASCα method
conv
of N , and differ for other molecules. is a simple procedure which permits to transform the
2
BJ model potential into a new AA-BJ potential which
satisfies the Levy-Perdewvirialrelationwithout increas-
IV. CONCLUSIONS ing the computational cost (in fact, this procedure may
be applied to model potentials of any structure without
All potentials discussed here have the same structure: reducing computational efficiency). However, the calcu-
theycompriseaSlaterpotentialplusamodeledresponse lated AA-BJ values show only a slight improvement on
term. Asaresult,thecomputationalcostforallmodelsis the BJ ones, except for the bond lengths, which are in
approximately the same, except that for the AA-PW91, excellent agreement with the KS x-only results.
GLLB, AA-m2, CEDA, LHF cases, where an additional
[1] R.T.SharpandG.K.Horton,Phys.Rev.90,317(1953). [10] S. Hirata, S. Ivanov, I. Grabowski, R.J. Bartlett, K.
[2] J.D. Talman and W.F. Shadwick, Phys. Rev. A 14, 36 Burke, and J.D. Talman, J. Chem. Phys. 115, 1635
(1976). (2001).
[3] J.B. Krieger, Y. Li, and G.J. Iafrate, Phys. Rev. A 45, [11] T.Heaton-Burgess,F.A.Bulat,andW.Yang,Phys.Rev.
101 (1992). Lett. 98, 256401 (2007).
[4] A.Gorling and M. Levy,Phys.Rev.A 50, 196 (1994). [12] A. Gorling, A. Hesselmann, M. Jones, and M. Levy, J.
[5] A.Gorling, Phys.Rev.Lett. 83, 5459 (1999). Chem. Phys.128, 104104 (2008)
[6] S.Ivanov,S.Hirata, and R.J. Bartlett, Phys.Rev.Lett. [13] R. Pino, O. Bokanowski, E.V. Luden˜a, and R. L´opez
83, 5455 (1999). Boada, Theor. Chem. Acc. 123, 189 (2009).
[7] W.YangandQ.Wu,Phys.Rev.Lett.89,143002(2002). [14] O. Gritsenko, R. van Leeuwen, E. van Lenthe, and
[8] W. Yang, P.W. Ayers, and Q. Wu, Phys. Rev. Lett. 92, E.J. Baerends, Phys. Rev.A 51, 1944 (1995).
146404 (2004). [15] F.DellaSala,andA.G¨orling,J.Chem.Phys.115,5718
[9] P. Mori-S´anchez, Q. Wu, and W. Yang, J. Chem. Phys. (2001).
123, 062204 (2005). [16] O.V. Gritsenko and E.J. Baerends, Phys. Rev. A 64,
5
TABLE III:Bond lengths (in angstroms) obtained from different x-onlymethods.
HF KS(x-only) BJ(vir) BJ(conv) AA-BJ
H2 0.734 0.734 – 0.732 0.734
FH 0.897 0.896 0.860 0.895 0.896
N2 1.065 1.065 1.066 1.065 1.065
CO 1.102 1.101 1.086 1.099 1.099
042506 (2001). [27] H. Ou-Yang,and M. Levy,Phys.Lett. 65, 1036 (1990).
[17] M.K.HarbolaandK.D.Sen,Int.J.QuantumChem.90, [28] H. Ou-Yang,and M. Levy,Phys.A 44, 54 (1991).
327 (2002). [29] J. C. Slater, Quantum Theory of Atomic Structure
[18] A.HolasandM.Cinal, Phys.Rev.A72,032504 (2005). (McGraw-Hill, New York 1960), Vol. II.
[19] N.Umezawa, Phys. Rev.A 74, 032505 (2006). [30] V.V. Karasiev, E.V. Luden˜a, Phys. Rev. A 65, 062510
[20] A.D. Becke, and E. R. Johnson, J. Chem. Phys. 124, (2002).
221101 (2006). [31] F.Tran,P.Blaha, andK. Schwartz,J.Phys.: Condens.
[21] V.N.StaroverovandG.E.Scuseria,J.Chem.Phys.125, Matter 19, 196208 (2007).
081104 (2006). [32] E.Engel,andS.H.Vosko,Phys.Rev.A47,2800(1993).
[22] A.F. Izmaylov, V.N. Staroverov, G.E. Scuseria, E.R. [33] F.DellaSala,andA.G¨orling,J.Chem.Phys.118,10439
Davidson,G.Stoltz,andE.Canc´es,J.Chem.Phys.126, (2003).
084107 (2007). [34] J. Kobus, L. Laaksonen, D. Sundholm, Comput. Phys.
[23] V.N.Staroverov,J. Chem. Phys.129, 134103 (2008). Commun. 98, 346 (1996).
[24] V.V. Karasiev, and E.V. Luden˜a, Phys. Rev. A 65, [35] M. Gru¨ning, O.V. Gritsenko, and E.J. Baerends, J.
032515 (2002). Chem. Phys.116, 6435 (2002).
[25] M. Levy and J. P. Perdew, Phys. Rev.A32 (1985) 2010. [36] R. van Leeuwen, and E. J. Baerends, Phys. Rev. A 49,
[26] Q.Zhao,M.LevyandR.G.Parr,Phys.Rev.A47(1993) 2421 (1994).
918.