An Efficient Algorithm for Density Functional Theory Simulation of Large Quantum Dot Systems

7 Pages·0.17 MB·English

by  Hong Jiang

AnEfficient Algorithm forDensity Functional Theory SimulationofLargeQuantum DotSystems
Hong Jiang,1,2,3 Harold U. Baranger,2,∗ and Weitao Yang1,†
1Department of Chemistry, Duke University, Durham, North Carolina 27708-0354
2Department of Physics, Duke University, Durham, North Carolina 27708-0305
3College of Chemistryand Molecular Engineering, Peking University, Beijing, China 100871
3
(Dated:February2,2008)
0
0 Kohn-Shamspin-densityfunctionaltheoryprovidesanefficientandaccuratemodeltostudyelectron-electron
2 interactioneffectsinquantumdots,butitsapplicationtolargesystemsisachallenge.Anefficientalgorithmfor
n thedensity-functionaltheorysimulationofquantumdotsisdeveloped, whichincludestheparticle-in-the-box
a representationoftheKohn-Shamorbitals,anefficientconjugategradientmethodtodirectlyminimizethetotal
J energy,aFourierconvolutionapproachforthecalculationoftheHartreepotential,andasimplifiedmulti-grid
3 techniquetoacceleratetheconvergence. Thenew algorithmistestedina2Dmodelsystem. Usingthisnew
1 algorithm, numerical studies of large quantum dots with several hundred electrons become computationally
affordable.

]
l PACSnumbers:
l
a
h
- I. INTRODUCTION form(FFT)methodcanbeusedtotakefulladvantageofthe
s
periodicityofcrystalstructures. Inprinciple,thePWmethod
e
m is validonlyfor periodicsystems, butaperiodicsystems can
Quantumdots(QD)are onekindof nano-deviceinwhich
betreatedbyintroducingthesuper-celltechnique.21Inrecent
. the motion of electrons is quantized in all three dimen-
t years, severalgroupshave been advocatingthe use of basis-
a sions through the lateral confinement of a high-mobility
m modulation-doped two-dimensional electron gas in a semi- freerealspacemethodsforelectronicstructurecalculationsof
finitesystems,inwhichthewavefunctionsarerepresentedin
- conductor hetero-structure.1,2,3 Both quantum interference
d realspace,andthekineticenergyoperatorisdiscretizedbya
and electron-electron interactions play important roles, and
n high-orderfinitedifference(FD)method.23Withagivenrep-
their interplay underlies many properties that are fascinat-
o resentation, there are various methods to solve the resultant
ingintermsofbothfundamentalmesoscopicphysicsandfu-
c
numericalKS equation. Roughlythey fallinto two different
[ ture technical applications. Various theoretical models have
types: methods that minimize the total energy directly, and
been developed to explain experimental discoveries and to
1 those thatsolve the KS equationin a self-consistentway. In
predict new properties. The Kohn-Sham (KS) density func-
v addition, for DFT simulations of finite systems, another im-
tionaltheory(DFT)method4,5,6 providesanaccuratenumeri-
6
portantissueisthecalculationoftheHartreepotential.
7 calmodeltostudyelectron-electroninteractioneffectsinQD
1 systems.7,8,9,10,11,12,13,14,15,16,17,18 In spite of its comparatively Aiming at modelling QD systems efficiently, we have de-
1 lowcomputationalcost,previousDFTcalculationsarelimited velopednewtechniquesinboththenumericalrepresentation
0 tosystemsinwhichtheelectronnumberisgenerallylessthan ofKS orbitalsand the solutionofthe KS equation. We note
3 a few tens. In many experimental cases, however, the elec- thatin QD systems, the wave functionsvanishat the bound-
0
tron numbers involved are more than several hundred. The aries, and in some cases eventhe hard-wallboundarycondi-
/
t main theoretical approaches in the regime of large number tion is used. For a function in a rectangular box with zero
a
m of electrons are statistical methods2,3,19,20 which are usually boundaryvalues,themostnaturalbasissetistheparticle-in-
based on some generalassumptionswhose validity, in many the-box (PiB) basis set. The kinetic energy operator is di-
-
d cases,isyettobejustified. Itisthereforedesirabletoexplore agonalin the PiB space and the transformationbetween real
n the large dot regime directly by using numerically more ac- and PiB space can be efficiently performedby the fast sine-
o curatemodelssuchasDFT.Thisimposesademandingcom- transform (FST) method, which is a variant of the FFT. For
c putationaltaskbecausetoobtainmeaningfulstatistics, many thesolutionoftheKSequation,wemodifiedTeter,Payneand
:
v calculationswith several hundredelectrons need to be done. Allan (TPA)’s band-by-bandconjugate gradientmethod24 to
i It is therefore compelling to develop more efficient numeri- getamoreefficientdirectminimizationapproach.
X
caltechniquesforDFTsimulationofQDsystems;thisisthe The outline of the paper is as follows. In the next sec-
r
a mainobjectofthecurrentstudy. tion, after a simple description of the KS-DFT method, we
Two key issues are involved in the numerical implemen- present the main components of our algorithm: (1) the PiB
tation of the KS-DFT method:21,22 (1) the numerical repre- representation of the KS equation; (2) a modified band-by-
sentationofwavefunctionsandtheKSHamiltonian,and(2) band conjugate gradient method for the direct minimization
the solution of the numerical KS equation. While local ba- ofthe KS totalenergy;(3)a Fourierconvolutionmethodfor
sis sets (mainlyGaussian-typeorbitals) dominate in conven- thecalculationoftheHartreepotentialthatwasdevelopedby
tional quantum chemistry of molecular systems, the plane MartynaandTuckerman25inthePWpseudo-potentialcalcu-
wave (PW) basis,21,22 combined with the pseudo-potential lations;and(4)asimplifiedmulti-gridtechniquetoaccelerate
method,iswidelyusedintheabinitioelectroniccalculations theconvergence.InSectionIII,thenewalgorithmistestedin
of various material systems, in which the fast-Fourier trans- a2DmodelQDsystemwithelectronnumberN =100. Sec-
2
tionIVsummarizesthemainresultsandconcludesthepaper. approximation(GGA)areavailable,ithasbeenshownthatthe
GGAresultsareclosetothosefromLSDAcalculationsinQD
systems.13 Intermsoftheimplementation,thecalculationof
II. METHOD V istrivialwhenthespindensitiesareinrealspaceasisthe
xc
caseinouralgorithm,butthecalculationofV requiresmore
H
A. Kohn-Shamspin-density-functionaltheory effortsaswillbeshownlater.
Considering the important role played by electron spin in
QD systems, we take the effect of spin polarization explic- B. PiBRepresentation
itlyintoaccountintheframeworkofKohn-Shamspin-density
functional theory (KS-SDFT).4,6 In KS-SDFT, the ground To simplify the notation, we will take 1-D systems as an
state energy of an interacting system with electron number example,thegeneralizationtohigherdimensionalcasesbeing
N andthetotalspinS in thelocalexternalpotentialVext(r) straightforward.Anyregularfunctionf(x)thatislocalizedin
iswritten asa functionalof spindensitiesnσ with σ = α,β thefiniteregion0<x<Lwithzeroboundaryvaluescanbe
denotingspin-upandspin-down,respectively, expandedas
E nα,nβ =T nα,nβ + n(r)V (r)dr 2 nπx
(cid:2) (cid:3) s(cid:2) (cid:3) Z ext f(x)= CnrLsin L (7)
1 n(r)n(r′) Xn
+ drdr′+E nα,nβ . (1)
2Z |r−r′| xc(cid:2) (cid:3) andtheexpansioncoefficientsCnare
(Effectiveatomicunitsareusedthroughthepaper:forGaAs-
2 L nπx
AlGaAsQDs,thevaluesare10.08meVforenergyand10.95 C = f(x)sin dx. (8)
nm for length.) Ts nα,nβ is the kinetic energy of the KS n rLZ0 L
non-interactingrefe(cid:2)rencesy(cid:3)stem whichhasthe sameground
state spin densityas the interactingone, andE nα,nβ is Integrating Eq. (8) numerically on a set of equally spaced
xc discretepoints x j∆x j L usingtheextendedtrape-
theexchange-correlationenergyfunctionalThesp(cid:2)indensi(cid:3)ties { j ≡ ≡ Nx}
nσ satisfytheconstraint nσ(r)dr = Nσ withNα = (N + zoidalformula26leadsto
2S)/2andNβ =(N 2RS)/2.
Assumingthattheg−roundstateofthenon-interactingrefer- 2 Nx−1 πjn √2L
C = ∆x f sin =F (9)
encesystemisnon-degenerate,thenon-interactingkineticen- n rL j N n N
ergyis givenby T nα,nβ = ψσ 1 2 ψσ , and Xj=1 x x
s i,σ i −2∇ i
thegroundstatespin(cid:2)density(cid:3)isunPiquel(cid:10)yexp(cid:12)ressed(cid:12)as (cid:11) withf f(x )and
(cid:12) (cid:12) j j
≡
Nσ
nσ(r)= ψσ(r)2, σ =α,β. (2) Nx−1 πjn
Xi | i | Fn ≡ Xj=1 fjsin Nx ≡FST{fj} (10)
Hereψσ arethelowestsingle-particleorbitalswhichareob-
tainedfirom where FST fj denotes the fast sine-transform of the data
{ }
f .
j
Hσ ψσ(r)=εσψσ(r), (3) { O}neofthekeyingredientsinouralgorithmistheactionof
KS i i i
thesingle-particleHamiltonianoperatoronwavefunctions,
withtheKSHamiltonianH definedas
KS
Hf(x)=Tf(x)+Vf(x), (11)
1
Hσ 2+V (r)+V [n;r]+Vσ[nα,nβ;r]. (4)
KS ≡−2∇ ext H xc
the efficiency of which is critical for the performanceof the
V [n;r] and Vσ[nα,nβ;r] are the Hartree and exchange- wholealgorithm.Whereasthepotentialenergyoperatorisdi-
H xc
correlationpotentials,respectively, agonalinrealspace,thekineticenergyoperatorisdiagonalin
PiB space. Wave functions can be transformed between the
n(r′) twospacesefficientlybythefastsinetransform. Inouralgo-
V [n;r] d3r′, (5)
H ≡ Z r r′ rithm,wavefunctionsareindiscreterealspace fj ,andthe
| − | application of the potential energyoperatoris t{here}foretriv-
δE nα,nβ
Vσ[nα,nβ;r] xc . (6) ial, V f = V f , whereV is the value ofthe potential
xc ≡ δn(cid:2)σ(r) (cid:3) { j} { j j} j
atthepointx . Thekineticenergyoperatorisappliedtowave
j
functionsinPiBspace,
We have used the local spin-density approximation
(LSDA)4,6 for E , which is widely used for the modelling
xc ¯h2 2 ¯h2k2 nπx
of material systems. Although more accurate exchange- Tf(x) 2f(x)= F n sin , (12)
n
correlationfunctionalformssuchasthegeneralized-gradient ≡−2m∇ Nx Xn 2m L
3
withk =nπ/L,whichinthediscreteformbecomes where
n
T{fj}= N2 FST{Fn¯h22mkn2}. (13) γim = ζmhζ−im1 |ζζimm−i1 (16)
x i i
(cid:10) (cid:12) (cid:11)
ThePiBrepresentationformulatedaboveiscloselyrelated withγ1 = 0. Theconjugatevecto(cid:12)ris furtherorthogonalized
to the PW method. In fact, for finite systems with soft-wall tothepiresentband ψm andnormalized(Nisdenotedasthe
boundaries, the two representations are numerically equiva- | i i
normalizationoperator),
lent. Mathematically, however, they are different: The PiB
basissetisrealandthezeroboundaryconditionisimposedby ϕ′m =N(1 ψm ψm ) ϕm (17)
thebasissetitself,butinthePWcasethebasisfunctionsare | i i −| i ih i | | i i
complex and the wave functionsare forced to be zero at the Thenewwavefunctionforthei-thorbital ψm+1 isformed
i
boundariesby the external potential of systems under study. fromthelinearcombination (cid:12) (cid:11)
(cid:12)
ThePWmethodwillfailinthecaseofthehard-wallbound-
aryproblem,whichisquitecommoninthestudiesofquantum ψm+1 = ψm cosθ + ϕ′m sinθ (18)
i | i i min | i i min
dots,butthePiBmethodisstillvalid. (cid:12) (cid:11)
which is g(cid:12)uaranteedto remain normalizedand orthogonalto
all other orbitals. θ is obtained by minimizing the to-
min
C. Direct-minimizKaotihonn-cSohnajmugeaqtueagtriaodnientmethodforthe tϕal′menesrignyθ.aIsnaTfPuAn’cstaiolgnoorifthθmw,2i1t,2h4|θψi(θi)sid=eter|mψiminiedcobsyθth+e
| i i min
followingapproximatescheme:Thetotalenergyasafunction
In direct minimization approaches, the KS total energy ofθisapproximatedbyE(θ) Eavg+A1cos2θ+B1sin2θ;
≈
functional is minimized directly over orbital wave functions thethreeunknowns,Eavg,A1andB1,aredeterminedaccord-
underorthonormalconstraints. Wehavemadeseveralimpor- ingtothreepiecesofinformation: E(θ = 0), ∂E(θ) and
∂θ |θ=0
tantmodificationstoTPA’sband-by-bandconjugate-gradient E(θ =π/300).
scheme to obtain higher efficiency. Here we give an outline Herewe proposea moreefficientapproximateschemefor
of the algorithm, and emphasize the modifications we have thedeterminationofθ .ThederivativeofE(θ)withrespect
min
made. toθcanbeobtainedfrom
The basic idea of a band-by-band scheme is to minimize
∂E(θ)
the total energy over one band (or orbital) at a time, which, =2 ϕ′m H (θ) ψm cos2θ
compared to other conjugate-gradient schemes,27,28 has the ∂θ h i | KS | i i
following advantages: (1) much lower requirement for stor- −(hψim|HKS(θ)|ψimi−hϕ′im|HKS(θ)|ϕ′imi)sin2θ
age space, (2) simpler implementation, and (3) particularly (19)
in our algorithm, an efficient approximateline-minimization
scheme,aswillbeshown. AssumingH (θ) H (0),from ∂E(θ) =0weget
KS ≈ KS ∂θ
Firstthesteepestdescent(SD)vectorforthei-thorbitalat
them-thiterationiscalculatedfrom 1 B
θ = tan−1 (20)
min
2 A
|ζimi=1− |ψjihψj|(λmi −HKS)|ψimi (14) with
Xj6=i
  A= ψm H (0) ψm + ϕ′m H (0) ϕ′m (21)
withλm = ψm H ψm (tosimplifythenotation,weuse −h i | KS | i i h i | KS | i i
i h i | KS| i i
Dirac’s state vector notation and drop the spin index in the and
followingformulation). InTPA’salgorithm,theSDvectoris
preconditioned before it is used to build the conjugate vec- B =2hϕ′im|HKS(0)|ψimi. (22)
tor.Inourcalculations,however,itwasfoundthatalthoughin
Theunderlyingapproximationinourschemeissimilartothat
manycasespreconditioningdoesacceleratetheconvergence,
of TPA’s, but our scheme is much more efficient in terms of
itseffectisnotalwayspositive.Ontheotherhand,thecompu-
computationaleffort: InTPA’sscheme,ateachbanditeration
tationaloverheadinthepreconditioningstep,whichinvolves
thetotalenergymustbecalculatedtwice, i.e. E(θ = 0)and
anotherorthogonalizationprocessaswellastheactionofthe
E(θ =π/300);inourscheme,themosttime-consumingstep
preconditioningoperatoron the SD vector,can be expensive
istheactionofH on ϕ′m ,whichismuchfasterthanthe
forlargesystems. Asshownlater,byusingasimplifiedmulti- KS | i i
calculation of the total energy. Considering further the fact
gridtechnique,wecanachievefastconvergenceevenwithout
thatthetotalnumberofbanditerationscanbeverylarge,we
preconditioningoftheSDvectors.
The conjugate vector ϕm is then constructed as a linear note that an efficient line-minimizationscheme such as ours
combinationoftheSDve|ctoiriζm andthepreviousconjugate iscrucialtoreducethecomputationaleffort.
| i i InTPA’s algorithm,after thewave-functionis updatedac-
vector,
cordingto Eq. (18), the Kohn-ShamHamiltonianis updated
ϕm = ζm +γm ϕm−1 (15) immediately,whichinvolvesthereconstructionofV (r)and
| i i | i i i i H
(cid:12) (cid:11)
(cid:12)
4
V (r)accordingtothenewdensity.Thisisactuallyquiteex-
xc
pensiveforlargesystems. Ontheotherhand,weexpectthat
theKSHamiltonianwillnotexperiencelargechangesinside
theiterationsofasingleorbital. Soinouralgorithm,weup-
dateH aftereveryN banditerations,andtheoptimal
KS update
valueofN willbeexploredinthenextsection.
update
In each orbital, the proceduredescribed above is repeated
N times; the iterations are then started on the next or-
band
bital. After the wave functionsof all orbitals are updated in
thisway,thetotalenergyiscalculatedandiscomparedtothat
ofthepreviouscycletodetermineif thefinalconvergenceis
achieved. The main parameters in the algorithm are N
band
and N . Their effects on the performance of the algo-
update
rithmwillbetestedindetailinthenextsection.
FIG.1: Illustrationof theFourierconvolution methodforthecal-
culationoftheHartreepotentialinfinitesystems. Theoriginalgrid
wherethedensityandpotentialaredefinedisextendedintodoubled
D. CalculationofVH space. Thedensityintheextended regionissettozero. Imposing
theperiodicboundaryconditionforthisextendedgrid,theunphysi-
calinteractionbetweenneighboringsuper-cellscanbeavoided.
For finite systems, the simplest and perhaps most ineffi-
cientway to calculate V is by directnumericalintegration,
H
which is feasible only for small systems. Another widely- as
used approachis to solve the Poisson equationequivalentto
Eq. (5). Though the Poisson equation itself can be solved n(r) ifr Ω
with great efficiency, the calculation of boundaryvalues can n2L(r)=(cid:26) 0 othe∈rwise (24)
be quite expensive even by using efficient multipole expan-
siontechniques. Additionally,thePoissonsolverapproachis Imposing periodic boundary conditions to both the density
valid only for 3D systems; in the case of 2D systems, there andthe Coulombinteractionkernelin the extendedgrid, the
isnoPoissonequationequivalenttoEq. (5). Inrecentyears, potentialcanbecalculatedaccordingtotheconvolutiontheo-
several schemes have been proposed to extend the conven- rem,
tional Fourier convolution method to finite systems.25,29,30,31
InparticularwehaveincorporatedMartynaandTuckerman’s V (r)= n (k)v (k)eik·r (25)
H 2L c
method25 into our algorithm. Considering that the Martyna- Xk
Tuckermanmethod was developedmainly for the modelling
where n (k) and v (k) are respectivelyfinite Fourier inte-
of molecular and material systems within the plane-wave 2L c
gralsofthedensityandtheCoulombinteractionkernelinthe
pseudo-potentialframework,we willformulatethe approach
extendedgrid,
herewithsomedetail.
The calculation of the Hartree potential is straightforward
1
for periodic systems, but this is not the case for finite aperi- n (k) = n (r)e−ik·rdr, (26)
2L 2L
odicsystems. ThepotentialV (r)hastheformoftheconvo- Ω2L ZΩ2L
H
lution between the density and the Coulombinteraction ker- v (k) = v (r)e−ik·rdr (27)
nel, v (r) = 1/r, whichhasthefollowingsimplerelationin c Z c
c Ω2L
theFourierspace,26
whereΩ isusedtodenoteboththeextendgridanditsvol-
2L
V˜ (k)=n˜(k)v˜ (k), (23) ume(orareainthe2Dcase).
H c While n (k) can be easily obtained from the discrete
2L
wheref˜(k)referstotheFouriertransformoff(r). Eq.(23)is FouriertransformofitsrealspacevaluesbyFFT,thecalcula-
usefulonlywhenwehavetheanalyticalformofn˜(k)andcan tionofvc(k)ismuchmoreinvolvedbecauseofthesingular-
ityoftheCoulombinteractionkernelinrealspace.Thekeyto
performthe inverseFouriertransformofV˜ (k) analytically,
H Martyna-Tuckerman’sapproachistodecomposetheCoulomb
whichisnottrueformostcases wherethedensityisusually
interactionkernelintolong-andshort-rangeparts,
representedindiscreterealspace.
When applying the Fourier method to discrete finite sys- erf(αr) erfc(αr)
tems, periodic boundary conditions are always assumed. It vc(r)= r + r ≡vc(long)(r)+vc(short)(r) (28)
haslongbeenknownthattheunphysicalinteractionsbetween
neighboringsuper-cellscanbeavoidedbycalculatingV ina whereerf(x) and erfc(x) arethe errorfunctionandits com-
H
doublyextendedgrid.32 Inparticularfor2Dsystemsasillus- plement,respectively,andαistheparameterthatcontrolsthe
tratedinFig. 1,theoriginalL L grid(Ω)isextendedto effectivecut-offrange. ThefiniteFourierintegraloftheshort
x y
×
2L 2L (Ω ).Thedensityintheextendedgridisdefined range part can be well approximated by its infinite Fourier
x y 2L
×
5
transform, 100
HOFD, 5-Point
y
v(cshort)(k) ≡ Z vc(short)(r)e−ik·rdr Energ10-1 HPiOBFD, 13-Point
Ω2L v(short)(r)e−ik·rdr v˜(short)(k). Total 10-2
≈ Zwholespace c ≡ c ed 10-3
g
(29) ver
on10-4
C
which is analytically known in both 2D and 3D cases. The n
finiteFourierintegralofthelong-rangeinteractioncanbedi- or i10-5
rectlyobtainedfromthediscreteFouriertransformofitsreal Err
10-6
space values. In the practical implementation, v needs to
c 32 48 64 80
be calculated onlyonce atthe beginning. The calculationof Number of Grid Points in Each Direction
V (r) involves only two FFTs (one forward and one back-
H
ward), which makes this approach much more efficient than FIG.2: Convergenceofthetotalenergywithrespecttothenumber
methodsbasedonaPoissonsolver. ofgridpointsusing5-PointFD,13-PointFDandPiBrepresentation.
E. Asimplifiedmulti-gridtechnique pointsisNx = Ny = 64. Allthefollowingnumericalresults
come from calculationswith electronnumberN = 100 and
spin S = 0. For the exchange-correlationenergy, E , we
The multi-grid method is an efficient technique to ac- xc
useTanatarandCeperley’sparameterizedformoftheLSDA
celerate the convergence in various real space relaxation
functional.34 The convergence criterion is set as ǫ = 10−6,
approaches.26Thebasicideaisthatlow-frequencyerrorsare
whichcorrespondstoabout10−5meVinGaAs-AlGaASQD
easiertoeliminateinacoarsegridthaninafinegrid. Leeet
systems.
al.33proposedasimpleone-waymulti-gridtechniqueinfinite-
Considering that the FD method has been widely used in
differencereal space KS calculations. Wave functionsbeing
thenumericalmodellingofQDsystems,7,8,13,16wefirstmake
representedinreal-spaceinouralgorithm,asimilarsimplified
acomparisonbetweenFDandPiB.IntheFDrepresentation,
multi-grid (SMG) method can be implementedin a straight-
the second order derivative in the kinetic energy operator is
forwardway: theKSenergyfunctionalisfirstminimizedina
locallydiscretizedinrealspace
coarsegrid;theconvergedwavefunctions,afterinterpolation
andre-orthogonalization,aretakenastheinitialguessforthe
∂2 1 m
minimizationonthefinegrid.Withthesewell-preconditioned f(x)= C f +O(h2m) (31)
initialwavefunctions,theconvergenceinthefinegridcanbe ∂x2 h2 j=X−m j j
easilyattained.
wherehisthediscretizationstepandthecoefficientsC can
j
be obtained systematically for any m. Fig. 2 illustrates the
III. NUMERICALTESTS convergence of the total energy with respect to the number
of grid points using 5-point (m = 2), 13-point (m = 6)
InmostexperimentalQDsystems,theexcitationinthever- FDandPiBrepresentations.Thelower-orderfinitedifference
tical direction can be neglected, and a 2D model is a good scheme,m = 2, ispoorin termsofaccuracy,andconverges
approximation.Inthispaperwewillreporttestresultsonlyin slowlyasthegridsizeincreases. Thehigh-orderFDscheme,
a2Dsystem;theextensionto3Dsystemsisstraightforward, m = 6, doesimprovethe accuracyby two ordersof magni-
andthe highefficiencyofouralgorithmmakesitpossible to tude, butis still less accuratethan PiB for Nx = 32, 48and
explore large dots even in the 3D case. We test the perfor- 64.
mance of the new algorithm in a coupled quartic oscillator We check the accuracy and efficiency of our line-
potentialsystem, minimization scheme by comparing with TPA’s approach as
wellasthenumericallyexactBrent’slinesearchalgorithm.26
x4 In thiscase, we use N = 5 and N = 1 as recom-
V (r)=a +by4 2λx2y2+γ(x2y xy2)r , band update
ext (cid:18) b − − (cid:19) mended in TPA’s original algorithm.21 Fig. 3 plots θmin in
(30) thefirst25banditerationscalculatedfromthreeschemesre-
witha = 10−4,b = π/4,λ = 0.6andγ = 0.1. Considering spectively.Thevaluesofθ calculatedfrombothTPA’sand
min
thattheconvergencebehaviorinKScalculationsisgenerally ourmethodagreeverywellwithexactvaluesandtherelative
related to the effective interaction strength, which is charac- errorsarealwayssmallerthan1%. Butintermsofcomputa-
terizedbyr (πn)−1/2/a in the2D case wherenis the tionaleffort,ournewschemeismuchmoreefficientasargued
s B
≡
average density and a is the Bohr radius, we have chosen intheprevioussection.
B
the parameters in V so that the estimated r is about 1.5, To find the optimal N and N , we do the calcu-
ext s band update
which is close to experimental values. The calculations are lations with different values of N and N , and the
band update
doneinagridofsizeL = L = 50andthenumberofgrid results are shown in Fig. 4. With fixed N = 1, it is
x y update
6
0.5 tationaleffortinoneKScalculationasafunctionoftheorder
Ours oftheLagrangeinterpolationformulainbothtwo-andthree-
0.4 a TPA level SMG calculations. We see that a high-orderinterpola-
Exact
min tionschemeisusefultoimprovetheperformance.
θ
1
0.3
Two-Level
Three-Level
0.2
e
m0.8
%) 0.4 b Ti
or ( 0.2 PU
e Err e C0.6
v 0 v
Relati-0.2 Relati
0 5 10 15 20 25 0.4
FIG.3: Comparisonofθ calculatedbythreedifferentschemes.
min
0 2 4 6 8
Order of Lagrange Interpolation
(a) (b)
N =5 N =1
102 band update FIG. 5: Relative computational efforts for one KS calculation as
y Nband=10 Nupdate=5 afunctionoftheorderofinterpolationinbothtwo-andthree-level
al Energ 100 NNNbbaanndd===234000 NNuuppddaattee==1200 StaMkeGnacsaltchuelautnioitn.s. TheCPUtimeinthesingle-level calculationis
ot band
n T10-2
r i IV. SUMMARY
o
r
Er10-4
In this paper, we presented an efficient algorithm for the
KS-SDFTsimulationoflargequantumdotsystems. Themain
10-6
elementsof the algorithmare: (1) Wave functionsare repre-
0 5000 10000 0 1000 2000 3000
CPU Time(s) CPU Time(s) sentedinrealspace,andthekineticenergyoperatorisapplied
towavefunctionsbyfastsinetransform. (2)TheHartreepo-
FIG.4: ErrorsinthetotalenergyasafunctionofCPUtimeduring tentialis calculatedbyMartyna-Tuckerman’sFourierconvo-
the DMCG calculation for different Nband with Nupdate = 1 (a) lutionmethod.(3)ForthesolutionoftheKSequation,wein-
andfordifferentNupdate withNband = 20(b)withrespecttothe troducedseveralimportantmodificationstoTeteretal.’sband-
exacttotalenergythatiscalculatedusingafinergrid(Nx =80)and by-bandconjugate-gradientmethod.Amoreefficientapprox-
tighterconvergencecriterion(ǫ=10−7),
imateline-minimizationschemewasdeveloped;itwasfound
thatlargebanditerationnumberanda delayedupdateof the
KS Hamiltonian inside the band iterations increase the effi-
seenthatarelativelylargerN ismoreefficientthansmall
band ciencybyoneorderofmagnitude.(4)Asimplifiedmulti-grid
N . FixingN =20,N =20givesthebestper-
band band update techniquewasintroducedtoacceleratetheconvergence. The
formance. The combinationof large N and N re-
band update new algorithm has been used to study spin and conductance
ducesthecomputationaleffortbyalmostoneorderofmagni-
peak-spacing distributions in a 2D chaotic QD system with
tude.ThoughtheactualvaluesofoptimalN andN
band update electronnumberN upto200,fromwhichnewphysicalphe-
mayvaryfordifferentsystems,thebasicideademonstratedin
nomenawererevealed.18
thistestcalculationisbelievedtobeofgeneralsignificance.
We have implemented both the two-level (h and 2h) and
three-level(h,4h and 2h)SMG schemes. Insteadof usinga
3
sophisticated interpolation as in Ref. 33, we use a simpler
Acknowledgments
Lagrangepolynomialinterpolationmethod.Comparisonwith
the single-level calculation shows a quite obvious improve-
ment in computationalefficiency. To check the effect of the ThisworkwassupportedinpartbyNSFGrantNo. DMR-
interpolationaccuracy,in Fig. 5 we plotthe relativecompu- 0103003.
∗ [e-mail] 1 L. P. Kouwenhoven, C. M. Marcus, P. L. McEuen, S. Tarucha,
† [e-mail] R. M. Wetervelt, and N. S. Wingreen, in Mesoscopic electron
7
transport,editedbyL.L.Sohn,G.Scho¨n,andL.P.Kouwenhoven 18 H.Jiang, H.U.Baranger,andW.Yang, Phys.Rev.Lett.(2002),
(Kluwer,Dordrecht,1997),pp.105–214. accepted.
2 Y.Alhassid,Rev.Mod.Phys.72,895(2000). 19 D.UllmoandH.U.Baranger,Phys.Rev.B64,245324(2001).
3 I.L.Aleiner,P.W.Brouwer,andL.I.Glazman,Phys.Rep.358, 20 G.UsajandH.U.Baranger,Phys.Rev.B66,155333(2002).
309(2002). 21 M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, and J. D.
4 R.G.ParrandW.Yang,Density-FunctionalTheoryofAtomsand Joannopoulos,Rev.Mod.Phys.64,1045(1992).
Molecules(OxfordUniversityPress,NewYork,1989). 22 G.KresseandJ.Furthmu¨ller,Phys.Rev.B54,11169(1996).
5 R.O.JonesandO.Gunnarsson,Rev.Mod.Phys.61,689(1989). 23 J.R.Chelikowsky,N.Troullier,K.Wu,andY.Saad,Phys.Rev.B
6 R. M. Dreizler and E. K. U. Gross, Density Functional Theory 50,11355(1994).
: An Approach to the Quantum Many-Body Problem (Springer- 24 M.P.Teter,M.C.Payne,andD.C.Allan,Phys.Rev.B40,12255
Verlag,Berlin,1990). (1989).
7 A.Kumar,S.E.Laux,andF.Stern,Phys.Rev.B42,5166(1990), 25 G. J. Martyna and M. E.Tuckerman, J. Chem. Phys. 110, 2810
5166-5175. (1999).
8 M.Macucci, K.Hess, andG.J.Iafrate,Phys.Rev.B48, 17354 26 W.H.Press,B.P.Flannery,S.A.Teukolsky,andW.T.Vetterlin,
(1993). NumericalRecipes: TheArtofScientificComputing(Cambridge
9 D.JovanovicandJ.P.Leburton,Phys.Rev.B49,7474(1994). University,Cambridge,England,1989).
10 M.Stopa,Phys.Rev.B54,13767(1996). 27 I. Stich, R. Car, M. Parrinello, and S. Baroni, Phys. Rev. B 39,
11 S.Nagaraja,P.Matagne,V.Y.Thean,J.P.Leburton,Y.H.Kim, 4997(1989).
andR.M.Martin,Phys.Rev.B56,15752(1997). 28 M.J.Gillan,J.Phys.:Condens.Matter1,689(1989).
12 M.Koskinen,M.Manninen,andS.M.Reimann,Phys.Rev.Lett. 29 G.Onida,L.Reining,R.W.Godby,R.DelSole,andW.Andreoni,
79,1389(1997). Phys.Rev.Lett.75,818(1995).
13 I.-H.Lee,V.Rao,R.M.Martin,andJ.-P.Leburton,Phys.Rev.B 30 L.M.Fraser,W.M.C.Foulkes,G.Rajagopal,R.J.Needs,S.D.
57,9035(1998). Kenny,andA.J.Williamson,Phys.Rev.B53,1814(1996).
14 S. Bednarek, B. Szafran, and J. Adamowski, Phys. Rev. B 64, 31 M.R.Jarvis,I.D.White,R.W.Godby,andM.C.Payne,Phys.
195303(2001). Rev.B56,14972(1997).
15 I.Yakimenko,A.M.Bychkov,andK.-F.Berggren,Phys.Rev.B 32 R.W.Hockney,MethodsComput.Phys.9,135(1970).
63,165309(2001). 33 I.-H.Lee, Y.-H.Kim, and R.M.Martin, Phys.Rev. B61, 4397
16 M.Pi,A.Emperador,M.Barranco,andF.Garcias,Phys.Rev.B (2000).
63,115316(2001). 34 B.TanatarandD.M.Ceperley,Phys.Rev.B39,5005(1989).
17 K.HiroseandN.S.Wingreen,Phys.Rev.B65,193305(2002).