Density-gradient-corrected embedded atom method

0.95 MB·English

by  Gang Wu

Density Gradient Corrected Embedded Atom Method
Gang Wu,1, Gang Lu,1, Carlos J. Garc´ıa-Cervera,2 and Weinan E3
∗ †
1Department of Physics and Astronomy, California State University Northridge, Northridge, CA 91330-8268
2Mathematics Department, University of California Santa Barbara, CA 93106
3Department of Mathematics, Princeton University, Princeton, NJ 08544-1000
(Dated: January 29, 2009)
Through detailed comparisons betweenEmbeddedAtomMethod (EAM)andfirst-principlescal-
9 culations for Al, we find that EAM tends to fail when there are large electron density gradients
0 present. We attribute the observed failures to the violation of the uniform density approximation
0 (UDA) underlying EAM. To remedy the insufficiency of UDA, we propose a gradient-corrected
2 EAMmodelwhichintroducesgradientcorrectionstotheembeddingfunctionintermsofexchange-
correlation and kinetic energies. Based on the perturbation theory of “quasiatoms” and density
n
functional theory, thenew embedding function captures theessential physics missing in UDA,and
a
paves the way for developing more transferable EAM potentials. With Voter-Chen EAM potential
J
as an example, we show that the gradient corrections can significantly improve the transferability
9
of thepotential.
2

