Efficient and Accurate Linear Algebraic Methods for Large-scale Electronic Structure Calculations with Non-orthogonal Atomic Orbitals

1.2 MB·English

by  H. Teng

Efficient and Accurate Linear Algebraic Methods for Large-scale Electronic Structure
Calculations with Non-orthogonal Atomic Orbitals
H. Teng1,∗ T. Fujiwara1,2,† T. Hoshi3,2, T. Sogabe4,2, S.-L. Zhang5,2, and S. Yamamoto6,2
(1) Center for Research and Development of Higher Education,
The University of Tokyo, Bunkyo-ku, Tokyo, 113-8656, Japan
(2) Core Research for Evolutional Science and Technology,
Japan Science and Technology Agency (CREST-JST), Japan
(3) Department of Applied Mathematics and Physics, Tottori University, Tottori 680-8550, Japan
1 (4) School of Information Science and Technology,
1 Aichi Prefecture University, Nagakute-cho, Aichi 480-1198, Japan
0 (5) Department of Computational Science and Engineering,
2 Nagoya University, Chikusa-ku, Nagoya 464-8603, Japan and

n (6) School of Computer Science, Tokyo University of Technology,
a Katakura-machi, Hachioji, Tokyo 192-0982, Japan
J (Dated: January 26, 2011)

5
The need for large-scale electronic structure calculations arises recently in the field of material
2
physicsandefficientandaccuratealgebraicmethodsforlargesimultaneouslinearequationsbecome

