A novel multigrid method for electronic structure calculations

9 Pages·0.37 MB·English

by  M. Heiskanen

A novel multigrid method for electronic structure calculations
M. Heiskanen, T. Torsti, M. J. Puska, and R. M. Nieminen
Laboratory of Physics, Helsinki University of Technology, P.O, Box 1100, FIN-02015 HUT, FINLAND
(February 1, 2008)
Ageneralreal-spacemultigridalgorithmfortheself-consistent solutionoftheKohn-Shamequa-
1
tions appearing in the state-of-the-art electronic-structure calculations is described. The most im-
0
portant part of the method is the multigrid solver for the Schr¨odinger equation. Our choice is the
0
Rayleigh quotient multigrid method (RQMG), which applies directly to the minimization of the
2
Rayleigh quotient on the finest level. Very coarse correction grids can be used, because there is no
n
need to be able to represent the states on the coarse levels. The RQMG method is generalized for
a
thesimultaneoussolutionofallthestatesofthesystemusingapenaltyfunctionaltokeepthestates
J
orthogonal. The performance of thescheme is demonstrated by applying it in a few molecular and
6
solid-state systems described by non-local norm-conserving pseudopotentials.
1

]
i
c
s
- I. INTRODUCTION wave functions. In the RS methods periodic boundary
l conditions are not necessary. This leads to ease and ac-
r
mt One of the goals of computational materials science curacy in describing charged atomic clusters in contrast
is to calculate from first principles the various physical to PW methods requiring an artificialneutralizing back-
t. and chemical properties. This requires the solution of ground charge. Besides the above “physical arguments”
a
m the electronic and ionic structures of the materials sys- therearealsomethodologicalandcomputationalaspects
tem in question. The density-functional theory (DFT) favouringthe RSmethods. TheRSmethodsallowasys-
d- makes a huge step towards this goal by casting the un- tematic convergence control by increasing the grid (or
n tractable problem of many interacting electrons to that basis function) density. (The PW methods do also so by
o of noninteracting particles under the influence of an ef- adjustingthecutoffenergyoftheplanewaveexpansions.)
c fective potential1. The adiabatic approximation allows The so-called “order-N” methods12, the computational
[ one to separate the ionic degrees of freedom from those cost of which scales linearly with the number of elec-
1 of the electrons. However, in order to apply DFT in trons,requirelocalizedreal-spacewavefunctions leading
v practice one has to resortto approximationsfor electron naturally to the employment of RS methods13. The dis-
3 exchange and correlation such as the local-density ap- cretizations in the real-spacegrid can be made localand
3 proximation (LDA) or the generalized-gradient approxi- therefore parallelization can effectively use data decom-
2
mation(GGA). Moreover,inthe caseofsystemsconsist- positioninwhichdifferentreal-spaceregionsarehandled
1
ing of hundreds or more atoms it is still a challenge to with different processing units and the communications
0
1 solve numerically efficiently for the ensuing Kohn-Sham between processing units will be mainly short-ranged14.
0 equations. More specifically, our choice for the numericalmethod
/ The numerical solution of the Kohn-Sham equations is a multigrid scheme15,2. Several approaches employing
t
a is the concern of our present work. It deals with real- themultigridideawithinelectronicstructurecalculations
m space (RS) methods, in which the values of the dif- have appeared during recent years16,14,17,8,18. The main
- ferent functions are presented using three-dimensional idea of multigrid methods is that they avoid the critical
d
pointgrids,andthepartialdifferentialequationsaredis- slowing-down (CSD) phenomenon occuring when a par-
n
cretized using finite differences2,3. The RS methods, as tial differential equation discretized on a real space grid
o
c suggested by the name chosen, are contrasted with the is solved with a simple relaxation method such as the
: popular plane-wave (PW) schemes4,5. There are several Gauss-Seidel method. The discretization operators typi-
v
aspects favouring the RS methods over the PW meth- callyuseinformationfromaratherlocalizedregionofthe
i
X ods. Both of the methods are used in the context of gridat a time. Therefore the high frequency errorof the
r pseudopotentialsdescribingtheelectron-ioninteractions, lengthscaleofthegridspacingisreducedveryrapidlyin
a
but only the RS can easily be used in all-electron calcu- the relaxation. However, once the high frequency error
lations or with hard pseudopotentials of, i.e. first-row has effectively been removed, the very slow convergence
or transition metal atoms, because the RS grid can be of the low frequency components dominates the overall
refinedin a naturalwayin the ioncoreregions(compos- error reduction rate15, i.e. CSD occurs. In multigrid
ite grids6–8, adaptive coordinates9–11). Systems, suchas methods one stops the relaxation on a given (fine) grid
surfaces, containing different length scales are more eco- beforeCSDsetsinandtransferstheequationtoacoarser
nomically described in the RS than in the PW scheme grid (the so-called restriction operation) where the low-
because one needs not waste many grid points in the frequencycomponentscanbe solvedmoreefficiently. On
vacuum regions to describe the slowly varying tails of the coarsest grid the problem is solved exactly or as ac-
1
curately as possible, after which one interpolates (the can accelerate the convergence remarkably with respect
so-called prolongation operation) the correction to finer to the (above-described) V-cycle scheme in which one
grids, performing simultaneously relaxations in order to starts from the finest level. Fattebert8 used a multigrid
removethehigh-frequencyerrorsintroducedintheinter- method with a block Galerkin inverse iteration (BGII)
polation. and GMRES in the relaxations. In the method, the
The solution of the Poisson equation by multigrid current approximation is kept orthogonalagainst all the
methods is straightworward15. This is because the error nearby states during the multigrid cycle. The inverse it-
(or the correctionneeded) also obeys a Poissonequation eration converges for a given guess for the energy eigen-
and thus will be a smooth function to be presented and value towards the nearest eigenvalue. In order to solve
solved on the repeatedly coarser grids optimal to handle all the desired lowest eigenstates a good guess for the
thelowerfrequencies. Thesolutionofaneigenvalueprob- eigenvalue spectrum is needed in the beginning of iter-
lem, such as the Schr¨odinger equation, is a much more ations, but thereafter large computational savings may
complicatedtaskthanthatofthe Poissonequation. The be expected because explicit orthogonalizations are not
problem is no more linear because both the eigenfunc- needed (at least between well-separated states).
tionandtheeigenvaluehavetobesolvedsimultaneously. Asevereprobleminthe existingmultigridschemesfor
Then the error no longer obeys the same equation as the Schr¨odinger equation is often that the coarse grids
the solution. Also one has to solve for several eigenpairs cannot well approximate the solutions of the coarse grid
(eigenvalues and corresponding eigenvectors). Moreover, equations themselves. As a consequence the correction
the existence of both negative and positive eigenvalues fromcoarsegrids,nomatterhowaccuratelytheequation
makes the problem indefinite. This implies severe dif- issolved,maybeineffectiveincorrectingthefinegridso-
ficulties for many simple iterative methods which con- lutionandasaresulttheoverallprocessconvergesslowly.
verge only in the case of a positive definite iteration ma- Thereforeonehastorestricttotheuseofratherfinegrids
trix. In particular, it can easily be shown that using only and the convergence speed of the scheme is drasti-
Gauss-Seidelrelaxationfor the Schr¨odingerequationthe callylowered. Inthosemultigridmethods,whichusethe
high frequency components typically converge as in the potential also on the coarsegrids the size of the coarsest
caseofthePoissonequation,butthe lowfrequencycom- grid has been typically of the order of 31x31x318,17. Us-
ponents may diverge, although the divergence may be ingtheFASmethodcoarsergridsarepossibleatleastfor
slow19. MorecomplicatedmethodssuchasKaczmarzre- systems with a small number of eigenstates solved18. If
laxationareguaranteedtoconverge,butmayhaveclearly alargenumberofeigenstateshavetobe solvedproblems
inferior high frequency reduction rates, which are essen- may arise because the coarse grids may not be able to
tial for the overall speed of multigrid methods. Other representeigenstateswithmanynodesorthe orderingof
possible convergentmethods include GMRES20 which is the states may change between the successive grids. To
considerablymorecomplexthanGauss-Seidelrelaxation. bypasstheseproblemsinFAS,rathercomplicatedstrate-
A standard recipe for dealing with eigenproblems gies are needed23.
with multigrids is the full-approximation-storage (FAS) In order to avoid the coarse grid representation
method originally described by Brandt21. In FAS one problems we utilize the so-called Rayleigh Quotient
solvesforthe entireproblemonthecoarsegridsalsoand Multigrid (RQMG) method introduced by Mandel and
endsupinsolvingforaproperlymodifiedproblemsothat McCormick24. In this method the coarse grid relaxation
its solution can be used in correcting the fine grid solu- passesareperformedsothattheRayleighquotientcalcu-
tion. The FAS method may not be very straighworward latedonthefinegridwillbeminimized. Inthiswaythere
to implement for the Schr¨odinger equation. It is also is no requirement for the solution to be well represented
difficult to present some actual potential on the coarse ona coarsegridandthe coarsegridrepresentationprob-
levelsaccuratelyenough. However,somesuccesfulappli- lem is avoided. Mandel and McCormick24 introduced
cationsofFAShaveappearedinthecontextofelectronic the method for the solution of the eigenpair correspond-
structure calculations by Beck et al.22,18 and advanced ing to the lowest eigenvalue. We have generalized it to
strategies for FAS have been proposed23. the simultaneous solution of a desired number of lowest
Briggs et al.16,14 employ a multigrid method in eigenenergy states by developing a scheme which keeps
electronic structure calculations by linearizing the the eigenstates separated by the use of a penalty func-
Schr¨odinger problem and presenting the potential con- tional,Gramm-Schmidtorthogonalization,andsubspace
tribution onthe coarselevels by anerrorterm (residual) rotations. OurgeneralizationoftheRQMGmethodisan
only. Then on the coarse levels they solve effectively attractive alternative for large-scale electronic structure
for the Poisson problem. Ancilotto et al.17 modified the calculations.
method by Briggs et al. by shifting to a full multi- The Kohn-Sham equations have to be solved self-
grid (FMG) scheme and by solving on the coarse grids consistently, i.e. the wave functions solved from the
a problem including a local potential term. The idea of single-particle equation determine via the density (so-
FMG is to start the smoothing iterations from a coarse lution of the Poisson equation and the calculation of
grid. Then the interpolation to a finer grid provides a the exchange-correlation potential) the effective poten-
good initial guess of the solution. The FMG scheme tial for which they should again be solved. To approach
2
this self-consistency requires an optimized strategy so consists of an external potential V (r) due to ions (or
ion
that numerical accuracy of the wave functions and the nuclei in all-electron calculations), the Hartree potential
potential increase in balance, enabling the most effi- V (r) calculated from the electron density distribution,
H
cient convergence18. In order to avoid the divergence andtheexchange-correlationpotentialV (r). Intheex-
XC
ofthe self-consistencyiterations,the mixing ofthe input amples of the present work we use the norm-conserving
and output solutions is needed. For this feedback pro- non-local pseudopotentials for the electron-ion interac-
cedure sophisticated schemes25 and control strategies11 tions andthe local-densityapproximation(LDA) for the
have been presented. exchange-correlationenergy
The outline of the present paper is as follows. In Sec-
tion II we represent shortly the most important ideas of E [n(r)]= ǫ (n(r))n(r)dr, (6)
XC XC
the density functional theory. Section III is devoted for Z
numerical methods, the most important of which is the
and for the exchange-correlationpotential
Schr¨odinger equation solver developed. Also the strat-
egy for the self-consistency iterations is discussed. In dǫ
V (r)=ǫ (n(r))+n(r) XC . (7)
Section IV we demonstrate by the help of a couple of XC XC dn |n=n(r)
examples the performance of our scheme in calculating
the electronic structures of small molecules and solid- TheHartreepotentialissolvedfromthePoissonequa-
state systems described by pseudopotentials. Section V tion
summarises the work and gives outlines for the future
∇2V (r)=−4πn(r). (8)
developments. H
In practice, the electron density n(r) is substituted by
the totalchargedensity ρ(r), whichincludes the positive
II. THE KOHN-SHAM SCHEME
ionic (nuclear) charge neutralizing the system. In the
case of finite systems, Dirichlet boundary conditions are
In the Kohn-Sham method for electronic structure
used with the Coulomb potential values calculated using
calculations one solves for a set of equations self-
a multipole expansion. For periodic systems we fix the
consistently1. In the following, we present the equations
average Coulomb potential to zero and allow the peri-
inthespin-compensatedform. Inpractice,wehavemade
odic boundary conditions to result in the corresponding
the straightforwardgeneralizationusing the spin-density
convergedpotential.
functional theory. The set of equations reads as (atomic
The self-consistent solution of the above Kohn-Sham
units with h¯ =m =e=1 are used):
e equations leads to the ground state electronic structure
minimizing the total energy
1
− ∇2+V (r) Ψ =ǫ Ψ , (1)
eff i i i
(cid:18) 2 (cid:19) 1 1
E = Ψ∗(r) − ∇2 Ψ (r)dr+ V (r)n(r)dr
tot Z i (cid:18) 2 (cid:19) i 2Z H
Xi
N
n(r)= |Ψ (r)|2, (2) + Vion(r)n(r)dr+EXC+Eion−ion , (9)
i Z
Xi
where E is the repulsive interaction between the
ion−ion
ions(nuclei)ofthesystem. Insteadoftheself-consistency
Veff(r)=Vion(r)+VH(r)+VXC(r), (3) iterationsthesolutionoftheKohn-Shamproblemcanbe
found by minimizing directly the total energy with re-
n(r′) spect to the wave function parameters, e.g. plane-wave
V (r)= dr′, (4) coefficients4. However, Kresse and Furthmu¨ller5,25 have
H Z |r−r′|
found this scheme less efficient than the self-consistency
iterations.
δE [n(r)]
V (r)= XC . (5)
XC δn(r)
III. NUMERICAL METHODS
The first equation (1) is a Schr¨odingerequation for non-
interacting particles in an effective potential V (r). For A. Schr¨odinger equation solver
eff
finite systems the wave functions are required to van-
ish at the boundaries of the computation volume. In Inourrealspacemethodwestartfromaninitialguess
the case of infinite periodic systems the complex wave for the effective potential and initial wave functions gen-
functions have to obey the Bloch theorem at the cell erated by random numbers in grid points. The wave
boundaries. Theelectrondensityn(r)isobtainedfroma functions and the Hartree potential are updated alter-
sum over the N occupied states. The effective potential natingly towards self-consistency. The solution of the
3
Poisson equation is a standard task for the multigrid the Rayleigh quotient calculated on the fine grid will be
scheme. If a reasonable guess for the Coulomb poten- minimized. In this way there is no requirement for the
tial is not available, the FMG method will produce the solutionto be wellrepresentedona coarsegrid. Inprac-
solutionstartingfromrandomnumbersandrequiringthe tice,acoarsegridsearchsubstitutesthefinegridsolution
workwhichscaleslinearly asafunction ofthe size ofthe by
system (O(N)). During the Kohn-Sham iterations one
u′ =u +αIfd , (13)
can start from the present approximation of Coulomb f f c c
potential and update it with respect to the new charge
wherethesubscriptsf andcstandforthefineandcoarse
density by performing only a few V-cycles.
grids, respectively, and If a prolongation operator in-
Thesolutionofthewavefunctionsisamuchmorecom- c
terpolating the coarse grid vector to the fine grid. The
plicated task than that of the Poisson equation because
Rayleigh quotient to be minimized is then
onehastosolveaneigenvalueproblemwhichinthestate-
of-the-artelectronicstructurecalculationsmeansthe de- huf+αIcfdc|Hf|uf+αIcfdci =
termination of several hundreds of eigenpairs. For this huf+αIcfdc|Bf|uf+αIcfdci
purpose we have developed a scheme based on RQMG huf|Hfufi+2αhIfcHfuf|dci+α2hdc|Hcdci. (14)
method introduced by Mandel and McCormick24 for the huf|Bfufi+2αhIfcBfuf|dci+α2hdc|Bcdci
solution of the eigenpair corresponding to the lowest
The second form is obtained by relating the coarse grid
eigenvalue. We begin by reviewing the basic principles
operators, H and B , with the fine grid ones, H and
of RQMG. This is most easily done in the framework of c c f
B , by the Galerkin condition
the so-called coordinate relaxation method. Thereafter f
we go throughthe modifications made in order to simul- H =IcH If
taneously solve for several eigenpairs. c f f c (15)
B =IcB If,
Coordinate relaxation is a method of solving the dis- c f f c
cretized eigenproblem and the restriction operator Ic has to be the transpose
f
of the prolongation operator
Hu=λBu (10)
Ic = If T . (16)
by minimizing the Rayleigh quotient f c
(cid:0) (cid:1)
The key point to note is that when H u and B u are
hu|H|ui f f f f
. (11) provided from the fine grid to the coarse grid, the re-
hu|B|ui
maining integrals can be calculated on the coarse grid
itself. Thus one really applies coordinate relaxation on
Above,H andB arematrixoperatorschosensothatthe
the coarsegrids to minimize the fine level Rayleighquo-
Schr¨odinger equation discretized on a real-space point
grid with spacing h is satisfied to a chosen order O(hn). tient. This is a major departure from the earlier meth-
ods,whichtosomeextentrelyontheabilitytorepresent
In Eq. (11) u is a vector containing the wave function
the solution of some coarse grid equation on the coarse
values at the grid points. In the relaxation method, the
grid itself. Here, on the other hand, one can calculate
currentestimate u is replacedby itself plus a multiple of
the exact change in the Rayleigh quotient due to any
some search vector d
coarse grid change, no matter how coarse the grid itself
u′ =u+αd, (12) is. Thereisnoequationwhosesolutionwouldhavetobe
representable.
and α is chosento minimize the Rayleighquotient. This Thus,inthe Rayleighquotientminimizationmultigrid
leads to a simple quadratic equation for α. (Find the (RQMG) algorithm the coordinate relaxation passes on
minimumoftheexpression(14)belowwithrespecttoα. each level keep track of the integrals in Eq. (14). Actu-
In the case of a complex wave function one has to solve ally, on the finest level we use Gauss-Seidel relaxation,
fortherealandimaginarypartsofαfromacoupledpair which very effectively smoothens the errors of the wave-
ofquadraticequations.) Moreover,the searchvectord is length corresponding to the grid spacing. When calcu-
simply chosen to be unity in one grid point and to van- latingseveraleigenpairsGauss-Seidelrelaxationmayalso
ish in all other points. A complete coordinate relaxation workasaresidualminimizationmethod. Theideaisthat
pass is then obtained by performing the minimization at the coarsegrid-iterationswithGramm-Schmidtorthogo-
eachpointinturnandthesepassescanberepeateduntil nalizationcanprovidetheseparationoftheeigenstatesso
the lowest state is found with desired accuracy. wellthatthe subsequentfinestlevelrelaxationsconverge
Naturally, also the coordinate relaxation suffers from to the correct (nearest) eigenstates without orhogonal-
CSD because of the use of local information only in up- ization. This requires that the effect of the coarse-level
dating u in a certain point. In order to avoid it one smoothingsonthelow-frequencycomponentsofthesolu-
applies the multigrid idea. In the multigrid scheme by tionsovercomesthepossibledivergencetendencyofthese
MandelandMcCormick24thecrucialpointisthatcoarse componentscausedbytheGauss-Seidelrelaxationonthe
grid coordinate relaxation passes are performed so that finest level.
4
Moreover, we discretize the original equation sepa- future parallelizationoverthe eigenstates. Thus the cur-
rately on each grid (discretization coarse grid approxi- rent approximations are used for u , i = 1...k. More-
i
mation (DCA)) instead of using the Galerkin conditions over, changes in the u during a given relaxation sweep
i
of Eq. (15) This may in principle decrease the conver- are not used to update the penalty term in Eq. (17).
gence rate and force a limit to the coarsest possible grid This is sufficient, when the states are always ordered in
in order to avoid instability or divergence. However, we the same way, in the order of increasing eigenvalue. In
haveobservedthisDCAimplementationofRQMGtobe order to reduce computations, the B-innerproduct is ac-
quite stable and efficient. To avoid possible coarse-level tually used in calculating the penalty term integrals be-
instabilities occuring especially during the first few it- causethevaluesofBuarereadilyavailablefromthefiner
eration cycles we may recalculate the Rayleigh quotient level. The substitution (13) is introduced in the func-
whenever coarse grid corrections are interpolated to a tional (17)andtheminimizationwithrespecttoαleads
finer grid. Later when approaching the convergence the againto a quadratic equation. This time the coefficients
recalculation can be omitted. contain terms due to the penalty part.
For the matrix operators H and B we have used ei- On the finest level, we do not apply the minimization
ther high-order (O(h4) or higher) Mehrstellen or central of the penalty functional. The ideal situation would be
differencestencils(CDS)14,8. Theuseofhigh-ordersten- if a residual minimization method, such as the Gauss-
cilsreducesremarkablythedensityofgridpointsneeded. Seidel method, would keep the states calculated on the
The benefit of the Mehrstellen scheme is that more local coarselevels separated. We found out in practicalcalcu-
information is used. The scheme leads to controlled ac- lations that this is not true at least when the states are
curacyand convergencepropertiesand to more isotropic far from convergence. Therefore we have developed for
smoothing of the error in comparison with the use of the finest level a scheme which by employing Gramm-
CDS’s. Thelocalnatureenablesalsoamoreefficientpar- Schmidt orthogonalization and subspace rotation keeps
allel coding. As the prolongationoperator If we usually the eigenstates orthogonal. The subspace rotation is a
c
usetrilinearinterpolationandastherestrictionoperator method to find the most optimally separated eigenvec-
Ic its transpose, the so-called full-weighting operator,in tors from the approximative ones. The major steps of
f
which the coarse-grid values are chosen to be the av- the rotation are:
eraged values of the surrounding fine grid points. The (i)CalculationoftheHamiltonianmatrixelementsbe-
integrations are performed by the trapezoidal rule. tween the current states:
Next we consider the generalization of the RQMG
H¯ =hu |B−1H|u i. (19)
method to the simultaneous solution of several (N) mu- i,j i j
tually orthogonal eigenpairs. The separation of the dif- (ii) Calculation of the overlap matrix:
ferent states is divided into two or three subtasks. First,
in orderto make the coarsegridrelaxationsconvergeto- S¯ =hu |u i. (20)
i,j i j
wards the desired state we apply a penalty functional
TheuseofmatrixelementsofEqs. (19)and(20)leadsto
scheme. Giventhek lowesteigenfunctions,thenextlow-
eigenvectorsorthogonalinthedesiredEuclidiansense(I-
est, (k + 1)’th state is searched for by minimizing the
orthogonal)and not in the sense of the B-innerproduct.
functional
(iii) Diagonalization to find the optimal eigenvectors
hu |H|u i k hu |u i2 (u′k = jA¯k,juj) and corresponding eigenvalues (λk):
k+1 k+1 + qi i k+1 . (17) P
huk+1|B|uk+1i Xi=1 hui|uii·huk+1|uk+1i H¯i,jA¯j,k =λk S¯i,jA¯j,k. (21)
Xj Xj
The overlap integral in the penalty term is squared to
make the penalty positive definite. The denominator is In practice, we apply the approximation
requiredtomakethefunctionalindependentofthenorms
hu |Hu i
ofui, i=1...k+1. The minimizationofthis functional hui|B−1H|uji≈hui|ujihui|Buji. (22)
is equivalent to imposing the orthonormality constraints i j
against the lower k states, when qi →∞. By increasing The Gramm-Schmidt orthogonalization and the sub-
the shifts qi any desired accuracy can be obtained, but space rotations are organized so that the space of the
in order to obtain a computationally efficient algorithm eigenvectorsis first divided to small clusters correspond-
a reasonable finite value should be used, for example ing to close eigenvalues. The Gramm-Schmidt orthogo-
nalization is then performed for each cluster at a time
qi =(λk+1−λi)+Q, (18) so that its eigenvectors become orthogonal against the
eigenvectors of the clusters of lower eigenvalues. Then a
where Q is a sufficiently large positive constant. In our subspace rotationis performed within the states belong-
test calculations Q is of the order of Q=0.5...2 Ha. ing to the present cluster. The division to clusters re-
We minimize the expression (17) simultaneously for ducesremarkablythecostofthesubspacerotation. This
all N states. This simplifies the algorithmand enables a is because the costis proportionalto O(N3), where N is
5
thenumberofstatesrotated. Moreover,thesubspacero- level. The accuracyofthe wavefunctionsis controlledby
tation requires the calculation of matrix elements which calculating the norms of the residuals of the eigenstates
are more complicated than those for the simple Gramm- anditis finally improvedby adding more V-cyclesstart-
Schmidt orthogonalization. ing from the finest level. A certain accuracy is needed
According to our test calculations this subspace rota- in order to initiate self-consistency iterations which con-
tionscheme leadsquite effectively to I-orthogonaleigen- verge without large density oscillations. Then the new
states. This is seen as a convergence of the eigenvalue electron density and the ensuing effective potential are
problem within the numerical accuracy, i.e. the residu- calculated. Thenew potentialisnotdirectlyfedintothe
als of different eigenstates vanish. In order to prove this next iteration but it is mixed in this place, as well as
one has, in practice, to introduce potential shifts which laterbetweenthe self-consistencyiterations,withthe in-
reducethe number ofsignificantdigits inthe eigenvalues put potential of the iteration. We monitor the accuracy
so that the numerical accuracy of the eigenvalue does of the wave functions by calculating their residuals and
notpreventtoreachthenumericalaccuracyofthewave- require that the accuracy has improved from the previ-
function, i.e. the vanishing residual. The error in the ous iteration. Usually one V-cycle is sufficient for this,
eigenvaluescalesasthesquareofthe residual. Whenap- because the changes in the potential are small.
plying the subspace rotation it is important to completPeSfrag repla
ementsPSfrag repla
ements
thehighesteigenvaluecluster;otherwisetherotationmay
h h
become inefficient. Non self-
PoSnfsraisgNtreeonpntla
sseeomllfue-n
ttoisonnsistent solution
The orthogonalization needed scales as O(N3). For Non self-
onsistent solution
small systems of severaltens of eigenpairsthis is not yet
a problem. The algorithm is effective and the number h
of fine grid orthogonalizations remains quite plausible,
for example, in comparison with the conjugate gradient 2 h 2 h 2 h
searchofeigenpairsemployingonlythefinestgrid26. But 4 h 4 h
4 h
for larger systems with hundreds of states it will be the 8 h 8 h
bottleneck. One solution could be to rely on the finest 16 h 16 h
8 h
level only on a residual minimization method when the
initialstagesoftheiterationprocesshavebeenperformed
16 h
andthe solutionis clearlyona stable tSraeclkf-t
oownasrdisstceonn
-y iterations
vergence.
B. Strategy for self-consistency iterations
Self-
onsisten
y iterations
The Kohn-Sham problem has to be solved self-
P P
consistently. This means that an optimal strategy is
neededsothatcomputingtimeisnotwastedinthebegin-
ning of the self-consistency iterations to obtain unneces-
sarilyaccuratewavefunctions, becausethesewillchange :::
during the later iterations due to the changes in the po-
tential. Updating thepotential,including the solutionof
thePoissonequation,isamuchlesstime-consumingtask
than the update of all the wavefunctions. Therefore the
potential update can be performed frequently18.
FIG. 1. Strategy of self-consistency iterations. First, the
The examples of this paper are small-molecule and wavefunctions are solved nonselfconsistently using the full
bulk-solid systems described by pseudopotentials. The multigridmethodintheinitialpotentialcorrespondingtothe
strategy used is schematically presented in Fig. 1. Sim- superpositionofpseudoatoms. Thentheeffectivepotentialis
ilar strategies can certainly be applied in other kind of updated (this is denoted by P in the figure). The potential
Kohn-Sham calculations, for example in those employ- update amounts to calculation of the new electron density,
ing all-electron or jellium-type models. In the exam- the solution of the Poisson equation and calculation of the
ples of this work the initial electron density is the su- new exchange correlation potential. Next the wave-functions
perpositionof the pseudoatomdensities centeredaround are updated by one V-cycle. These two steps are repeated
given nuclear positions. From the superposition we cal- untilself-consistency has been reached.
culate the initial effective potential, where the wave-
functions are solved accurately enough using the full- An important point is also to find a proper balance
multigrid method. The FMG process is started from with respect to the pre- and postsmoothening sweeps on
random numbers for the wavefunctions on the coarsest the different grid levels. Typically, on the finest level
6
we made two pre- and postsmoothening sweeps whereas we have found that this may increase the convergence
on the coarser grids their number is four. Actually, this rate and stability of the method.
means that on the finest level four successive sweeps are The CO molecule is placed diagonically in the center
2
done if the potential is not updated. A potential up- of a cubic computation volume of the size of (12.6 a )3.
0
date is always preceded and followed by two immediate Experimentalbondlengthsareused. Dirichletboundary
smoothening sweeps. conditions are used so that the potential values outside
the cube are obtained from a multipole expansion of the
charge density. The point mesh used is 633, giving the
IV. TESTS gridspacing h = 0.20a . The Mehrstellen discretization
0
by Briggs et al.14 is used.
We test the performance of our scheme by calculating Inthiscalculationweusedamixingscheme,wherethe
theself-consistentelectronicstructuresofaCO2molecule new effective potential Viin+1 is obtained from the input
aswellasthatofperfectbulkSilatticewithasupercellof and output potentials according to
64Siatoms. Theformersystemisanexampleoftheem-
Vi+1 =(1−κ)Vi +κVi . (25)
ploymentofDirichletboundaryconditionsandtheuseof in in out
“hard” pseudopotentials whereas the latter system rep-
resents the use of periodic boundary conditions and a
supercell size typical in electronic structure calculations
for point defects in solids. 100
The ions are described by pseudopotentials of the V)
Kleinman-Bylander form27, (e
E
V (r)= V (|r |) NC10−2
ion ion,loc a E
Xa R
E
1
+a,Xn,lm ∆Vlam (cid:12)∆Vion,l(ra)ulm(ra)(cid:11) DIFF10−4
× ∆Vion(cid:10),l(ra′)u(cid:11)lm(cid:12) (r′a) , (23) GY
R
(cid:10) (cid:12)
where ∆Vlam is a normalization fac(cid:12)tor, ENE10−6
(cid:10) (cid:11)
∆Va = u (r )∆V (r )u (r )d3r, (24)
lm lm a ion,l a lm a
Z
(cid:10) (cid:11) 0 5 10 15
SCF−ITERATION
and ra = r−Ra, ulm are the atomic pseudopoten- FIG. 2. Convergence of the total energy for the
tialwavefunctions ofangularandazimuthalmomentum CO2-moleculeusingdirectmixingwithdifferentvaluesofthe
quantum numbers (l,m), from which the l-dependent feedback parameter κ; κ = 0.4 (solid line), κ = 0.5 (dashed
ionic pseudopotentials Vion,l(r) are generated using the line), κ=0.6 (dashdotted line) and κ=0.7 (dotted line). A
Troullier-Martins scheme28. The ion core is assumed horizontallinehasbeenaddedtoindicatethechemicalaccu-
to be spherically symmetric. ∆V (r) = V (r) − racy of 1 meV.
ion,l ion,l
V (r) is the difference between the l-component of
ion,loc
the ionic pseudopotential and the local ionic potential. The convergence of the self-consistency iterations em-
We have chosen the s-component of the pseudopotential ploying the strategy described above (Fig. 1) is shown
as the local component. in Fig. 2. The deviation of the total energy from the
Because the functions u (r) are short-ranged, oper- convergedvalue is given as a function of self-consistency
lm
ating on the wave-function by the nonlocal parts of the iteration steps performed. The zeroth iteration is a full-
pseudopotentialisinpracticeamultiplicationbyasparse multigridsolutionforthewavefunctionsintheinitialpo-
matrix. The numerical work required to compute this tential. Two V-cycles are performed on the finest level
scales as the square of the number of atoms in the sys- at this point. The effective potential obtained from the
tem,whereasintheconventionalreciprocal-spaceformu- output electron density was mixed with the initial po-
lation the work scales as the cube of the system size. tential using the feedback κ = 0.4. Next, at iteration
The advantage of implementing the nonlocal pseudopo- one,the wave-functionsare relaxedin this new potential
tentials in real space has been noted also in the context using one V-cycle. From this point on, the four curves
of plane-wave methods29. in the figure give the convergence with different values
Inthepreviousmultigridimplementationsofthepseu- of the feedback parameter κ. One V-cycle per one self-
dopotential method17,14, the nonlocal parts have only consistency iteration step is done. A wide range of val-
beenemployedonthefinestgrid. Itis,however,straight- ues for κ gives satisfactory convergence indicating a ro-
forwardtoimplementthemalsoonthecoarselevels,and bust behaviour for the scheme. The accuracy of 1 meV,
7
whichissufficientinpracticalcalculations,isreachedaf- degenerate eigenenergies with an accuracy better than 1
ter three or four V-cycles. The implementation of the meV. The convergence towards to the self-consistent so-
non-local parts of the pseudopotential on the coarse lev- lutionoccurssimilarlyasfortheCO moleculeinFig. 2.
2
elsisfoundtospeeduptheconvergenceespeciallyinthis Thus, the convergence process seems to be independent
region. From Fig. 2 we obtain an average convergence of the size of the system.
rate of approximately one decade per self-consistency it-
eration. This is of the same order as those reported by
Wang and Beck18 in their FAS scheme or by Kresse and
Furthmu¨ller25 in their pseudopotential scheme employ- V. SUMMARY AND OUTLOOK
ing self-consistency iterations. The convergence rate of
one decade per self-consistency iteration is better than
In this work we have generalized the RQMG method
that obtained by Ancilotto et al.17 in the FMG scheme
introducedbyMandelandMcCormick24 forthesimulta-
and much better than the rate reached in the linearized
neous solutionofa desirednumber oflowesteigenenergy
multigrid scheme by Briggs et al.14.
states. This approach can be viewed as belonging to a
thirdgroupofmultigridmethods,inadditiontoFASand
the techniques where the eigenproblemis linearized(e.g.
inverse iteration. In principle, one can use arbitrarily
coarse grids in RQMG, whereas in the other multigrid
methods one has to be able to representallthe stateson
the coarsest grid.
We have demonstrated the feasibility of the method
byelectronicstructurecalculationsfortheCO molecule
2
andbulk Si describedby pseudopotentials. Our strategy
fortheself-consistentsolutionconsistsofafull-multigrid
solution for the wave-functions in the initial potential,
andsubsequentself-consistencyiterations. Lessthanfive
V-cycles are generally sufficient for practically sufficient
accuracy. The cpu-times requiredfor the FMG andSCF
steps are roughly equal.
FIG. 3. Valence electron density in the (110)-plane ob-
We have applied the method also in two-dimensional
tained in the Γ-point calculation for the 64-atom supercell
of bulk Si. The area of the figure corresponds to the extent problems for quantum dots employing the current-spin-
of the supercell. density functional theory, in three-dimensional cylin-
drically symmetric systems, and also for calculation of
We have solved for the electronic structure of perfect positron states in solids.
Si lattice described by a supercell of 64 Si ions. The We believe that our method will eventually compete
lattice constant of 20.38 a used is the equilibrium value with the standard plane-wave methods for electronic
0
obtainedinaplane-wavecalculation,withwhichwehave structure calculations. However, some straightforward
comparedourreal-spaceresults. The first Brillouinzone programming is still required. For calculations, where
is sampled in this test using the Γ-point only. The point the optimization of the ionic structure is necessary, the
mesh used for the wave-functions is 643, giving the grid Hellmann-Feynman forces will be implemented. In or-
spacing h = 0.32 a . For the densities and potentials der to remove the spurious dependence of the total en-
0
we use a finer grid of 1283 points. The other numerical ergy on the position of the atoms with respect to the
parameters and the iteration strategy are the same as gridpoints,Fourier-filteringofthepseudopotentialsisre-
in the CO test. The resulting valence electron density quired. Complexwavefunctionsforanyk-pointareeasily
2
on the (110)-plane is given in Fig. 3. The area of the implemented, and are already in use in two-dimensional
figure corresponds to the extent of the supercell. One geometries.
notes that exactly the same features are reproduced at Parallelization over k-points can be done easily. One
theequivalentpointsindifferentregionsofthe supercell. only needs to communicate the electron density and ef-
Thismeansthatafully convergedresulthasbeenfound. fective potential at the end of each V-cycle. During
We have compared the results of our real-space code to the RQMG V-cycle, the states are all relaxed simulta-
thoseobtainedusingtheplane-wavemethod. Theenergy neously and independently of each other. Therefore par-
cutoff, 18Ry, of the plane-waveexpansionwas chosenso allelizationoverstates is naturalandeasy to implement.
that it results in a real-space point mesh of 643, i.e. it However,forlargersystemsthe Gramm-Schmidtorthog-
is the same as in our real-space calculation. The widths onalization becomes very inefficient in a state-parellel
of the valence band and band gaps obtained by the two code. The most efficient and yet straightforward choice
methods agree with an accuracy of 3 meV. In the case isreal-spacedatadecomposition,whereeachprocessoris
of degenerate eigenstates the real-space code results in mapped to a specific region of space.
8
ACKNOWLEDGMENTS gakis, E. Kaxiras, Modelling and Simulation in Materials
Science and Engineering, in print (cond-mat/0006183).
12S.Goedecker, Rev.Mod. Phys. 71, 1085 (1999).
M.H. acknowledges stimulating discussions with Prof.
13J.-L. Fattebert and J. Bernholc, Phys. Rev. B, 62, 1713
S. F. McCormick and Prof. T. L. Beck. T.T. acknowl-
(2000).
edges financial support by the Vilho, Yrj¨o and Kalle 14E.L.BriggsandD.J.SullivanandJ.Bernholc,Phys.Rev.
Va¨is¨al¨a foundation. This research has been supported
B, 54, 14362 (1996).
by the Academy of Finland through its Centre of Excel- 15A.Brandt, Math. Comp. 31, 333 (1977).
lence Programme (2000 - 2005). 16E.L.BriggsandD.J.SullivanandJ.Bernholc,Phys.Rev.
B, 52, 5471 (1995).
17F. Ancilotto, P. Blandin, and F. Toigo, Phys. Rev. B, 59,
7868 (1999).
18J.WangandT.L.BeckJ.Chem.Phys.112, 9223(2000).
19F. F. Grinstein and H. Rabitz and A. Askar, J. Comput.
Phys., 49, 4423 (1983).
1For a review, see R. O. Jones and O. Gunnarsson, Rev. 20Y. Saad and M. H. Schultz, SIAM J. Sci. Stat. Comput.,
Mod. Phys. 61, 689 (1989). 7, 856 (1986).
2T. L. Beck, Rev.Mod. Phys. 72, 1041 (2000). 21A.Brandt,S.F.McCormickandJ.W.Ruge,SIAMJ.Sci.
3T. A. Arias, Rev. Mod. Phys. 71, 267 (1999). Stat.Comput. 4, 244 (1983).
4M.C.Payne,M.P.Teter,D.C.Allan,T.A.Arias,andJ. 22T.L.BeckandK.A.IyerandM.P.Merrick,Int.J.Quan-
D. Joannopoulos, Rev.Mod. Phys. 64, 1045 (1992). tum Chem. 61, 341 (1997).
5G. Kresse and J. Furthmu¨ller, Phys. Rev. B, 54, 11169 23S.CostinerandS.Ta’asan, Phys.Rev.E,52,1181(1995).
(1996). 24J.MandelandS.McCormick,J.Comput.Phys.80(1989)
6D. Bai and A. Brandt, SIAMJ. Sci. Stat. Comput. 8, 109 442.
(1987). 25G. Kresse and J. Furthmu¨ller, Comput. Mat. Sci. 6, 15
7E. J. Bylaska, S. R. Kohn, S. B. Baden, A. Edelman, R. (1996).
Kawai, M. E. G. Ong, and J. H. Weare,in Proc. of the 26A.P.SeitsonenandM.J.PuskaandR.M.Nieminen,Phys.
7th SIAM Conference on Parallel Processing for Scientific Rev.B, 51, 14057 (1995).
Computation (San Francisco, 1995) p. R5471. 27L. Kleinman and D. M. Bylander, Phys. Rev. Lett., 48,
8J.-L. Fattebert, J. Comput. Phys. 149, 75 (1999). 1425 (1982).
9F. Gygi and G. Galli, Phys. Rev.B, 52, 2229 (1995). 28N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993
10N.A.ModineandG.ZumbachandE.Kaxiras,Phys.Rev. (1991).
B, 55, 10289 (1997). 29R. D. King-Smith, M. C. Payne and J. S. Lin, Phys. Rev.
11U.V.Waghmare,H.Kim,I.J.Park,M.Modine,P.Mara- B 44, 13063 (1991).
9