PACSnumbers: 31.15.xv,61.50.Ah,62.20.F-
]
l
l
a
h I. INTRODUCTION ity of EAM is attributable to the uniform background
- density approximationof EAM embedding function. We
s
find that EAM fails whenever there are large gradients
e Atomistic simulations have become an increasingly
m powerfultoolin materialsresearchanda worthypartner of electron density in the system. We overcome the de-
ficiency of the uniform density approximation(UDA) by
. oftheoryandexperiment. Among the greatmanyatom-
at isticmodels,theEmbeddedAtomMethod(EAM)1,2 has proposingadensitygradientcorrectedEAMmodelwhich
m emerged as one of the most successful and versatile ap- incorporates the gradients of the valence electronic den-
sity in the embedding function. Specifically, we intro-
- proaches, representing the mainstay of empirical atom-
d istic simulations. To date, EAM has been applied to duceadditionaltermsintotheembeddingfunctionwhich
n correspond to the density gradient corrections to the
a variety of material systems, such as liquids, metals
o exchange-correlation and kinetic energy contributions.
and alloys, semiconductors, ceramics, polymers, nano-
c MotivatedbythePerdew-Burke-Ernzerhof(PBE)19Gen-
[ structures, and composite materials. Examples of prob-
eralized Gradient Approximation (GGA) of DFT and
lems that EAM has studied include structure, energet-
1 ics and dynamics of lattice defects,3,4,5 elastic response the perturbation theory of “quasiatoms”20, the present
v andphonons,6,7,8 fractureandplastic deformation,9,10,11 model applies to an inhomogeneous background density
4
9 surface and surface growth,12,13,14,15 thermodynamics andhasthecorrectlimitingbehaviorastheexactenergy
5 properties,16 and phase transitions,17 etc. The applica- functions. As a consequence, it extends the applicability
ofEAMandpavesthe wayfor developingmoretransfer-
4 tionsofEAMsimulationshavebeenreviewedinRef. 18.
. ThesuccessandpopularityofEAMareaconsequenceof able potentials.
1
0 its sound theoretical foundation - the density functional
9 theory (DFT) and its simple analytical expression. The
0 former assures that the essential physics be captured by II. FAILURES OF THE UNIFORM DENSITY
v: EAM and the latter endows EAM with excellent numer- APPROXIMATION
i ical efficiency, in par with pair potentials.
X
Despite its great success, EAM suffers from a major First, we demonstrate that the failure of EAM can be
ar deficiency - the lack of transferability. Most of EAM linked to the presence of large gradients of electron den-
models are only reliable in regimes for which they were sity by comparing EAM with first-principles DFT cal-
parameterized; beyond the regimes of parametrization, culations. We establish this fact in bulk Al for which
thereliabilityofEAMpotentialsquicklydeteriorates. As EAM is supposed to work very well. Several excellent
aresult,thepredicabilityofEAMisoftenquestionablein EAM potentials21,22,23 exist and they are used for com-
defect systems and in non-equilibrium conditions where parisons. Both elastic properties and stacking fault en-
relevant physical quantities are not known accurately a ergy of Al are calculated. We compute the cohesive en-
priori and hence not included in the parametrization of ergy per atom and the stress tensor as a function of the
the potentials. As to all empirical models, the lack of rightCauchy-GreendeformationtensorC (i,j =1,2,3)
ij
transferability of EAM is an indication that some theo- for a primitive unit cell of bulk Al. There are six inde-
retical approximations of EAM model are not generally pendent C elements, with the diagonal and off-diagonal
valid. elements varying from -0.28 to 0.28 and -0.18 to 0.18,
In this paper, we show that the lack of transferabil- respectively, giving rise to different elastic deformations.
2
In order to provide precise comparisons, we utilize the
stoparresperegsreindtmaefitnheodh2ig4h--daimnoevneslioanlaglomritehsmhvtehrayteaffillcoiwenstluys. E(eV)0.15 (a) ˚3(eV/A)0.5 (b)
Specifically,wesample483,201datapointsinsparsegrids ∆ 0 σk 0
whichcorrespondtoaregulargridwith65pointsineach −0.∆2Ω00.20.40.6 0.4 0.2 0∆V−0.2−0.4 k∆ −0.∆2Ω00.20.40.6 0.4 0.2∆0V −0.2−0.4
tsooooyfuftmttdahamlfie.traessMtitxpr-yopodrriioneimnoftcsvetienhfprrsel,oeiomsbsnyyass4ntto8edaf3mk,CE2i,ni0Ajwg1Maeatnodfcduva2ra6l4tcn5,hu56teal6ar≈g7terif7oeoo.dn5rfus×w.tche1he0itc1huh0enpwdnoeeuirncmltyasbirnrienygr ∆Ω−0000....02462− 000.7.0.65.40.40.30.2−0.20.1 0 0.10.2(c) 0.4 00001....2468∆Ω−0000....02462− 0.06.0.540.40.30.2−00.1.2 0 0.2(d) 0.4 000...246
The first-principlesDFT calculationsare basedonthe ∆V ∆V
0.6 0 0.6
piamglaepnl(eeV-mwAeaSnvPtee)da2n6i,d2n7P.trhWoejeeVcutieosnernAPauBAgEmb-GeinnGtietAido-wWSitaimhveualmahteiigtohhnopdPl2aa5ncakes-- ∆Ω00..024 −0.4−0.2 0 (e) −0.5∆Ω00..0241110.8 0.40.2 (f) 001..68
weragvyeacnudtosfftreesns.erTgyheokf-p3o6i0ntseVaretosaombptlaeidnarcecloiarbdliengetno- −0.2− −0.60.4 −0.2 ∆0V 0.2 0.4 −1 −0.2− 110.40.6 −0.2 ∆0V 0.2 0.4 00..24
Monkhorst-Pack method with the k-point spacing less 0.6 0.8 0.6 0.8
tpAEhlAsoayMnfeod0rp.0Etoo2At5esM2npt˚Aeieca−dals1l,c.uudpAleatvtGhieoelaonupcss,oesndwiavenebryesgmmeEnpercalcoeroyiloneftgshstoirhefeae0nc.wd1alicedAuVedllayaitsmiuoessnm2e1sd-,. ∆Ω−000...0242 0.30.40.40.2 0.1 0 (g) 0.1 0000...246∆Ω−000...0242 00..650.40.0.420.3 0.1 (h) 000...246
−0.4 −0.2 0 0.2 0.4 −0.4 −0.2 0 0.2 0.4
Mishinet al.22,andVoter andChen23. The firsttwo po- ∆V ∆V
tentialswereconstructedbyfitting to bothexperimental
FIG. 1: The difference for the elastic properties calculated
andfirst-principlesdata,whileVoter-Chenpotentialwas
by EAM and VASP. (a), (c), (e), and (g) show the cohesive
fitted only to experimental data. Although some EAM energy differences in eV/atom; (b), (d), (f), and (h) show
potentials21 were constructed with different motivations the stress differences in eV/˚A3. (a), (b), (c), and (d) are
from the original EAM model, they all use the UDA in calculated by Ercolessi-Adams potential. (c) and (d) are the
effect. For the convenience of presentation, we define corresponding contour plots of (a) and (b), respectively. (e)
∆V = V V0 and ∆Ω= Ω Ω0, where V (V ) and Ω (Ω ) and(f)arecalculatedbyMishinpotentialand(g)and(h)are
V−0 −Ω0 0 0 calculated by Voter-Chen potential. (e), (f), (g), and (h) are
are the volume and solid angle of the deformed (unde-
contourplots.
formed) unit cell. The solid angle is defined relative to
the basis vectors of the unit cell. ∆V and ∆Ω charac-
terize the volumetric and non-volumetric deformation of
the unit cell. of electronic density increases. For a large compression,
the electron density gradient could become too large for
In Figure 1, we present the difference in the cohe-
the UDA of EAM to be valid. Indeed, we find that the
siveenergy(∆E =E −E )andthestresstensor
EAM VASP
density gradient increases considerably in compressions,
(k∆σk=kσ −σ k). Overall, we find excellent
EAM VASP
withthemaximumvalueofthegradientsrisingfrom0.38
agreement between the first-principles and EAM results
˚A 4 for the perfectlattice to 0.54˚A 4 for 40%compres-
over a wide range of deformations - a remarkable feat − −
sion. On the other hand, the expansion of the lattice
of EAM. Furthermore, the errors are insensitive to the
reduces the density gradient, and thus does not violate
solid angle, ∆Ω, which reflects the delocalized nature of
the UDA.
the metallic bonds in Al, and thus justify the use of an
angular - independent model in Al. On the other hand, To further the argument, we perform additional cal-
we find that the EAM errors depend very sensitively on culations for the generalized stacking fault energy (γ)
the change of volume, ∆V; in particular, the EAM val- surface, which along with elastic constants, determines
ues deviate significantly from the first-principles results the plastic behavior of materials. We have carried out
forlargecompressions. Theerrorsinenergycanreachas 293 energy calculations for the entire γ-surface with the
high as 1 eV/atom for 40% compression. This dramatic sparse grid representation. The supercell consists of 9
differencecannotbeaccountedforbythefittingerrorsof layers in the h111i direction for both EAM and VASP
EAM because all three potentials show exactly the same calculations.
behavior. Moreover,the shortestinteratomic distance in The γ-energyalong[12¯1] and[101]directions is shown
thecompressedunitcellis2.1˚A,whichisstillwithinthe in Figure 2a and 2b. Again, overall agreement between
fitting range of the potentials. For example, the fitting the three EAM potentials and VASP is good. However,
rangeofbondlengthisfrom2.0˚Ato6.3˚AinMishinpo- in the neighborhood of the [101] unstable stacking fault
tential and 2.0 ˚A to 5.6 ˚A in Ercolessi-Adams potential. andthe run-onstackingfault (the lasttwo entriesin Ta-
The results suggestthat the errorscome from the model bleI),themagnitudeoftheenergyerrorissignificant. In
itself. particular,thelargesterrorofEAMoccursattherun-on
When the interatomicdistance decreases,the gradient stacking fault in which the atoms in the two neighbor-
3
Ercolessi-Adams (c) III. DENSITY GRADIENT CORRECTION
Mishinetal. 3 MODEL
Voter-Chen
Al
VASP
GCEAM-VC 2
2m) 2 (a) 1 Having established the importance of the density gra-
(J/1.5 [1 dient, we propose a gradient corrected model which
SFenergy0.15 Al −0111]coord cTZoahureledmmbpoado2te0elniotsniabla“lysqeudiamsoipnartotovhmeesp”thi.oenSetteroratinntsgfeetwraoabrl.kilitohyfaSvotefotsEthAoawMndn.
G 00 0.5 [12¯1]1 1.5 Al −2inates that by using a perturbation expansion for an inhomo-
2m)0.8 (b) −3(A)˚ geneous background density, the embedding energy of a
J/0.6 “quasiatom” can be expressed rigorously as a function
(
y −4 of the background density and its gradient. Based on
nerg0.4 Al the “quasiatoms” theory, we introduce three additional
SFe0.2 −5 terms which account for the gradient corrections to the
G
00 0.2 0.4 0.6 0.8 1 −2 −1 0 1 2 exchange,correlationandkineticenergycontributionsto
[101] [11¯2]coordinates(˚A)
the embedding energyof EAM.In this context,the orig-
inal embedding function of EAM can be regardedas the
UDAtotheembeddingenergy. Specifically,thecorrected
FIG.2: γ-energyalong(a)[12¯1]and(b)[101]directionsfrom
embedding function becomes
VASP and EAM calculations. The horizontal axis is in the
unit of the Burgers vector of Al (2.86 ˚A). The contours of
valence electronic density and the density gradient for the
“run-on”configurationareshownin(c). Thearrowrepresents Fi(ρ¯i,si) = F0(ρ¯i)+F˜C(ρ¯i)g(si)
thedirection and magnitude of thedensity gradient. +F˜ (ρ¯)h(s )+F˜ (ρ¯,s ), (1)
X i i G i i
whereρ¯ ≡ ρat(R )isthebackgrounddensityatatom
i j ij
ing (111) layers are right on top of each other, result- j=i
ing in large density gradients. The valence electronic i,andρajt isP6the densitycontributionfromatomj. Rij =
density and its gradient in the “run-on” configuration is R~ −R~ and R~ stands for the atomic coordinates. Al-
i j i
presented in Figure 2c. Noted that the maximum gradi-
ebnletaonfdvatlhenecreuenl-eocntrostnacdkeinnsgityfa,u|glt|smaisxcoofmbpotahratbhlee utonstthae- t(cid:12)(cid:12)(cid:12)tohtoaulgdhetnh(cid:12)(cid:12)(cid:12)seitybaρc(k~rg)rowuhnedredeρn(s~ri)ty=ρ¯i isρnajott t~rhe−sRa~mje a,sththeye
scuorgrgeessptoinngditnhgatvatlhuee foafilluarregseocfomEAprMesscioannsb(e∼in0d.5ee˚Ad−a4t)-, arecloselyrelatedandthegradienPtjsofb(cid:16)o(cid:12)(cid:12)thdens(cid:12)(cid:12)i(cid:17)tiesare
(cid:12) (cid:12)
well defined. In particular, we can define a dimension-
tributed to large density gradients, irrespective of the
specific atomic configurations. Furthermore, the data in lessbackgrounddensitygradientsi assi ∝ ∇R~iρ¯i /ρ¯4i/3,
the brackets of the Table shows a general trend that the (cid:12) (cid:12)
magnitudeofEAMerrorsincreasesasthemaximumden- where ∇R~iρ¯i is the amplitude of backgr(cid:12)(cid:12)ound d(cid:12)(cid:12)ensity
sity gradient increases. Again the shortest interatomic gradien(cid:12)t. In p(cid:12)ractice, si can be approximated by its lo-
distance in Table I is 2.33 ˚A, which is within the fitting calaver(cid:12)(cid:12)age: s(cid:12)(cid:12)≃hs i∝ 1 ∂ρajt(Rij) . F (ρ¯)isthe
range of bond length for the EAM potentials. i i [ρ¯i]4/3j=i ∂Rij 0
UDA embedding function andP6F˜(cid:12)(cid:12)is the gr(cid:12)(cid:12)adientcorrec-
G(cid:12) (cid:12)
tion to the kinetic energy. The leading term of F˜ is of
G
TABLEI: Fault vectorsandenergies forfourstackingfaults the von Weizs¨acker form28, and can be approximated as
obtained from VASPand EAMcalculations. All energies are F˜ (ρ¯,s)=K˜ (ρ¯)k(s). Here K˜ resembles the Thomas-
in mJ/m2. |g|max denotes the maximum gradient of valence FeGrmi kinetic0energy29,30 and k0(s) = λ 1+k11s2+k12s4s2.
electronic density calculated by VASP, and is in ˚A−4. PBE, 01+k21s2+k22s4
λ , k , k , k and k are undetermined parameters.
EA, Mishin and Voter stand for the results calculated by 0 11 12 21 22
For exchange and correlation energy corrections, we
VASP,Ercolessi-Adams, Mishin et al., andVoter-ChenEAM
adopt the functional form of PBE-GGA due to its sim-
potential, respectively. The errors of EAM results are pre-
sentedinbrackets. Rmin (˚A)isthenearestneighbordistance plicity. F˜C andF˜X inEquation1correspondstothecor-
in thecorresponding configurations. relationandexchangeenergyofthelocaldensityapprox-
imation (LDA) of DFT and g(s ) and h(s ) are the cor-
i i
Vector PBE EA Mishin Voter |g|max Rmin responding gradient corrections. The explicit forms F˜
C
1/6[12¯1] 111 121(10) 157(46) 87(-24) 0.398 2.86 andF˜ canbefoundinstandardreferencesofLDA31,32.
X
1/10[12¯1] 184 132(-56) 190(6) 118(-66) 0.399 2.74 Weassumespindegeneracyherealthoughthespinpolar-
1/4[101] 564 663(99) 603(39) 455(-109) 0.507 2.47 izationcanbeconsideredeasilyandcouldbeusefulinthe
1/3[12¯1] 1208 1667(459) 1354(146) 1096(-112) 0.581 2.33 developmentofspin-dependentEAMpotentialsformag-
neticmaterials. Inaddition,werequirethatthemodified
embedding functions have the same limiting behavior as
4
the exact functions: the self-interstitial energies, stacking fault energies, etc.
which involve high density gradient configurations.
g(s)F˜ +h(s)F˜ ∝s2F˜X, s→0 . (2) The GCEAM-VC potential takes the general form of
C X (→−F˜C +κ0F˜X, s→∞ EAM model. The cohesive energy of a system can be
h i written as
We chooseg(s)=− s4 andh(s)= κ0s2 ,which
g0+g1s2+s4 1+αs2
satisfytheaboveconditionsalthoughotherformsofg(s) 1
E = ϕ (R )+ F (ρ¯,s ), (4)
ij ij i i i
and h(s) can also be used. 2
Over all, there are six functions, F˜ , F˜ , K˜ , ρat(R), i,Xj(6=i) Xi
C X 0
F0(ρ), and ϕ(R), which are to be fitted. The first three where embedding function Fi(ρ¯i,si) is expressed in Eq.
are new terms and in conjunction with g(s), h(s) and (1).
k(s), they represent the gradient corrections to the em- TheparameterfittinginGCEAM-VCfollowsthesame
bedding function. procedure of the Voter-Chen potential23, but with two
Since the UDA functions, K˜0(ρ¯), F˜C(ρ¯), and F˜X(ρ¯), modifications. Thefirstmodificationisthatthepairwise
have the same functional forms as their DFT/LDA interaction now is taken the form of:
counterparts29,30,31,32, we have
K˜ (ρ¯)∝ρ¯5/3. (3a) ϕ(R) = ϕ (R)+ϕ (R), (5)
0 1 2
2
ϕ1(R) = DM 1−e−αM(R−RM) −DM,
1 1
F˜C(ρ¯)=ρ¯(c1+c2rs)ln 1+ βrp+1 , rs ∝ ρ¯1/3, ϕ (R) = C RhM n. i
(cid:18) s (cid:19) 2 R
(3b) (cid:18) (cid:19)
Here, ϕ (R) is a Morse potential used in the original
F˜ (ρ¯)∝ρ¯4/3, (3c) 1
X Voter-Chen potential. ϕ (R) is added to account for
2
the repulsive interaction at short distance, ϕ(R→0)→
Replacing the proportional sign ‘∝’ in all the above
∞ when C > 0. However, if one prefers to use fewer
equationswithanequalsign‘=’,onecanfitthe gradient
parameters,ϕ (R)canbeignoredwithoutworseningthe
corrected EAM (GCEAM) potential by introducing 13 2
results (see the discussions of Fig. 4).
additional parameters. These additional parameters are
The second modification is that we do not fit the di-
c , c , β, p, g , g , h , h , λ , k , k , k , and k .
1 2 0 1 0 1 0 11 12 21 22 atomic molecular data. Instead, the force constants of
Amongthem,sevenparameters(k ,k ,k ,k ,g ,g
11 12 21 22 0 1 bulkAlwithdifferentlatticeconstantsarefittedbecause
and α) are introduced through the gradient corrections.
accurateforceconstantsgiverisetoaccuratephonondis-
κ andλ canbeabsorbedintotheembeddingfunctions.
0 0 persions,and hence, accurate thermalproperties such as
thermal capacity and conductivity. The force constants
are fitted for several lattice constants, including 0.9a ,
IV. EXAMPLE: GRADIENT CORRECTED 0
0.95a ,a ,1.05a ,and1.1a ,wherea istheequilibrium
VOTER-CHEN POTENTIAL 0 0 0 0 0
lattice constant of bulk Al. It is found the GCEAM-VC
potentialgivesgooddescriptionfor the diatomicproper-
Inthissection,weapplythegradientcorrectionstothe
tieswithoutfittingthem(SeeTableII).Thisisanexam-
Voter-Chen(VC)potential,andtheresultantpotentialis
ple of improved transferability of the GCEAM model.
termed as GCEAM-VC. It is important to mention that
TheembeddingfunctionofbothVCandGCEAM-VC
the corrections are not constrained in any way by the
potentials is determined by fitting the equation of states
specific form of the EAM potential. We choose the VC
(EOS) to the universal EOS of Rose et al.33. Because
potential because its simplicity - it has only five param-
the universal EOS does not agree exactly with the DFT
eters, much fewer than other EAM potentials, such as
(VASP)values(see Fig. 3),itis inevitablethatbothVC
Mishin and Ercolessi-Adams potentials. As a result, VC
andGCEAM-VCpotentialswoulddeviatefromDFTre-
isnotasaccurateastheotherEAMpotentials. However,
sultsforlargelatticeexpansionsorlargeinteratomicdis-
thesimplicityoftheVCpotentialrendersmoretranspar-
tances. However,asshownlater,thegradientcorrections
entphysics,andfreesusfromintensiveparameterfitting
can improve significantly the description of high density
-whichisnottheemphasisofthepresentpaper. Thegoal
gradient configurations involving lattice defects.
of the article is to illustrate the importance of the den-
The density function of the VC potential is given as
sity gradientcorrections in improving the transferability
of EAM, rather than to generate the best possible EAM ρ(r)=r6 e β2r +29e 2β2r , (6)
− −
potentialfor Al. Had we startedfromanEAMpotential
withmoreparameters,wewouldhavegottenevenbetter β needstobefitted. W(cid:2)ekeepthesamesm(cid:3) oothnesscon-
2
test results for Al. Nevertheless, even with the VC po- ditions for the pairwise interaction, atomic density and
tential,thegradientcorrectionscansignificantlyimprove EOS function in GCEAM-VC as in the VC potential23,
5
0 1
0
−0.2 V)
Ecoh−0.4 R)(e−−21
E(a)/−0.6 ϕ(−3 VGoCtEerA−MC−hVenC
−0.8 Roseetal. −41 1.5 2 2.5 3 3.5 4 4.5 5 5.5
VASP
R(˚A)
−1 0.8 1 1.2 1.4 1.6 1.8 2 2.2 0.8
a/a0 Voter−Chen
0.6 GCEAM−VC
FIG.3: Thecohesiveenergyperatom offcc Alasafunction R)0.4
of thelattice parameter (thescaled equation of states). E ρ(
coh 0.2
isthecohesiveenergyforequilibriumlatticeconstanta0. The
Roseetal.’suniversalEOSisrepresentedbytheblueline,and 0
0 1 2 3 4 5
theVASPvalues are represented byred open circles. R(˚A)
1 Voter−Chen
GCEAM−VC
with the cutoff radius R of these functions to be fit-
cut R)
ted. Thus, there are seven parameters, DM, αM, RM, ρ˜(0.5
C, n, β and R to be fitted before applying the gra-
2 cut
dient corrections. Overall, there are twenty parameters
0
in GCEAM-VC potential, including thirteen parameters 0 1 2 3 4 5
R(˚A)
associated with the gradient correctionterms. The opti-
mized values of all the parameters are presented in Ta-
FIG. 4: The pair interaction function ϕ(R) and atomic den-
ble II.
sity ρ(R)of VC(bluedash) and GCEAM-VC(red solid) po-
tentials.
TABLEII:OnesetofoptimizedparametersforGCEAM-VC.
Length and energy unit are in ˚A and eV,respectively.
Although GCEAM-VC gives a better result for the
c1 -0.45618 g1 -1.60375 k12 -58.09505 RM 1.34426 meltingtemperaturethanVC,thedeviationfromtheex-
c2 -1.03903 h0 0.41370 k21 -10.15397 C 0.43750 perimentalvalueisstilllarge. Thisisduetothefactthat
β 1.17295 h1 0.00327 k22 26.82479 n 4.06038 the melting process is associated with long-range inter-
p 0.29559 λ0 -0.52795 DM 4.00963 β2 3.46742 actions,whereasthe density gradientcorrectionstend to
be short-ranged. Therefore the gradient corrections are
g0 1.63667 k11 157.54365 αM 2.05403 Rcut 5.56250
notexpectedtohavesignificanteffectonmeltingtemper-
ature. This is not an intrinsic problem of the GCEAM
Thepairinteractionϕ(R)andtheatomicdensityρ(R) model because one could improve the melting tempera-
of GCEAM-VC potential are plotted in Fig. 4 in com- turebyfittingmoreaccuratelythelong-rangetailsofthe
parison with the VC potential. It is found that the pair UDA functions, e.g., φ(R), ρat(R), and F (ρ).
0
interaction of GCEAM-VC changes very little from that InFig. 5,wecomparethephonondispersionsbetween
ofVC.AlthoughtheatomicdensityfunctionofGCEAM- GCEAM-VC and VC, against the experimental data. It
VC potential appears to be rather different from that of is found that GCEAM-VC predicts much better results
VC,thisturnsoutnottobethecase. Usingthefactthat than VC.
Eq. 4 is invariant under the transformation:
10
ρ(R)→tρ(R),F(ρ,s)→F(ρ/t,s/t),
8 GCEAM-VC
Voter-Chen
we can define a scaled atomic density ρ˜(R) = z) 6
ρ(R)/max(ρ(R)). ρ˜(R) is plotted in the bottom panel H
T
of Fig. 4 and one finds little difference between the ν ( 4
GCEAM-VC and VC atomic density functions. There-
2
fore,weconcludethatallimprovementstotheVCpoten-
0
tial come from the gradient corrections, i.e., the model G [q00] X K [qq0] G [qqq] L
itself.
Some important properties predicted by VC and FIG. 5: The phonon dispersion curves for Al. Red lines are
GCEAM-VC potentials are collected in Table III. From calculated with GCEAM-VC potential, blue dash lines are
Table III, it is found that the GCEAM-VC improves calculated with VC potential, and open circles are experi-
the overall performance of the VC potential, especially mental data taken from Ref. 40.
for high density gradient configurations, such as self-
interstitials. The table clearly demonstrates the success Furthermore,with the help ofthe sparse-gridmethod,
and improved transferability of GCEAM-VC. we calculate the entire γ-surface using first-principles
6
(a) 0.12
(c)
[101]0000....2468 0.005.05 00.0.010..505025.25000.000000...50.13.2..4133.555450.30.02.1 0.00.505 0.10.35 0.20.000.40.5.000305..50.1.1202.515 0000....1234 [101]0000....2468 0.008.00.60 6 0.00.140.0600..1122 000..1.000406..00082.08 00.06.060.060.012.00.18200..01080.000...160162 00000.....000012468
0 0
0 0.5 1 1.5 0 0.5 1 1.5
[1 2 ¯1] [1 2 ¯1]
(b) 0.14 (d)
[101]0000....2468 0.00.202 0.040.00006..1.00000..8210.068 0.04 0.020.040.006.100..001.22008.002.04 000000......00001124682 [101]0000....24680.0400..006040.0.062 0.00.022 0.040.06 0.02 0.000..6000640.00.4.00002.40.0420.02 0.02 00000.....0000012345
0 0
0 0.5 1 1.5 0 0.5 1 1.5
[1 2 ¯1] [1 2 ¯1]
FIG. 6: The absolute error between the γ-surface from EAM potentials and first-principles calculations ˛˛EGEASFM−EGVASFSP˛˛. In
(a), (b),(c) and (d),the EAM potentials are Ercolessi-Adams, Mishin, VCand GCEAM-VC, respectively.
VASP method. From the γ-surface, one can derive the fact that the modified embedding functions can be fac-
properties of all {111}-type dislocations in Al41. The tored by a ρ-dependent term and an s-dependent term,
absoluteerrorsbetweentheγ-surfacedeterminedbyvar- the analytical expression of force remains simple - it has
ious EAM potentials and VASP calculations are shown several additional terms that are of similar complexity
in Fig. 6. The projection of the γ-surface along two of that of standard EAM. To compute these additional
special orientations is plotted in Figure 2a and 2b. It terms, the GCEAM needs to perform extra calculations
is found that the GCEAM-VC potential yields the most of which the most time-consuming part is the second
accurateresultfor overallγ-surface. Here weshould em- derivativesofthechargedensitywithrespecttodistance.
phasizethat GCEAM-VCdoes notfitany stackingfault By using cubic spline interpolations, these calculations
configurations. In contrast,Ercolessi-Adamsand Mishin canbemaderatherefficientandasaresult,theGCEAM
et al. potentials both have included γ-energies in their force calculation takes less than twice of the CPU time
fitting database, and yet their results are not as good as of the standard EAM. The code package for calculat-
GCEAM-VC. This is an important success of GCEAM ing the energy and force with GCEAM-VC potential is
potentialintermsoftransferability. Moreover,GCEAM- availablevia the WorldWide Web42 orviae-mailatwu-
VCpotentialgivesmuchmoreaccuratestackingfaulten- [email protected].
ergy near the “run-on” configuration ( √3,0 , √3,1 ,
3 3
and 5√3,1 in Fig. 6). These results(cid:16)confir(cid:17)m(cid:16)that in(cid:17)-
6 2 V. DISCUSSION AND CONCLUSION
deed(cid:16)the gra(cid:17)dient corrections are crucial for describing
high density gradient configurations, such as the “run-
Finally, it is instructive to relate the present correc-
on” stacking faults. On the other hand, the errors of tions to other EAM models2,43,44. In the original EAM
GCEAM-VCappearatconfigurationswhereinteratomic
model, the electron correlations arising from the inho-
distance is greater than that of a perfect lattice. These
mogeneous background density are largely ignored. The
errors are not the intrinsic problem of the gradient cor-
goalofthepresentmodelistocapturethemissingcorre-
rection model, but rather due to the fitting strategy of
lations by taking into considerationof density gradients.
the Voter-Chen potential.
Apart from the inhomogeneity of the density, the corre-
Finally, it is useful to mention that the force calcula- lation effect also manifests itself in small molecules and
tion in GCEAM maintains the comparablenumericalef- clusters - a well-known fact in quantum chemistry that
ficiency with the standard EAM models. Thanks to the motivated the development of GGAs. By introducing
7
a PBE GGA-like correction to the exchange-correlation
TABLE III: Properties of Al predicted by VC and and
part of the embedding energy, the present model could
GCEAM-VC potentials in comparison with experimental
improve the description of the correlation effect. The
and/or ab initio data. ∗Fitted properties.
modifiedEAM(MEAM)anditsmultistatevariantstrive
Experimental Voter-Chen toimprovethetransferabilitybymakingthebackground
GCEAM-VC
or ab initio (Ref. 23) density ρ¯ angular and reference-state dependent. How-
ever, since they are based on UDA, the MEAM model
Lattice properties:
does not treat the electron correlations adequately. As
a0(˚A)∗ 4.05a 4.05 4.05
a result, it cannot deal with small clusters accurately
E0(eV/atom)∗ -3.36b -3.36 -3.36 as documented in the literature45. The charge transfer
B(GPa)∗ 79c 79 79.5
EAM (CT-EAM) also recognizes the importance of the
c11(GPa)∗ 114c 107 113 correlation effect. However it addresses the problem by
c12(GPa)∗ 61.9c 65.2 62.5 introducing a reference-state (and its charge) dependent
c44(GPa)∗ 31.6c 32.2 32.9 background density ρ¯. Since the present model consid-
Diatomic Properties: ers exchange-correlationenergy explicitly, it can achieve
De(eV) 1.60d 1.54 1.61 the same goalof CT-EAMwith a simpler function form.
Moreover,onecouldincorporateMEAManditsvariants
Re(˚A) 2.47d 2.45 2.52
intothepresentmodelbymakingthebackgrounddensity
Phonon frequencies:
ρ¯ in Eq. (1) as angular, reference-state and/or charge-
νL(X)(THz) 9.69e 8.55 9.62 dependent if so desired.
νT(X)(THz) 5.80e 5.20 5.36 In conclusion, we have performed detailed EAM and
νL(L)(THz) 9.69e 8.86 10.1 first-principles calculations of Al for elastic deformation
νT(L)(THz) 4.19e 3.70 3.82 and generalized stacking fault energy. We find that al-
νL(K)(THz) 7.59e 6.87 7.69 though EAM models reproduce well the first-principles
νT1(K)(THz) 5.64e 4.80 5.00 results for most cases, they tend to fail when the elec-
νT2(K)(THz) 8.65e 7.76 8.69 tron density gradients become substantial. We attribute
the failures of EAM to the violation of UDA underlying
Vacancy:
the existing EAM models. To remedy the deficiency of
Ef(eV) 0.68f 0.63 0.65
v UDA, we propose an improved EAM model which con-
Self-interstitial:
siders explicitly the gradient corrections to the embed-
Ef(O )(eV) 2.73 2.10 2.41
I h dingfunctionintermsoftheexchange-correlationenergy
EIf(Td)(eV) 3.08 2.55 2.87 and the kinetic energy. We show that the gradient cor-
Ef([111]dumbell)(eV) 2.97 2.48 2.81 rectedmodelcansignificantlyimprovethetransferability
I
Ef([110]dumbell)(eV) 2.76 2.12 2.38 of EAM, and represents a new direction for developing
I
Ef([100]dumbell)(eV) 2.53 2.02 2.30 more transferable EAM potentials.
I
Melting temperature:
Tm(K) 933.6 593.5±10 672.5±10
Acknowledgments
aReference34.
bReference35.
cReference36.
TheresearchatCaliforniaStateUniversityNorthridge
dReference37.
eReference38. was supported in part by DOE grant DE-FC02-
fReference39. 06ER25791andNSFgrantDMR-0611562. CJGC’swork
was funded by an NSF CAREER award. We thank Art
Voter for valuable comments on Voter-Chen potential.
∗ electronic mail: [email protected] 5 A.D. LeClaire, J. Nucl.Mater. 69/70, 70 (1978).
† electronic mail: [email protected] 6 M.S.DawandR.D.Hatcher,SolidStateCommun.56,697
1 M.S. Daw and M.I. Baskes, Phys. Rev. Lett. 50, 1285 (1985).
(1983); Phys. Rev.B 29, 6443 (1984). 7 J.S. Nelson, M.S. Daw and E.C. Sowa, Phys. Rev. B 40,
2 M.I.Baskes, Phys.Rev.Lett. 59,2666 (1987); Phys.Rev. 1465 (1989).
B 46, 2727 (1992). 8 J.S. Nelson, E.C. Sowa and M.S. Daw, Phys. Rev. Lett.
3 S.M. Foiles, M.I. Baskes and M.S. Daw, Phys. Rev.B 33, 61, 1977 (1988).
7983 (1986); ibid.37, 10378 (1988). 9 M.I. Baskes and M.S. Daw, in: Fourth International Con-
4 J.B. Adams, S.M. Foiles and W.G. Wolfer, J. Mater. Res. ference on the Effect of Hydrogen on the Behavior of Ma-
4, 102 (1989). terials, Jackson Lake Lodge, Moran, WY, eds. N. Moody
8
and A. Thompson (The Minerals, Metals, and Materials 27 G. Kresse and J. Furthmu¨ller, Comput. Mater. Sci. 6, 15
Society, Warrendale, PA, 1989). (1996); Phys. Rev.B 54, 11 169 (1996).
10 R.G.Hoagland,M.S.Daw,S.M.FoilesandM.I.Baskes,in: 28 C.F.vonWeizs¨acker,Z.Phys.96,431(1935);D.A.Kirzh-
Atomic Scale Calculations of Structure in Materials, eds. nits, Sov.Phys. JETP 5, 64 (1957).
M.S.DawandM.A.Schliiter(Materials ResearchSociety, 29 L.H.Thomas,Proc.Cambridgephilos.Soc.23,542(1926).
Pittsburgh, PA, 1990). 30 E.Fermi, Z.Phys. 48, 73 (1928).
11 R.G. Hoagland, M.S. Daw and J.P. Hirth, J. Mater. Res. 31 J.P.PerdewandA.Zunger,Phys.Rev.B23,5048(1981).
6, 2565 (1991). 32 J.P.PerdewandY.Wang,Phys.Rev.B45,13244(1992).
12 S.P. Chen, A. Voter and D.L. Srolovitz, Phys. Rev. Lett. 33 J.H. Rose, J.R. Smith, F. Guinea, and J. Ferrante, Phys.
57, 1308 (1986). Rev.B 29, 2963 (1984).
13 S.P. Chen, D.J. Srolovitz and A.F. Voter, J. Mater. Res. 34 C. Kittel, Introduction to Solid State Physics (Wiley-
4, 62 (1989). Interscience, New York,1986).
14 T. Ning, Q.Yu and Y.Ye, Surface Sci. 206, L857 (1988). 35 HandbookofChemistryandPhysics,editedbyR.C.Weast
15 K. Takahashi, C. Nara, T. Yamagishi, and T. Onzawa, (CRC, Boca Raton, FL, 1984).
Appl.Surf. Science 151, 299 (1999). 36 G. Simonsand H.Wang, Single Crystal Elastic Constants
16 S.M.FoilesandM.S.Daw,Phys.Rev.B38,12643(1988). and Calculated Aggregate Properties (MIT Press, Cam-
17 S.M.FoilesandJ.B.Adams,Phys.Rev.B40,5909(1989). bridge, MA, 1977).
18 M.S.Daw,S.M.Foiles,andM.I.Baskes,Mater.Sci.Rep. 37 K.P. Huber and G. Hertzberg, Constants of Diatomic
9, 251 (1993). Molecules (Van Nostrand Reinhold,New York,1979).
19 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. 38 R.Stedman and G. Nilsson, Phys.Rev. 145, 492 (1966).
Lett. 77, 3865 (1996). 39 H.-E.Schaefer,R.Gugelmeier,M.Schmolz,andA.Seeger,
20 M.J.StottandE.Zaremba,Phys.Rev.B22,1564(1980). Mater. Sci. Forum 15-18, 111 (1987).
21 F.ErcolessiandJ.B.Adams,Europhys.Lett.,26(8),583 40 R.Stedman and G. Nilsson, Phys.Rev. 145, 492 (1966).
(1994). 41 G. Lu, N. Kioussis, V.V. Bulatov, and E. Kaxiras, Phys.
22 Y.Mishin, D.Farkas,M.J.MehlandD.A.Papaconstan- Rev.B 62, 3099 (2000).
topoulos, Phys. Rev.B 59, 3393 (1999). 42 http://wugaxp.com/Documents/GCEAM-VC.tar.gz.
23 A.VoterandS.Chen,Mat.Res.Soc.Symp.Proc.82,175 43 M.I. Baskes, S. G. Srinivasan, S. M. Valone, and R.G.
(1987). Hoagland, Phys.Rev.B 75, 094113 (2007).
24 H.-J. Bungartz and M. Griebel, Acta Numerica 13, 147 44 S.M.ValoneandS.R.Atlas,Philos.Mag.86,2683(2006).
(2004). 45 M.I. Baskes, J.S. Nelson and A.F. Wright, Phys. Rev. B
25 G. Kresse and D. Joubert, Phys.Rev. B 59, 1758 (1999). 40, 6085 (1989).
26 G. Kresse and J. Hafner, Phys.Rev.B 48, 13115 (1993).