Real space approach to electronic-structure calculations

fI) Pergamon

Solid State Communications, Vol. 94, No. I, pp. 5-8, 1995 Elsevier Science Ltd Printed in Great Britain 0038-1098/95 $9.50+.00 0038-1098(95)00006-2

REAL SPACE APPROACH TO ELECTRONIC-STRUCTURE CALCULATIONS Eiji Tsuchida and Masaru Tsukada Department of Physics, Faculty of Science, Cniversitv of Tokyo, Hongo 7-3-1, Bunkyo-ku, Tokyo 113. Japan

(Received 13 October 1994; accepted in revisedform 17 December 1994 by H. Kamimura) We have applied the Finite Element Method to the self-consistent electronicstructure calculations of molecules and solids for the first time. In this approach all the calculations are .performed in "real space" and the use of non-uniform mesh is made possible, thus enabling us to deal with localized systems with ease. To illustrate the utility of this method, we perform an all-electron calculation of hvdrogen molecule in a supercell with LDA approximation. Our method is also applicable to mesoscopic systems. Keyword: D. electronic states

I. Introduction

variables. Then the wavefunctions are interpolated in each clement as follows.

Electronic-structure calculations of materials have a long history and various methods have been developed for that purpose. One of the most popular methods for electronic-structure calculations of condensed matter systems today is the use of plane-wave basis sets and pseudo-potentials (PWPP) with local density approximation (LDA) [1]. This method usually has enough accuracy for quantitative discussion, and reasonable efficiency to treat large systems. However, for localized systems such as the first row elements or transition metals, the use of PWPP leads to a large number of plane-waves, thus making the calculation almost impossible. One way to avoid this difficulty is to use the ultrasoft pseudo-potentials [2]. In this method, the constraint of norm-conservation is relaxed, and the number of plane-waves necessary for the calculation is considerably reduced. However, the deficit charge must be restored at every step, and the calculation becomes complicated. Furthermore, for systems with large vacuum regions such as molecules or slab model of surfaces in supercell geometry, even the ultrasoft pseudo-potential does not significantly reduce the number of plane-waves. So an efficient method to deal with such systems is desirable.

8

'lj;(l: ,y,Z)lelement =

L?/Js

'W"

(1)

8=1

(2) (3) (4)

(5)

(6) (7) (8)

(9) where hI, h 2 , h3 denote the length of the three sides of the element andVJI, ... , 1/J8 denote the value of the wavefunction on each vertex of the element (Fig.l ). Note that :,us equals 1 at the s-th vertex and 0 at the other vertices. In this way, the wavefunctions are approximated by continuous piecewise polynomial functions. Here we put the values of the wavefunction of the i-th band on every vertex in the unit cell in some appropriate order and write it as 1j~,.

II. Details of the Calculation Here we propose another way to solve above difficulties. In our method, all the calculations are performed in rea} space, and non-uniform meshes are used. First. the unit cell is divided into ortho-rhombic elements in real space, and the values of the wavefunctions on eacl. vertex of the elements are taken as the basic

5

ELECTRONIC-STRUCTURE CALCULATIONS

6 z

Vol. 94, No. I

where n (0 ) is the reciprocal-space representation of the density and Z; and ·fi denote the charge and the position of each nucleus in the unit cell resp ectively. Here we define a pote ntial in real space as

8

7