greatly important. We investigate the generalized shifted conjugate orthogonal conjugate gradient
]
l method,the generalized Lanczos method and the generalized Arnoldimethod. They are thesolver
l
a methods of large simultaneous linear equations of one-electron Schr¨odinger equation and maps the
h whole Hilbert space to a small subspace called the Krylov subspace. These methods are applied
- to systems of fcc Au with the NRL tight-binding Hamiltonian (Phys. Rev. B 63, 195101 (2001)).
s
We compare results by these methods and the exact calculation and show them equally accurate.
e
m The system size dependence of the CPU time is also discussed. The generalized Lanczos method
and the generalized Arnoldi method are the most suitable for the large-scale molecular dynamics
.
t simulations from theview point of CPU time and memory size.
a
m
PACSnumbers: 71.15.-m,02.70.Ns,71.15.Pd
-
d
n I. INTRODUCTION solving simultaneous linear equations
o
c
[ (z−H)x=b, (1)
Inrecentyears,moleculardynamics (MD) simulations

1 withelectronicstructurecalculationsinnano-scalestruc-
where H is self-adjoint or real symmetric matrix, b is
v tures have attracted much attention. One needs a large
a given vector, z = ε + iη, ε is an energy parameter
8 size of systems of severalhundred thousands atoms with
6 afewhundredpico-seconds(ormorelongertime)process andη isaninfinitesimallysmallpositivenumber,respec-
7 tively. Solutions of Eq.(1) relate to the standard eigen-
in order to investigate characteristics of nano-scale sys-
4 valueproblem(ε−H)x=0. We developedthe subspace
temssuchasphenomenaofcompetitionbetweendifferent
1. physical principles or phenomena of the multi-physics, diagonalizationmethodandtheshiftedconjugateorthog-
0 onalconjugate gradient(COCG) method.13–16 Then the
e.g. energy competition between the strain field and
1 chemical bonds.1–4 Several requirements for large-scale methods were applied to the fracture propagation and
1 surface formation in Si crystals with the tight-binding
MD simulationwith electronic structure calculationsare
v: contradictory to each other, e.g. total energy accuracy Hamiltonian basedon an orthogonalbasis set.1,2 On the
Xi vs. largersystemsizeorlongerphysicaltimeofprocesses. other hand, since its Hamiltonian is described by the
tight-bindingHamiltonianbasedonanon-orthogonalba-
r There are several approaches for large-scale MD sim- sisset,theproblemofthe formationofAumultishellhe-
a ulations;5 (a) the Fermi operator expansion,6 (b) the lical nanowires was solved by the exact diagonalization
divide-and-conquer method,7 and (c) the minimization method.3,4
method (the density matrix minimization8 or the wave-
Development of efficient linear algebraic methods has
function minimization9). Another classification may be
been, so far, mainly based on the orthogonal basis
the one according as the basis set of wavefunctions; (a) sets.13,14,17–19 However,localizedbasiswavefunctionsare
the plane wavebasis set and switching between the real-
generallynon-orthogonalanditismuchdesirabletogen-
space and k-spacerepresentation,10 and (b) localizedor-
eralizethemethods tothecaseofanon-orthogonalbasis
bitals11 or tight-binding basis set.12 Computation with
set. The most popular strategy of the generalized eigen-
“massively parallel machine” is also an important issue.
value problem (represented by the non-orthogonal basis
Animportantaspectisdevelopmentofnovelalgebraic set) would be the transformation to the standard eigen-
algorithm for extra-large scale systems. The most gen- value problem.19 Our target in the present paper is to
eral and important algorithm may be the linear algebra solve simultaneous linear equations with self-adjoint or
2
real symmetric matrix S; and the orthogonality relation is expressed as
(zS−H)x=b, (2) ∗
hψ |ψ i = w(α) w(β)S (6)
α β i j ij
X
which relates to the generalized eigen-value problem ij
(εS −H)x = 0. We will investigate efficient methods ≡ (w(α),w(β)) =δ , (7)
S αβ
of solving Eq.(2) with a complex energy variable z when
the matrix size of H and S is huge. Severalalgebraic al- where w(α) = (w(α),w(α),······)t. We call the repre-
gorithmswillbediscussedanddirectlyappliedtoatight- 1 2
sentation(w(α),w(β)) the“S-product”andtherelation
bindingHamiltonianbasedonnon-orthogonalatomicor- S
Eq. (7) the “S-orthogonalization”of basis vectors w(α).
bitals in large-scale systems.
When ψ (r) satisfies the Schr¨odinger equation
The structure of the present paper is as follows. In α
Section II, the idea of non-orthogonal atomic orbitals
and physical properties (e.g. the band energy, the lo- Hˆψα(r)=εαψα(r), (8)
cal/partial density of states, numbers of occupied elec-
tron states, the chemical potential et al.) are summa- coefficients {w(α)} should be elements of an eigen-vector
i
rized. Sections III, IV and V explain three different al- ofasimultaneouslinearequationintheφ-representation;
gorithms of large-scale linear equations, i.e. the gen-
eralized shifted conjugate orthogonal conjugate gradi- H w(α) =ε S w(α) (9)
entmethod(GsCOCG),thegeneralizedLanczosmethod ji i α ji i
X X
i i
and the generalized Arnoldi method which generate the
Krylov subspace from the whole Hilbert space. In these or, in matrix-vector form,
sections, numerical examples are presented by using the
NRL tight-binding Hamiltonian. The generalized Lanc- Hw(α) =ε Sw(α). (10)
α
zos method becomes applicable to actual large systems
with a high accuracy if one use the modified Gram-
Matrices H = (H ) and S = (S ) are self-adjoint in
ij ij
Schmidtreorthogonalizationtomaintaintheorthogonal-
φ-representation.
ityofgeneratedbasisvectors. InSectionVI,wecompare
theCPU-timesofeachalgorithmanddiscusstheapplica-
bility to large-scale electronic structure calculations and
MD simulations. Section VII is conclusions. The exam- B. Green’s function and local/partial density of
pleswithoutreorthogonalizationinthegeneralizedLanc- states represented by the non-orthogonal basis set
zosmethodareshownanddiscussedinAppendixA. Ap-
pendix B is devoted to discuss the consistency between The Green’s operator Gˆ is defined as
the total energy and force.
Gˆ(z)={(ε+iη)ˆ1−Hˆ}−1, (11)
II. THEORETICAL BACKGROUND whereˆ1istheidentityoperatorandz =ε+iη. Elements
of the Green’s function matrix can, then, be written as
A. Non-orthogonal basis set and
S-orthogonalization G (z) = hφ |Gˆ(z)|φ i={S(zS−H)−1S} (12)
ij i j ij
1
We define two sets of wavefunctions, {φ (r)} and = S wα∗ wα S . (13)
{ψα(r)}, where {φi(r)} is the non-orthogonali (normal- Xk,l iknXα k z−εα l o lj
ized) basis set (e.g. atomic orbitals and ‘i’ denotes an
atomic site and energy level), and {ψα(r)} is the or- The local (partial) densities of states is expressed in the
thonormalized basis set. Then the overlap matrix Sij φ-representation as follows:
and the Hamiltonian matrix H are defined as
ij
1
D (ε) = − Im(G(z)S−1) . (14)
S = hφ |φ i= φ∗φ dr, S =1 , (3) ij π ij
ij i j Z i j ii
The normalization of the Green’s functions and the lo-
H = hφ |Hˆ|φ i= φ∗Hˆφ dr , (4)
ij i j Z i j cal/partial density of states (DOS) is then
where Hˆ is the Hamiltonian operator. The orthonormal 1 ∞
− dεImG (z)=hφ |φ i=S , (15)
basis set {ψα(r)} can be expanded in terms of {φj} as (cid:16) π(cid:17)Z−∞ ij i j ij
∞
ψα(r) = φi(r)wi(α) (5) Z dεDii(ε)=1. (16)
X −∞
i
3
C. Total band energy and Green’s function Moreover, any physical property can be expressed by
using the density matrix as
1. Density matrix and energy density matrix
2
hXi = − dεf(ε) X Im(S−1G(z)S−1)
In the simulation process, the density matrix ρij and (cid:16) π(cid:17)Z Xij ij ji
the energy density matrix π appear repeatedly in the
ij = 2 X ρ . (29)
ij ji
calculation of the Mulliken charge, the total energy and
X
ij
forces,20 whose definition may be
TheexpressionsEqs.(27)and(28)andalso(29)aresat-
1 wi(α)∗wj(α) isfied not only in the whole Hilbert space but also in
ρ = − Im dε f(ε ) (17)
ij (cid:16) π(cid:17) Z X α z−εα themappedsubspaceinwhichweconstructapproximate
α
eigen-states.
= f(εα)wi(α)∗wj(α), (18) Now we have obtained three different expressions
X
α Eqs.(21)∼(23) for N and Eqs.(26)∼(28) for E .
tot tot
1 w(α)∗w(α) Theseexpressionsnormallygivedifferentvalues,because
πij = (cid:16)− π(cid:17)ImZ dεεXα f(εα) iz−εαj (19) awpepurosuxaimllyatueseeigfiennit-setavateluseisnotfhtehemeanpepregdysiuntbesrpvaacl,e.ηFaonrd-
= f(εα)εαwi(α)∗wj(α), (20) tunately, if the formula ijρijHji = ijSijπji is satis-
X fied, the consistency betPween the totaPl band energy and
α
the force can be kept as shown in Appendix B.
where f(ε) is the Fermi-Dirac function f(ε) = {1 +
exp((ε−µ)/τ)}−1, where µ and τ are the chemical po-
tential and temperature. III. GENERALIZED SHIFTED COCG METHOD
WedevelopedtheshiftedCOCGmethodforlarge-scale
2. Physical property
linearequations(1).14,15,21 Itwasshownthattheconver-
gencebehaviorcanbemonitoredbyobservingthebehav-
The chemicalpotential µ shouldbe determinedby the ior of the “residual norm”. The shifted COCG method
equation for the total electron number Ntot: is generalized for Eq.(2) in this section.
N = 2 dεf(ε)D (ε) (21)
tot Z ii
Xi A. Definition of the problem
2 w(α)∗w(α)
j i
= Xijα Sij(cid:16)− π(cid:17)ImZ dεf(εα) z−εα (22) equTahteioneigiseneq-vuaivluaelenptrotobltehme socfatstteartiniognparryoblSecmh;r¨odinger
= 2 Sijρji, (23) (zˆ1−Hˆ)ψ(r)=χ(r), (30)
X
ij
wherez =ε+iηandεisanenergyparameterofincident
where a factor “2” is the spin degeneracy.
waves. The wavefunction ψ(r) is expanded by the set of
The total band energy of the system is given as
non-orthogonalatomic orbitals {φ };
j
occ
Etot = 2 εα =2 εαf(εα) (24) ψ(r)= φj(r)xj(z). (31)
Xα Xα Xj
2
= − Im dεεf(ε)(G(z)S−1) , (25)
π Z ii Substituting Eq. (31) into Eq. (30), one obtains general-
X
i ized linear equations
wherethesummation occrunsovertheoccupiedstates.
α (zS−H)x(z)=b, (32)
This equation can bePexpressed by the density of states,
the density matrix or the energy density matrix as
where the j-th component of the vector b is b =hφ |χi.
j j
The solution of the linear equation x(z) is then
E = 2 dεεf(ε)D (ε) (26)
tot Z ii
Xi x(z)=(zS−H)−1b=S−1G(z)S−1b (33)
= 2 ρ H (27)
ij ji with a help of Eq. (12). By setting a vector b as
X
ij
= 2 Sijπji. (28) b=ej =(0,0,··· ,0, 1 ,0,0···)t, (34)
X
ij j
|{z}
4
we can get the corresponding solution xj(z) as whereB† is the Hermitianconjugatematrix ofB andb∗
is the complex conjugate vector of b. The actual algo-
xj(z)=S−1G(z)S−1e . (35)
j rithms may be as follows. Under the initial conditions
The product of e and xj (not S-product), the i-th el-
i x =p =0, (45)
ement of a vector xj, is identical just to the energy- 0 −1
component of the density matrix ρij(ε); r0 =b, (46)
α =1, β =0, (47)
1 −1 −1
ρ (ε) = − Im et·xj(ε+iη) , (36)
ij π h i i and a definition r′ = S−1r , we evaluate the follow-
0 0
which relates to the local DOS as ing equations for the ‘seed’ energy σs iteratively for
n=0,1,2,···:
D (ε) = S ρ (ε). (37)
ii ik ki
Xk pn = rn′ +βn−1pn−1,
Then the density matrix and the energy density matrix α = (rn′,rn′)S ,
are given by the integrations of ρij(ε) as n (pn,S−1(A+σsS)pn)S
x = x +α p , (48)
n+1 n n n
ρ = dεf(ε)ρ (ε), (38)
ij Z ij rn+1 = rn−αn(A+σsS)pn,
r′ = S−1r ,
π = dεεf(ε)ρ (ε). (39) n+1 n+1
ij Z ij (r′ ,r′ )
β = n+1 n+1 S.
It should be noticed here that there is no quantities of n (rn′,rn′)S
eigen-energiesintheKrylovsubspaceand,weshoulduse
Importantpointhereisouruseofr′ =S−1r . Inactual
the calculation procedure throughρij(ε) rather than the procedure,we employ a form r =nSr′ at eanch iteration
calculation of Eqs.(18) and (20). Their resultant values n n
step by CG method. Since the overlap matrix S is real
depend on the interval of energy mesh-points for the en-
symmetric positive definite and sparse, the convergence
ergy integration and a fictitious finite value of η.
of CG iteration can be fast.
The basic theorem of the Krylov subspace is the in-
B. Generalized shifted conjugate orthogonal variance of the subspace under an energy shift σ. The
conjugate gradient (GsCOCG) method other very basic theorem is the collinear residual23
1
For non-orthogonal basis set, we can generalize the rσ = r . (49)
n πσ n
shifted COCG procedure, named the generalized shifted n
COCG (GsCOCG) method.22 The linear equations of
Owing to these theorems, once we solve the set of equa-
the‘seed’energyσ andthe‘shift’energyσ,respectively,
s tions for the ‘seed’ energy σ , we can obtain the results
s
are written down as
foranyshiftenergyσ onlybyscalarmultiplications. The
(S−1A+σ 1)x=S−1b, (40) recurrence equations for shift energies are given (all the
s
quantities are denoted by the superscript σ), with initial
(S−1A+σ1)x(i) =S−1b, (41) values πσ =πσ =1, as follows;
−1 0
where the matrix A is defined as
β
πσ = 1+α (σ−σ ) πσ+ n−1α (πσ−πσ )
A=zrefS−H (42) n+1 n n s o n αn−1 n n n−1
(50)
with an arbitrary reference energy z = ε +iη, 1 is
ref ref
the unit matrix and b = e . The seed energy and the
j and
‘shiftenergy’aregivenasε =ε +σ andε=ε +σ.
s ref s ref
Following the procedure of the shifted COCG xσ = xσ+ασpσ (51)
method,15,16 we try to find iterative n-th solutions x n+1 n n n
n
in the Krylov subspace defined
with
Kn(S−1A+σsI,S−1b) πσ
ασ = n α ,
=Span{S−1b,S−1AS−1b,(S−1A)2S−1b, n πσ n
n+1
··· ,(S−1A)nS−1b}. (43) πσ 2
βσ = n−1 β ,
This yields the residual vector r′ = S−1b− (S−1A + n−1 (cid:16) πnσ (cid:17) n−1
n
σS1)xn to be22 pσ = 1 r′ +βσ pσ .
n πσ n n−1 n−1
r′ ⊥K ((S−1A+σ 1)†,b∗), (44) n
n n S
5
Partial densities of states are shown in Fig. 1 for a C. Residual norm and convergence behavior
system of Au 864 atoms by NRL tight-binding Hamilto-
nian,24 in comparison with those by the exact calcula-
tions. In order to see the behavior of the peak positions
and the tail of the peaks, the figures are drawn in the
logarithmic scale. Two lines of GsCOCG and the ex-
actcalculationoverlapeachotheralmostcompletelyand
one can recognize an excellent agreement between the
two different calculations.
10
(a) GsCOCG
Exact
1
(a)
0
0.1
0.01
) 0.001
n
i 0
p
S -4
/ 10
gy (b) GsCOCG (b)
er Exact
n 1
E
/
1
(
s 0.1 FIG. 2: (Color on line) Convergence behavior of residual
e
t norms for a system of Au 256 atoms by NRL tight-binding
a
t Hamiltonian.24 Spectrumextendsbetween−0.5Ryto1.5Ry.
S 0.01
The inset in (b) show the total density of states D(ε) =
f
y o PiDii(ε),whereweuseafiniteimaginarynumberintheen-
t 0.001 ergyandtheprofileisofdensespikypeaks. (a)Residualnorm
i
ns krn(s,j) k at several energy points ε=−0.5, 0.0, 0.5, 1.0 Ry
De 10-4 for s-orbital. (b) Average residual norms Rn(j) with different
l (c) GsCOCG three seed energies (−0.5, 0.0, 0.5Ry) for s-orbital and they
a
i Exact all overlap with each other.
t
ar 1
P
0.1 TheusefulcharacteristicpropertyofGsCOCGmethod
is the capability of monitoring the norm of residual vec-
0.01 tors.14 The residual vectors for the seed and shift equa-
tions with an energy ε (with b=e and σ =ε −ε )
k j k k ref
arer(s,j) andr(k,j),respectively,andthe‘mapped’resid-
0.001 n n
ual vectors for the seed and shift equations r′(s,j) =
n
-4 S−1r(s,j) and r′(k,j) = S−1r(k,j). We usually need only
10 n n n
-1 -0.5 0 0.5 1 1.5 2 elements of the density matrix among near-sited orbital
pairsconnectedbynon-zeroelementsoftheHamiltonian
Energy ( ε ) Ry
or overlap matrices and the convergence monitoring is
necessary for these components.14 Therefore, in order to
FIG.1: (Coloronline)Partialdensityofstatesforasystemof monitor the convergence behavior, we adopt the “resid-
Au 864 atoms by NRL tight-bindingHamiltonian,24 normal- ual norm” defined as
ized to unity. (a) s-orbitals, (b) p-orbitals and (c) d-orbitals.
Comparison isforGsCOCG andtheexactcalculation, which Hij6=0
are almost identical to each other. Parameters in GsCOCG kr′(s/k,j) k2≡ |et·r′(s/k,j)|2. (52)
calculations are η =10−3Ry, τ =5η. The energy interval of n X i n
mesh-points is 10−4Ry. See Figs. 4, 5 and 7 for comparison. i
Furthermore, since the residual norm is different among
different energy points, the average quantity (“average
6
residual norm”) should be defined as14
1 Nene
R(j) ≡ kr′(k,j) k2
n N n
ene X
k
1 Nene 1
= kr′(s,j) k2 , (53)
n N |π |2
ene X k
k
where N is the number of energy points.
ene
The convergencebehaviorofthe residualnormfordif-
ferentseedenergiesisshowninFig.2a. Theconvergence
at the energy of the low DOS is very fast, because the
eigen-statecanbeconstructedbyasmallnumberofbasis
states. Theconvergenceoftheaveragednormisshownin
Fig.2b, whichconfirmsnumericallythe factthatthe av-
erageresidualnorm(andallthephysicalquantities)does
not depend sensitively on the choice of a seed energy.
D. Seed-switching technique
When one choosesa seed energyin anenergy rangeof
rapidconvergence,the spectra atmajorityenergypoints
have not been converged yet and one should restart the
calculationwithanewseedenergyasseeninFig.2a. The
mostdesirable seedenergy may be the one of the largest
(partial) DOS because the convergence at these point is
the most slowest.
However,even if one chooses a starting seed energy in
thehighestDOSregionandtheresidualnormattheseed
energy reaches the convergence criterion, it often hap-
pensthattherestillremainseveralenergypoints/regions FIG.3: (Coloronline)Residualnormattheseedenergiesfor
asystemofAu864atomsbytheNRLtight-bindingHamilto-
wheretheresidualnormhasnotbeensmallenough. For-
nian,24 in s-orbitals (a) and d-orbitals (b). (a) Seed-switches
tunatelytheshiftingenergydoesnotneedanyadditional
to0.0Ryfrom−0.5Ryat36-thstepandto−0.299Ryat376-
heavy computational task but several scalar manipula-
th step in the s-orbital case. (b) Seed-switches to 0.110 Ry
tion as Eq.(50). Because of this property of shifting en-
from−0.5Ryat28-thstepandto−0.034Ryat1002-ndstep
ergy, a choice of a seed energy σ can be arbitrary. As
s in the d-orbital case. Once the calculation using one seed
shown in Fig.2b, even if we start with an improper seed is converged and full convergence hasnot been achieved, one
energy and switch a seed, the total iteration times for shouldchoosethenextseedandcontinuethecalculation. The
desired convergence over whole energy range is not very gray lines show the residual norm by energy shift r′σ before
n
different. Theseed-switchingisveryefficienttechniqueto seed-switching.
avoid restarting the calculation from the beginning with
a new seed energy.16,25 One chooses a new seed energy
σnew andcancontinuethecalculationwithoutdiscarding
s
the information of the previous calculation with the old method practically useful and also the use of Eqs.(18)
σ byusingtheshiftproperty. Figure3showsthebehav- and (20) gives very efficient and accurate results.
s
ior of the residual norms in the seed-switching process.
A. Generalized Lanczos process
IV. GENERALIZED LANCZOS PROCESS AND
DENSITY OF STATES
First we define a matrix H as
The three-term recursive relation used in GsCOCG H≡S−1H, (54)
method leads us to the generalization of the Lanczos
method.18,19,26 As far as we know, the generalization of which is not self-adjoint but still satisfies the quasi-
the Lanczos method was presented first in Ref. 18. In Hermitian property in the S-product;
this section, we will stress that the generalized Gram-
Schmidt reorthogonalization process makes G-Lanczos (v,Hu) =(Hv,u) . (55)
S S
7
We can construct the three-term procedure of the Lanc- double summation of atomic sites and orbitals and it
zos process (n=0, 1···) as18,19 consumes a long CPU time. On the contrary, the cal-
culation of the density of states by Eq. (14) costs less
Hun = anun+bn+1un+1+bnun−1, (56) CPU time. The computational efficiency will be dis-
cussed later in Section VI. Here in this subsection, we
where
show several evaluated values, the density of states, the
a = (un,Hun) , integrated density of states as functions of energies for a
n S
system of gold 864 atoms of fcc structure described by
b2 = ((H−a )un−b un−1,(H−a )un−b un−1) ,
n+1 n n n n S the tight-binding Hamiltonian constructed by Mehl and
un+1 = {(H−a )un−b un−1}/b Papaconstantpolous.24
n n n+1
Severalevaluatedvaluesandconsistencybetweenthem
withconditionsb =0,b ≥0andthenthevectors{um}
0 n aresummarizedinTableI. Theparametersinthegener-
satisfy the S-orthogonality
alizedLanczosmethodareN =50,theconvergencecrite-
(un,um) =δ . (57) rion δ = 10−6Ry in the inner CG process of r′ = S−1r.
S nm
In GsCOCG, the imaginary small energy η = 10−3Ry,
ThisprocesswecallthegeneralizedLanczos(G-Lanczos) the total number of the energy integration mesh-points
process (method). It is well-known that the orthogonal- is 3,000, and the convergence criterion δ = 10−6Ry in
ity relation is brokenfor largern in the Lanczos method the inner CG and outer iteration procedures. The dif-
and this is also the case here. We adopt the modi- ference of the calculated total energy is of the order of
fied Gram-Schmidt reorthogonalization process in order 10−2 ∼ 10−3Ry. The scale of the band energy is 1 Ry
to keep the S-orthogonality. (See the results without and the relative error may be of 10−3. The bold num-
the modified Gram-Schmidt reorthogonalization process bers in Table I are a set of consistent values in eachcase
in Appendix A.)
We then stop the Lanczos process up to n = N and
assume um =0 (m=N+1,N+2,···). This procedure 1
constructs the Krylov subspace (a-1) LanEcxzaocst n) 0.9 (a-2)
pi 0.8
S
KN(H,b)=(b,Hb,H2b,··· ,HNb) (58) acy/ 00..67
er 0.5
and the matrix H is transformed in this subspace to a en 0.4
mveacStttroairxrstiounfngaiwntrittidhhieaagKnornaytaluolrvfaoslrumbba.sspisacue0an=deeja0c,hovneectgoernecroartrees- gy/Spin) Lanczos rbital-deg 0001...123
soprbointdaslst(oLoCrAuthmOo)n⇒;ormϕmal(irze)d=liXnjeaφrjc(orm)ubmjin.ation of ato(m59ic) y of States (1/Ener (b-1) Exact nsity of States (1/o 00000000........23456789 (b-2)
Thenormalizedeigen-statesinthegeneratedKrylovsub- sit De 0.1
space isψdαen(ro)ted=byN ϕn(r)Q(nα) = φj(r)wjα, (60) Partial Den (c-1) LaEncxzaocst Integrated 0000 1....6789 (c-2)
X X 0.5
n=0 j
0.4
which satisfies the Schr¨odinger equation 0.3
0.2
0.1
H Q(α) =ε Q(α), (61) 0
nm m α n -1 -0.5 0 0.5 1 1.5 2
Xm Energy ( ε ) Ry
where H =(un,Hum) =(un)tHum.
nm S FIG.4: (Coloronline)Partialdensityofstates(pDOS),nor-
malized to unity,and integrated density of states (IDOS)for
a system of Au 864 atoms by NRL tight-binding Hamilto-
B. Numerical test with NRL Hamiltonian for fcc nian.24 ComparisonisbetweenthoseofthegeneralizedLanc-
Au zosmethod(redsolidlines)andthoseoftheexactones(blue
solid lines). (a-1) and (a-2) :pDOS and IDOS for s-orbitals.
Chemical potential µ can be evaluated by using (b-1) and (b-2) :pDOS and IDOS for p-orbitals. (c-1) and
Eqs. (21)∼(23) in the generalized Lanczos method. Cal- (c-2) :pDOS and IDOS for d-orbitals. η=5.0×10−3Ry.
culation of the Green’s function uses Eq.(13) having a
8
TABLEI:ThegeneralizedLanczosprocessappliedtoasystemof864atomsoffccAuandcomparisonwiththatbyGsCOCG.
The Hamiltonian is the NRL tight-binding form.24 The chemical potential and the total energy are in Ry unit and N is
atom
thenumberof atoms in thesystem.
G-Lanczos GsCOCG
µby Eq.(21) µ by Eq.(22) µ byEq.(23) µ byEq.(21)
µ 0.3010 8669 0.2926 6528 0.2870 0985 0.3006 0832
Eq.(21) 5.5000 0001 5.4809 9705 5.4674 0401 Eq.(21) 5.5000 0000
N /2N Eq.(22) 5.5180 6152 5.4999 9998 5.4865 0739 —
tot atom
Eq.(23) 5.5316 2252 5.5135 2184 5.4999 9999 —
Eq.(26) -0.1436 6686 -0.1492 6754 -0.1531 8294 Eq.(26) -0.1436 4084
E /2N Eq.(27) ρ by Eq.(17) -0.1311 6499 -0.1364 7104 -0.1403 4132 ρ by Eq.(38) -0.1251 7255
tot atom
ρ by Eq.(18) -0.1316 2401 -0.1369 4155 -0.1408 2023
Eq.(28) π byEq.(19) -0.1311 6490 -0.1364 7094 -0.1403 4123 π byEq.(39) -0.1436 4084
π byEq.(20) -0.1316 2391 -0.1369 4146 -0.1408 2014
(for each equation determining the chemical potential) H, as19,27
and the combination of Eqs. (17) and (19), and that of
Eqs.(18)and(20)areconsistentpairsofdata. Thiscon- ln+1 = Hun (62)
sistency between the density matrix and energy density n
matrixiscrucialforconsistencybetweenthetotalenergy kn+1 = ln+1− um(um,ln)S (63)
X
and force (Appendix B). m=0
kn+1
InGsCOCGmethod,calculationofEqs.(38)and(39) un+1 = . (64)
needs ε-integration of ρij(ε) and the numerical integra- (kn+1,kn+1)S1/2
tioncausesacertainerrorintheintegrationofthetailof
Thisis the Arnoldiprocessandwecallit the generalized
the spectra. Once we reduce the value of η and increase
Arnoldi (G-Arnoldi) method. The generalized Arnoldi
the ε-points of the integration, this discrepancy can be
method generates the Krylov subspace
reduced. In the generalized Lanczos method, the com-
bination of Eqs. (23), (18) and (20) is the best scheme, K (H;b) = Span{b,Hb,H2b,··· ,HNb}
N+1
since we do not use the numerical energy integral.
= Span{u0,u1,u2,···uN}. (65)
Figure 4 shows the partial density of states and inte-
grated density of states for s-, p-, and d-orbitals. Those Thegeneratedvectorumcorrespondstoorthonormalized
by the generalized Lanczos method are compared with LCAO
exact results. The peak positions and the tails over-
ϕm(r)= φ (r)um (66)
lap excellently and we can conclude that the generalized j j
X
Lanczos method can be powerful and convenient tool. j
One must notice that the exact results are more smooth
as in G-Lanczos method. An eigen-function ψα(r) =
in the central energy region, since the correct density of
ϕn(r)Q(α) satisfies the Schr¨odinger equation
statesismoredenseinthisenergyregion. Inotherwords, n n
P
the exactresultsreserveentireprofilewith864×9peaks
H˜ Q(α) =ε Q(α) (67)
but, on the contrary, those by the generalized Lanczos nm m α n
X
method reserve the profiles with 51(= N +1) peaks or m
thedensityofstatesisexpressedasapolynomialfunction
where H˜ =hϕm|Hˆ|ϕni=(um)tHun and H˜ is an up-
of energy of the order 51 in the present calculation. mn
per Hessenberg matrix. We can say that this procedure
isakindofgeneralizationofthesubspacediagonalization
of the Krylov subspace developed before.13
B. Numerical test with NRL Hamiltonian for fcc
V. GENERALIZED ARNOLDI PROCESS AND Au
DENSITY OF STATES
We summarize, in Table II, several evaluated values
A. Generalized Arnoldi process and consistency between them, in comparison with the
resultsofGsCOCG.Thesystemisof864atomsoffccAu
by NRL tight-binding Hamiltonian.24 The bold numbers
WecanconstructtheKrylovsubspace,startingwitha inthetableareasetofconsistentresultsineachcalcula-
naturalbasis u0 =e , by using the Hamiltonian matrix tionofthechemicalpotential. N =50forthegeneralized
j0
9
TABLEII:ThegeneralizedArnoldiprocessappliedtoasystemof864atomsoffccAuandcomparisonwiththatbyGsCOCG.
The Hamiltonian is the NRL tight-binding form.24 The chemical potential and the total energy are in Ry unit and N is
atom
thenumberof atoms in thesystem.
G-Arnoldi GsCOCG
µby Eq.(21) µ by Eq.(22) µ byEq.(23) µ byEq.(21)
µ 0.3039 0679 0.2928 7198 0.2864 6136 0.3006 0832
Eq.(21) 5.5000 0002 5.4797 1790 5.4658 1773 Eq.(21) 5.5000 0000
N /2N Eq.(22) 5.5176 0887 5.5000 0001 5.4865 0981 —
tot atom
Eq.(23) 5.5311 6653 5.5135 1949 5.4999 9999 —
Eq.(26) -0.1423 0073 -0.1482 6563 -0.1522 3179 Eq.(26) -0.1436 4084
E /2N Eq.(27) ρ by Eq.(17) -0.1300 1076 -0.1351 4434 -0.1389 5498 ρ by Eq.(38) -0.1251 7255
tot atom
ρ by Eq.(18) -0.1304 6635 -0.1356 1105 -0.1394 2996
Eq.(28) π byEq.(19) -0.1300 1076 -0.1351 4434 -0.1389 5498 π byEq.(39) -0.1436 4084
π byEq.(20) -0.1304 6635 -0.1356 1105 -0.1394 2996
1 exactcalculation,whichone couldmake smallerwithin-
(a-1) Exact n) 0.9 (a-2)
pi 0.8 creasing a dimension N of the Krylov subspace.
S
y/ 0.7
c 0.6
a
er 0.5
en 0.4 VI. COMPARISON AMONG GSCOCG,
Spin) al-deg 000...123 G-LANCZOS MEMTEHTOHDOADND G-ARNOLDI
Energy/ (b-1) Exact (1/orbit 00 1..89 (b-2) A. Convergence
1/ s 0.7
( e
rtial Density of States (c-1) Exact egrated Density of Stat 000000000 1.........789123456 (c-2) wnNveGneoe-roTrLiymhushaaganestvhchcwecdzesuioismtmrmshhaaoeotadlwnerlxisffnrmiGimeoeest-rnuAihtedmhnerouotnfcadiosoltteleinhndefrvedoiaoermtnerKmiegnoereneaty(nrhcelgcshoo.tiegvdieee.ssbpvsδe.aeeuhnsq=bGadutsvhas1puiCel0aons−cirOttte6ieooClRrfimaGnNyttu)hiG.om+cehnIser1nCsttemohOswFiohadihCdgaluel.Gver2aieesr,,l
Pa Int 0.6 convergence region. One should use the same iteration
0.5 criteria both in the inner CG and outer procedures in
0.4
GsCOCG method. It sometimes happens that the resul-
0.3
0.2 tant DOS shows an un-physical behavior, e.g. negative
0.1 values of DOS, if one stops the iteration steps before
0
-1 -0.5 0 0.5 1 1.5 2 enough convergence in GsCOCG method. On the other
Energy ( ε ) Ry
hand,G-LanczosandG-Arnoldimethodsnevergivesuch
un-physicalDOSevenifonestopsatsmalliterationstep
FIG. 5: (Color on line) Partial density of states, normalized
because of the expression of Eq.(13). Furthermore, the
to unity,and integrated density of states. Comparison is be-
first N moments are preserved correctly in the energy
tweenthoseofthegeneralizedArnoldimethod(redsolidlines)
and those of theexact ones (bluesolid lines), for a system of spectra of G-Lanczos method. In the spectrum of the
Au864atomsbyNRLtight-bindingHamiltonian.24 (a-1)and presentmodelbytheexactcalculation,weobserveabout
(a-2) :pDOS and IDOSfor s-orbitals. (b-1) and (b-2) :pDOS fortyprominentpeaksandthenweuseN =50inthecal-
andIDOSforp-orbitals. (c-1)and(c-2):pDOSandIDOSfor culations of G-Lanczos and G-Arnoldi procedures. This
d-orbitals. η=5.0×10−3Ry. ispresumablythe reasonwhythe peak positionsandde-
tailed profiles in the spectra of G-Lanczos method show
excellent agreement with those of the exact calculation.
G-Lanczos method needs the Gram-Schmidt reorthogo-
Arnoldiprocess. DataofGsCOCGarethesameasinTa- nalization and also it is necessary to have enough con-
ble I. The calculated values of the total energies agree vergence in the inner CG procedure. G-Arnoldi method
with those by the generalized Lanczos method and the does not need such reorthogonalization procedure since
overalldifference is less than 1% as shown in Tables I. one solve the eigen-value problem in that subspace.
Figure 5 shows the partial density of states and the Incondensedmatters,the widthofthe valenceand/or
integrateddensity of states as functions ofenergies. The conduction bands W may be of the order of 1Ry. Then,
peak positions are deviated slightly from those by the when the number of atoms is N , we can estimate
atom
10
1
TABLE III: Chemical potential and the total energy (in Ry (a)
unit) of systems of 864 and 256 atoms of fcc Au24 by three Exact 864 atoms Lanczos 864 atoms
Exact 256 atoms Lanczos 256 atoms
different method. The values of ρ in the G-Lanczos and G-
Arnoldi methodsare evaluated by Eq.(18). s 0.5
e
G-Lanczos (864 atoms) (256 atoms) Stat
µ Eq.(22) 0.2926 6528 0.2875 4133 f
o
Etot/2Natom GE-Aq.r(n2o7l)di -(08.61436a4to7m10s4) -(02.51635a8to7m63s)9 nsity -00.5 0 0.5 1
e 1
µ Eq.(22) 0.2928 7198 0.2886 2149 D (b)
E /2N Eq.(27) -0.1351 4434 -0.1346 6195 d Exact 864 atoms Arnoldi 864 atoms
tot atom e Exact 256 atoms Arnoldi 256 atoms
GsCOCG (864 atoms) (256 atoms) rat
g
µ Eq.(21) 0.3006 0832 0.2947 1535 nte 0.5
E /2N Eq.(26) -0.1436 4084 -0.1433 1925 I
tot atom
0
-0.5 0 0.5 1
the separation of each energy level as of the order of Energy (Ry)
W/(9×N ). Presumably10%ofthisseparationwould
atom
beanenoughaccuracyintheenergyscale. Inourpresent FIG. 6: (Color on line) Comparison of IDOS of d-orbitals
case, with N ≈ 1,000, the convergence criterion can for systems of Au 256 and 864 atoms by NRL tight-binding
atom Hamiltonian.24 (a) G-Lanczos (solid line) and exact calcula-
be chosen as 0.1× W/(9 ×N ) ≈ 0.1× 1/9,000 ≈
atom
10−5Ry. We also observed that the maximum iteration tion (chain line) for 864 atoms (red) and 256 (blue) atoms.
(b) G-Arnoldi (solid line) and exact calculation (chain line)
steps are almost the same for the convergence criterion
for 864 atoms (red) and 256 (blue) atoms. The agreement
δ = 10−5 and δ = 10−6 in GsCOCG and G-Lanczos
betweentheresultsbyG-LanczosorG-Arnoldimethodsand
methods. This is the reasonwhy we choose δ =10−6Ry.
those bythe exact calculation is excellent.
We compare the results of the chemical potential and
the total band energy of systems of 256 and 864 atoms
in Table III. The difference of the chemical poten-
tial µ is of the order of 10−3Ry, and that of the to- methodfortheNRLHamiltonianoffccAusystemof256
tal band energy E is of the order of 10−4Ry. The and 864 atoms.24 The total number of orbitals equals to
tot
level separationin 256 atom system can be estimated as nine times of the totalnumber ofatoms (1 s, 3 p’s and 5
1.0/(9×256)≈4×10−4Ry andthatin864atomsystem d’s). We use, in the inner CG process of the generalized
10−4Ry. The difference in the chemical potential of two Lanczos (G-Lanczos) method, the convergence criterion
systems of different sizes is due to the difference of total δ = 10−6Ry. Two numbers in the row of the CPU time
numberoflevelswhichchangesthe value ofthe chemical are referred to those of N = 50 and N = 100, respec-
potentialsensitively. Ontheotherhand,thedifferenceof tively,forG-LanczosandG-Arnoldimethods,thoughthe
E is just the quantity related to the overall spectrum resultsofN =100almostcoincidewiththoseofN =50.
tot
and we can see an excellent convergenceof the results of ForGsCOCG(withshifted3,000energypoints),thedata
N =50. shown here are those of δ = 10−6Ry both for in the in-
Figure 6 shows the actual convergence behavior of ner and outer iteration processes. The repeated time of
IDOS by G-Lanczos and G-Arnoldi method. The agree- the inner CG process (S−1x part) in GsCOCG and G-
ment between the results by G-Lanczos or G-Arnoldi Lanczos method is 10 ∼ 11. (Repeated time of 25 ∼ 27
methods and those by the exact calculation is excellent is needed for δ =10−18Ry.)
both for systems of 864 atoms and that of 256 if we The system size dependence of the CPU time is lin-
adopt N = 50. The most apparent difference appears ear for the generalized Arnoldi (G-Arnoldi) method and
intheIDOScurvesofexactcalculationofsystemsof256 thegeneralizedLanczos(G-Lanczos)method,bilinearfor
atoms and 865 atoms in the mid energy region, and the GsCOCG method and cubic for the exact calculation.
calculated results by our present methods present this The generalizedArnoldi method is extremely efficient in
difference with complete fidelity. electronic structure calculations of extra-large systems
with several hundred thousands atoms.
B. CPU times
C. Applicability to large-scale electronic structure
We summarize, in Tables IV, the CPU times (by calculations and MD simulations
using single CPU of the standard workstation) for (s-
orbitals) D of the generalized Lanczos and the gener- In the exact calculation and GsCOCG method, calcu-
ii
alized Arnoldi methods with Eq.(14), that of GsCOCG lations of physical properties, such as the density ma-
with Eq.(37), and that of the exact diagonalization trix, the energy density matrix, and chemical potential,