V (il=-~Tt L:Z; (L: ~ eXP (io . (r-fm). ~ cell ; G

6

01'0

(14)

5

-:'>.

y

4

3

In practice, the above sum is calculate d with Ewald 's method [3]. Using this potential, the pot ential energy can be written as

.

.••.••............ E pot

x

Fig.I.

+ E pat + Ehartree + En + E ewald

1

celt

V (T)n (T)

E pat =

(10)

sr.

(15)

This can again be calculat ed in each element , and the result is the linear form of the density. So we can write with the use of a constant vector il,

An element and its vertices.

T he density corresponding to the wavefunction is approxima ted in the same way and rep resented by ii . Here t he values of the density is calculate d from the square sum of the wavefunctions. In order to conserve t he total charge, re-normalizat ion of the density is performed. T hough the present treat ment of the density is not sat isfacto ry, it causes no serious tro ubles such as inst abilities. so we adopt it for the ti me being. In this approach, the tot al energy of the syst em is represented in th e following way. First , we confine ourselves to the periodic systems to handle the long range Har tr ee pote ntial efficient ly. Furthermor e, for simplicity we consider only real wavefunctions and bare coulomb potentials, though the exten sion to complex wavefunctions and separable norm-conserving pseudo-pote ntials is st raightforward. T he to tal energy per unit cell is given by E!atal = E k i n

=

'e

(16)

ii,

v

where each element of is calculated as the sum of several integrals of the form f e/emen! V (T)ws (r') dr. In t he same way, the Hartr ee energy can be written as

1

'P (T) n (T) tli',

(17)

n (6) ?= (j2 exp (i G· T) .

(18)

Ehartree =

41r

cel l

with 'P (1') =

G1'0

Note that 'P(r') dep ends on the desity in cont rast with V (T). As a result , 'P(T) must be calculated every time the density is changed. Th e direct use of the right hand side of (18) is costly, because the fast Fourier tra nsform (FFT) is necessary to obtain t he values of n (0) for all C; . So we calculate '{J (T) by minimizing I ['{J] ==

1 {- IVrpl I

cell

2

2

(n (r"') - n) · '{J(T) } di~

(19)

with the conjugate gradient met hod under the contraint, in th e LDA approximation. Rydb erg unit s are used throughaut t he paper. Th e explicit form of th e kinet ic energy is

(11) T his integral can be calculate d exactly in each element for th e interpolated expression, and the result is the quad ratic form of th e values of wavef unctions at the vert ices. So, with th e aid of a constant mat rix G , we can write "' !~ ~ (12) Ekjn = L.... 1/1; . G · tP;· We now go on to th e pot ential energy. After the divergences of the long ran ge Coulomb potentials of nuclei and electrons ar e cancelled, the potenti al energy is expressed as E pa!

= -81r

?= n~~) . (L: Z;exp (i 0 . fi) ), G1'O

'

1

'{J (r')

di = 0,

where 7! denot e the average of the density. It can be easily shown that 'P (r"') given in Eq.(18) minimizes the functional (19) under the same const raint . T his technique for calculating Hartree potential is ado pted in [4] for a different reason. Here we approximate 'P (T) in the same way as t he density, and write it as 1ft . In t his method, cj need not be so exact when the tot al energy is far from convergent. Accordingly, when the density is changed, only a few iterations arc necessary to get th e new $ . Similar idea on this procedu re is described in

[5] . Now we can calculate the Har tree energy and the result is the bilinear form of ii and 1ft . Thus we can write with a constant matrix F , E hartr ee = 4Tt! Ift ·

(13)

(20)

cell

F . ii .

(21)

T he rest of the total energy, E xc and E ewald can be easily handl ed, and will not be given here.

Vol. 94. No. I

7

ELECTRONIC-STRUCTURE CALCULATIONS ~~~ ~'---- ~ '-..'-..

.'

<,

.

<,

:

<,

......... j

....

r.. ...·····t·········· :

:

"

:

<,

.

: e-
<,

T"

~

./

_: : - 1/

:{

. . .... .

...

<,

<,

<,

<,

<,

<,

<, <, <, <,

<,

<,

'----I'--....

<, <, <, <, <,

I

'----I'--....

'----f"-.....

'----:---..... K K,

~~

!

!

.....

F ig.3. In practice, we used t he mesh twice as dense as thi s figur e. T he mesh is taken a pproximat ely logarithmic. Fi g.2. Hyd rogen mole cul e in a cubic supercell (12 a.u . on a side) . T he calc ulations were performed wit h an eight h of t he supercell. Now t hat the t otal ene rgy is expressed as t he fun ction of {'l/7;} , we go on to min imize it with t he conjugate gradient method [1] under const ra ints

J1/J,(T) 1h(T) dr= s.;

(22)

or in our notat ion , t

w-;. F . 1/;j- = 8,J'

Then the constrai nts become IPi '

F'

7,

' 1iJj

• = Vij ,

IV. Discussion

(23)

Here F is a spa rse a nd sy mmet ric ma trix whose diagonal eleme nt is prop ortion al to the volume of the element t o which the vert ex belongs. So if we apply t he conj ugat e gra dient method di rectly when we are workin g wit h non-unifo rm meshes, t he conjugate gradient vector is destroyed in the orthogonalizat ion proced ure, and the minimum of t he total energy is never rea ched . In ord er to avoid this difficulty, we t ransfo rm the wavefun crion s as (24) 1f~; = T ' 1/~i for 11 i.

t t;

lecu le in a cubic supercell as illust rated in Fig.2. An all-electron calculat ion was perfor med with non-uniform meshes shown in Fig.3. At t he boundary between eleme nts with different size, th e wavefunct ions a re const rained to be conti nuous. T he calculat ed density is shown in Fig A. The values of t he calcula ted bond lengt h a nd vibratio nal frequency are given in Table I and compared with t hose from exper iments a nd other t heor ies. T he agreement is sat isfact ory.

T he advant ages of t his method a re as follows.

(i) Since all calc ulat ions are perform ed in real space, no FFT is necessa ry in contrast wit h PWPP. (ii) We can use non-uniform mesh es. In t he te rminology of PWPP, we can effectively change th e cut off-energy locally. Accor dingly, we can expre ss localiz ed orbitals without a ny serious difficul ti es. Fur thermore, we can red uce the num ber of variables cons idera bly when we ca lculate the elect ronic-struct ur e of surfaces in a supercell configura tion as illustra ted in Fig .5. Simila r idea is introd uced in Ref [4, 7], in which plan e-wave is used wit h adati ve Riem anni an metric. (iii) In th is meth od ,

(25)

wit h

(26) We choose t he matrix T t hat makes F' nearl y t he uni t ma trix . Anot he r possib le soluti on is to use t he unconstrained minimization [6], in which no explicit orthogo-nalizat ion is requi red . Not e t hat t his transformat ion is di fferent from t he usual pre-condition ing, which is det ermi ned by t he for m of the t otal energy.

0 . 40 0. 3 5 0 .30 0. 2 5 0 . 20 0 . 15 0.10 O. 0 5

III. Example To demonst rate t he util it y of th e pre sent met hod , we calculat ed t he equilibrium st ruc t ure of hyd rogen mo-

Fig.4. The calculat ed den sit y of hydrogen molec ule. T he sing ula rity at t he nucleus is well described.

ELECTRONIC-STRUCTURE CALCULATIONS

8

Vol. 94, No. I

Table 1. The other theory is from Ref [10], in which LDA calculations are performed with Gaussian orbitals. The same form for t xc is adopted in these calculations [11, 12].

This work Other theory Experiment

Bond Length (a.u.) 1.46 1.45 1.40

Vibrational Frequency (em-I) 4424 4277 4400

(iv) Orbitals localized in given regions of space [8) and their overlap matrix are naturally treated in this approach. (v) Our method is much easier to implement than PWPP. One clear limitation of this method is that it will be ineffective if we perform dynamical simulations with non-uniform meshes, which lead to the Pulay-like forces

[9]. Vacuum

V. Conclusion

atom

Surface

atom

Bulk

Fig.5. In the vacuum region, only a small number of elements will suffice. the effect of the potential of the nucleus appears only in an integral form, so the divergence of the bare coulomb potential at the origin does not cause any difficulty.

In conclusion, we have applied the Finite Element Method to the self-consistent calculation of electrons successfully for the first time. Since this method has several advantages over PWPP, it will be considerably favorable for the calculations of some systems. Extensions of this method to higher order interpolation formulas and separable pseudo-potentials are under way. Acknowledgement - One of the authors (E.Tsuchida) would like to thank Katsuyoshi Kobayashi and Jun Yamauchi for fruitful discussions, and Shigenobu Kimura for providing useful subroutines. The numerical calculations were performed on HP-9000.

References [1] For a review, see M.C.Payne et ol, Rev.Mod.Phys. 64, 1045 (1992) and references therein. [2] David Vanderbilt, Phys.Rev.B 41, 7892 (1990) [3] See, e.g., J.M.Ziman, Principles of the Theory of Solids (Cambridge University Press, 1972)

[7] F.Gygi, Phys.Rev.B 48, 11692 (1993) [8] G.Galli and M.Parrinello, Phys.Rev.Lett. 69, 3547 (1992) [9] P.Pulay, Mol.Phys. 17, 197 (1969)

[4] A.Devenyi et ai, Phys.Rev.B 49, 13373 (1994)

[10] G.S.Painter and F.W.Averill, Phys.Rev.B 26, 1781 (1982)

[5] R.Car and M.Parrinello, Solid State Com. 62, 403 (1987)

[11] D.M.Ceperley and B.J.Alder, Phys.Rev.Lett. 45, 566 (1980)

[6] F.Mauri, G.Galli, and R.Car, Phys.Rev.B 47,9973 (1993)

[12] S.H.Vosko, L.Wilk and M.Nusair, Can.J.Phys. 58, 1200 (1980)