Direct Approach to Density Functional Theory: Heaviside−Fermi Level

J. Phys. Chem. 1996, 100, 7903-7910

7903

Direct Approach to Density Functional Theory: Heaviside-Fermi Level Operator Using a Pseudopotential Treatment Donald J. Kouri*,† and Youhong Huang‡ Department of Chemistry and Department of Physics, UniVersity of Houston, Houston, Texas 77204-5641

David K. Hoffman§ Department of Chemistry and Ames Laboratory, Iowa State UniVersity, Ames, Iowa 50011 ReceiVed: NoVember 13, 1995; In Final Form: February 29, 1996X

A new method is presented for calculating electron densities in nonperiodic, polyatomic systems using Cartesian coordinates in three dimensions. The method blends a direct approach to spin-density-functional theory and ˆ ) rather than solving Schro¨dinger eigenvalue problems, uses (1) the “Heaviside-Fermi level operator” h(EF - H (2) the distributed approximating functional for discretization and interpolation, (3) a multigrid iteration procedure for accelerating the convergence, (4) a separable, nonlocal form of pseudopotential, and (5) a fast method for solving Poisson’s equation in nonperiodic systems. Example calculations of the electronic structure for the Ne atom and the C2 and O2 dimers are presented.

I. Introduction Significant progress has been made in the last decade in computational studies of density functional theory (DFT)1-5 in fields ranging from physical chemistry to solid-state physics. DFT1,2 focuses on the density of a many-electron system as the fundamental quantity to be calculated rather than the wave function. In their molecular dynamics (MD) simulation approach, Car and Parrinello3 combined DFT and MD in a Lagrangian formalism to treat the motions of electrons and ions simultaneously. This ab initio MD scheme has been used in the simulation of liquid water.6 Teter, Payne, and Allan7 applied the conjugate gradient method to minimize the Kohn-Sham (KS) energy functional for a fixed ionic configuration directly and efficiently. This highly efficient method has been applied to a dynamics calculation of phonon modes in silicon.8 An important ingredient in these successful methods and calculations is the use of pseudopotentials9-12 and plane wave bases. The pseudopotential method, pioneered by Phillips9 and Heine and Cohen,10 allows one to eliminate the need to include atomic core states in a calculation and to replace the strong electronnucleus potential with a much weaker valence electron-core ion effective potential. The plane wave basis allows the use of an efficient FFT algorithm in 3-D Cartesian coordinates for periodic systems as well as polyatomic systems using the supercell technique5 (which has, however, certain limitations). To avoid the use of supercells and a plane wave basis, Chelikowsky, Troullier, and Saad13 have devised a finitedifference-pseudopotential (FDP) method to calculate electronic structure in real space, which is particularly suited to localized nonperiodic systems and parallel computing. Further development of the FDP method has been carried out to optimize the use of the grids,14 to improve the efficiency of the FDP * Author to whom correspondence is to be addressed. † Supported in part under a grant from the Energy Laboratory, University of Houston. Acknowledgement is made to the donors of the Petroleum Research Fund administered by the American Chemical Society, for partial support of this research. ‡ Supported in part under National Science Foundation Grants CHE9403416 and NSF Advanced Scientific Computing Fellow under Grant ASC93-10235. § The Ames Laboratory is operated for the Department of Energy by Iowa State University under Contract No. 2-7405-ENG82. X Abstract published in AdVance ACS Abstracts, April 15, 1996.

S0022-3654(95)03339-9 CCC: $12.00

scheme,15,16 and to apply it to ab initio molecular dynamics simulations.17 Recently, we have proposed a direct approach scheme20 to calculate the electron density in DFT in real space that bypasses solving the KS equation. The key ingredient of this approach is the use of the Chebychev polynomial expansion of the Heaviside-Fermi energy step operator representing the electron density function. This approach results in an iteration procedure which yields a minimization or direct “downhill” search process for obtaining the minimum energy of the electron system. The distributed approximating functional (DAF) method21-24 was used to discretize the Hamiltonian of the system in the real space. This approach was successfully applied to do calculations of the ground states of two atomic systems, He and Ne. In this paper, we generalize this direct approach to deal with spindensity-functional theory (SDFT)1,2,25 so that general polyatomic systems can be treated. The efficiency of the direct approach is further improved by a multigrid iteration technique. As a step toward enabling the calculation of electronic structure efficiently in a large system, while using 3-D Cartesian coordinates, we develop an accurate, separable nonlocal pseudopotential form which is free of spurious states. We also develop a fast method for solving Poisson’s equation for nonperiodic systems. These methods are presented in section II. Example applications illustrating our approach are presented for the Ne atom and the C2 and O2 dimers in section III. We present a brief discussion and our conclusions in section IV. II. Methods A. A Direct Approach to Calculating the Electron Density. We begin by outlining the spin-density-functional theory.1,2,25 In the Kohn-Sham formulation of the density functional theory,2 the nondegenerate ground state of a manyelectron system is determined solely by the electronic density, and this density is obtained by minimizing the energy functional

E[F] ) Ts[F] + ∫V(b) r F(b) r db r+ 1

∫∫ 2

F(b) r F(b′) r |b r - b′| r

db r db′ r + Exc[F,ξ] + ∑ R,β

ZRZβ RRβ

(1)

Here, F(r b) is the electron density and Ts[F] is the kinetic energy © 1996 American Chemical Society

7904 J. Phys. Chem., Vol. 100, No. 19, 1996

Kouri et al.

functional of this density. The function V(r b) is the static potential arising from the charged particles other than electrons, and Exc[F,ξ] is the exchange and correlation energy of an interacting electron system which is treated using a local spin density (LSD) approximation.25-27 Here, ξ denotes the spin polarization, ξ ) (F+ - F-)/F, where F( are spin-up (+) and spin-down (-) electron densities, respectively. The quantity ZRZβ/RRβ is the potential of the interacting ions (nuclei), R and β. Writing the electronic density in terms of single-electron orthonormal wave functions

and

φEF,scale ) cos-1(EF,scale) EF,scale ) Hσ,scale )

Nσ

F(b) r ) ∑∑|ψσ,i(b)| r

2

(2)

[

1 r ) - ∇2 + V(b) r + Hσψσ,i(b) 2 F(b′) r

]

r ψσ,i(b) r db′ r + Vxc,σ(b) ∫ |br - b′| r ) σ,iψσ,i(b) r

(3)

b) is the exchange correlation potential where Vxc,σ(r

Vxc,σ(b) r )

δExc[F,ξ] δFσ

(4)

In our direct approach scheme, the electron density of eq 2 is expressed as28,20

r r F(b) r ) ∑〈b|h(E F - Hσ[F,ξ])|b〉

(5)

σ

where h(EF - Hσ[F,ξ]) is the Heaviside step function operator and EF is the Fermi energy (energy of the highest occupied state). Thus, h(EF - E), which we name the Heaviside-Fermi level operator, equals one if E is less than or equal to EF and is zero otherwise. Unlike in the KS equation, eq 3, the electron density F appears on both sides of eq 5. The electron density F that satisfies the nonlinear eq 5 minimizes the KS energy functional and is thus the ground state electron density of the system. A standard approach, which is equivalent to solving the KS equation, uses the eigenstates of Hσ[F,ξ] with a given F to calculate h(EF - Hσ[F,ξ]) and to solve eq 5 self-consistently. We instead expand h(EF - Hσ[F,ξ]) directly in terms of polynomials of Hσ[F,ξ]. It will be shown that, with the appropriate scaling of Hσ[F,ξ], the polynomial expansion scheme constitutes a “downhill” minimization process. Chebychev polynomials are used to expand h(EF - Hσ[F,ξ])20,29 according to Nc

h(EF - Hσ) ) lim ∑an(EF) Tn(Hσ,scale)

(6)

Ncf∞n)0

1 a0(EF) ) - (φEF,scale - π) π 2 sin(nφEF,scale), n > 0 πn

Hσ - H h ∆H Emax + Emin 2

∆H )

Emax - Emin 2

(7)

(8)

Here, Emax is the maximum energy of Hσ. In a conventional scheme, Emin is chosen to be the minimum energy of Hσ and the Chebychev polynomial expansion of eq 6 yields a uniformly convergent solution for all energies of Hσ. However, we notice that the iteration procedure described here for h(EF - Hσ[F]) is equivalent to finding the lowest energies (eEF) of the system. This minimization can be made more efficient in this iteration scheme by a particular choice of Emin. First, EF is set to a value estimated to be larger than the Fermi energy of the system. Then, Emin is chosen to be smaller than EF but larger than the Fermi energy. In this case, the scaled Hamiltonian Hσ,scale in eq 8 possesses eigenvalues between 1 and -1 for energies larger than the Fermi energy and eigenvalues below -1 for energies smaller than the Fermi energy. Thus, the eigenvalues for energies lower than the parameter EF have a magnitude larger than one. The Chebychev polynomial expansion of eq 6 converges exponentially faster for the energies of Hσ below the Fermi energy than for the energies above EF because

Tn(x) )

{

|x| e 1 cos(ncos-1(x)), -1 -cosh(ncosh (-x)), x < -1

(9)

for the Chebychev polynomial. As the result, this iteration procedure yields an automatic minimization or downhill search for the minimum energy of the system. Most importantly, since each iteration step lowers the total energy of the system, eq 5 can be solVed by updating Hσ[F,ξ] “continuously” with a new density F after eVery one, or a small number, of iterations in eq 6. Instead of using the discretized coordinate representation basis, |r bj〉, an alternative basis is selected for iteratively calculating the density F satisfying eq 6.20 The state vector, |ζσ,i〉, defined as

|ζσ,i〉 ) h(EF - Hσ[F,ξ])|χi〉

(10)

is, in fact, an eigenstate of h(EF - Hσ[F,ξ]), with eigenvalue one, since

h(EF - Hσ[F,ξ])|ζσ,i〉 ) h(EF - Hσ[F]) h(EF Hσ[F,ξ])|χi〉 ) h(EF - Hσ[F,ξ])|χi〉 ) |ζσ,i〉

where

an(EF) ) -

∆H

H h)

σ i)1

where σ indicates the spin (+ or -) and Nσ is the number of electrons with spin σ and minimizing the energy functional subject to the orthonormality constraint, one obtains the KohnSham (KS) equation (in au)

EF - H h

(11)

Here |χi〉 is some trial wave function. The state, |ζσ,i〉, can be the null vector if |χi〉 contains no energy components below the energy EF. For a Hermitian Hamiltonian Hσ|F,ξ], a complete basis set for h(EF - Hσ[F,ξ]) can thus be generated by using eq 10 and N linearly independent trial wave functions |χi〉 with nonzero components in the spectrum of Hσ below EF. Thus,

Direct Approach to DFT

J. Phys. Chem., Vol. 100, No. 19, 1996 7905

h(EF - Hσ[F,ξ]) is a projection operator, having eigenvalues of one or zero. In order to determine the expansion coefficients in an expansion of the Heaviside step function h(EF - Hσ[F,ξ]) in terms of the |ζσ,i〉’s, the latter must be orthonormalized. Then we can write Nσ

h(EF - Hσ[F,ξ]) ) ∑|Φσ,i〉〈Φσ,i|

(12)

i)1

where the |Φσ,i〉’s are orthonormal vectors obtained from |ζσ,i〉’s, and the electron density is given by Nσ

F(b) r ) ∑∑|Φσ,i(b)| r 2

(13)

σ i)1

a separable nonlocal pseudopotential can be used to ameliorate the need for highly refined grid meshes near the nuclei. B. Distributed Approximating Functional Representation of the Kinetic Energy Operator. The distributed approximating functional (DAF)21-24 has been constructed to provide a “well-tempered” approximation everywhere to an L2 wave function, which yields comparable accuracy for the function and its derivatives as opposed to being constrained to reproduce the function exactly on a set of grid points. The particular DAF used here is constructed using “Hermite functions” (the product of Hermite polynomials times their generating function). These have the attractive feature of producing a banded matrix representation of the kinetic energy operator in a discretized coordinate representation. In particular, the Cartesian coordinate DAF representation of a function, f(x), in 1-D may be written as

In the Chebychev polynomial expansion treatment, |ζσ,i〉 is given by Nc

|ζσ,i〉 ) lim ∑an(EF)|ηn,i〉

(14)

Ncf∞n)0

by eqs 10 and 6. The dynamical basis |ηn,i〉 is calculated recursively, according to

fM(x) ) ∫-∞IM(x - x′) f(x′) dx′ ∞

where the DAF kernel, IM(x - x′), is the following approximation to the Dirac δ function

IM(x - x′) ) e-(x-x′) /2σ 2

|η0,i〉 ) |χi〉

2

M/2 (-1/4)n

1

x2πσ

∑ n)0

H2n

n!

( ) x - x′

x2σ

(17)

|η1,i〉 ) Hσ,scale|η0,i〉 l |ηn,i〉 ) 2Hσ,scale|ηn-1,i〉 - |ηn-2,i〉

(16)

(15)

Although the |ζσ,i〉 generated after one or several iterations (Nc) may not be an eigenstate of h(EF - Hσ[F,ξ]), the expectation value of the energy, (〈ζσ,i|Hσ[F,ξ]|ζσ,i〉)/〈ζσ,i|ζσ,i〉, is guaranteed to be lower than that of the initial trial wave function, (〈χi|Hσ[F,ξ]|χi〉)/〈χi|χi〉. This is because, with Hσ,scale scaled as described above, the relative components of the eigenstates with Ej < Emin have grown exponentially during the Nc applications of the Hamiltonian. Thus, the new density F(r b) can be calculated using eq 13 after orthonormalizing N of the |ζσ,i〉’s and is then used to generate a new Hσ[F,ξ]. This scheme, in practical calculations, is equivalent to solving eq 6 by varying F subject to a minimization condition. This differs greatly from the KS procedure when the minimization requires first explicitly solving the KS eigenvalue equations before using the resulting oneelectron orbitals to generate a new Hamiltonian. The efficiency of this direct approach method depends critically on the initial conditions (initial wave functions |χi〉’s and the initial density F). Since each iteration lowers the energy of the system, the closer the initial conditions are to those of the ground state of the system, the faster the iteration converges. Another factor affecting the efficiency is the singular interaction between electrons and ions. A large number of basis functions or grid points are required to describe this interaction when an electron is close to the nucleus. Thus, under these conditions, Emax becomes larger and the rate of convergence becomes slower. Also in this direct approach, a new KS Hamiltonian is generated frequently and the need to recalculate the electrostatic potential can hinder the efficiency of the method. This can be dealt with by introducing a fast method of calculating the electrostatic potential. Furthermore, the efficiency of the direct approach can be improved by a multigrid iteration method, and

Here, H2n is the even Hermite polynomial and exp[(x - x′)2/ 2σ2] is its generating function. This DAF is highly peaked about x ) x′ due to the Gaussian factor. That the error, |fM(x) - f(x)|, is controlled as a function of x is easily understood by examining the DAF in momentum space.30 Since the DAF is a convolution, its Fourier transform has a simple product form

˜f M(k) ) ˜IM(k) ˜f (k)

(18)

Here, the momentum space DAF is given by

M/2

˜IM(k) ) e-(1/2)σ k ∑ 2 2

n)0

( ) 1 2 2 σk 2

n

n!

(19)

This is a filter function, which eliminates the unimportant highmomentum components of the Fourier space function ˜f(k). In the limit M f ∞, the sum in eq 19 equals exp(σ2k2/2) and ˜IM(k) ) 1, the coordinate space DAF becomes the Dirac δ function, and eq 16 becomes an identity. For finite M and k2 < 2/σ2, ˜IM(k) is very nearly equal to 1 (this region of k2 is called the “DAF plateau”); it equals 1 exactly only for k2 ) 0. Therefore, any function, along with its derivatives, whose Fourier transform lies mostly under the DAF plateau (termed a “DAF class function”) is well approximated by eq 16. The DAF is also localized in momentum space due to the Gaussian factor. It can be shown that DAF in eq 16 is the exact identity for (M + 1)th degree polynomials. The action of the kinetic energy operator, -(1/2m) ∂2/∂x2, on the function f(x) is therefore expressed as

-

∞ 1 ∂2 f(x) ) ∫-∞KM(x - x′) f(x′) dx′ 2m ∂x2

(20)

where the DAF approximation to the kinetic energy operator is

7906 J. Phys. Chem., Vol. 100, No. 19, 1996

KM(x - x′) ) -e-(x-x′) /2σ 2

2

M/2 (-1/4)n

1

∑ 3 n)0

4mx2πσ

n!

Kouri et al.

( )

H2n+2

x-x

x2σ

(21)

On a uniform discrete grid, the discretized DAF forms of eqs 16 and 20 are

fM(x) ) ∑IM(x - xk) f(xk)∆x

(22)

k

and

-

1 ∂2 2m ∂x2

f(x) ) ∑KM(x - xk) f(xk)∆x

(23)

k

respectively, where xk is the kth point on the grid and ∆x is the (uniform) grid spacing. Thus, the kinetic energy operator in the KS Hamiltonian is represented by eq 21 in Cartesian coordinates. DAFs can also be used for interpolation (by appropriate choice of parameters) as is done in a later subsection. C. An Efficient Multigrid Iteration Procedure. The efficiency of this direct approach method depends crucially on the choice of the initial trial wave functions |χi〉. To calculate the electron density, one convenient choice of the trial wave functions is provided by hydrogen-like eigenfunctions with particular symmetries. However, in most cases, these trial wave functions still differ a lot from the ground states of the system. In order to prepare improved trial wave functions, we employ a multigrid iteration procedure. At the first grid level, the hydrogen-like eigenfunctions evaluated on a relatively coarse grid are used as |χ1i 〉. The direct approach, including the use of a pseudopotential, is efficiently carried out using this grid due to the rather small number of grid points involved. After the minimum energy is attained, the resultant wave functions, |Φ1i 〉, are DAF fit onto a finer grid, whose values are to be used as the trial wave functions for the next minimization. These trial wave functions, |χ2i 〉, are much closer to the ground states than are the |χ1i 〉. The same direct approach is repeated with this finer grid. This procedure can be applied repeatedly on even finer grids until the result for the ground state is converged to the desired precision. The transition from the coarse grid to the next finer grid is done conveniently by using the DAF method. In particular, the DAF approximation in 1-D is given by n-1 n-1 (xk )∆xn-1 χin(xjn) ) ∑IM(xjn - xn-1 k ) Φi

the strong electron-nucleus potential with a much weaker valence electron-core ion effective potential. The construction of a pseudopotential from first principles is designed to produce the following: (1) pseudo valence eigenvalues, which are the same as the true valence eigenvalues (calculated from an allelectron (AE) calculation), and pseudo valence wave functions, which are smooth (nodeless) inside the core and which agree identically with the AE valence eigenfunctions outside the core, and (2) the correct boundary condition at the surface of the core, which is achieved by matching the logarithmic derivatives of the pseudo wave functions with those of AE eigenfunctions at the corresponding energies. Hamann, Schlu¨ter, and Chiang (HSC)11 introduced a first-principles, norm-conserving pseudopotential which not only satisfies the above two requirements but also generates consistent first-energy derivatives of the logarithmic derivatives. Since these logarithmic derivatives contain all the information about the total potential within surface of the core, the core portion of HSC pseudopotential is transferable to a variety of chemical environments where the valence eigenstates are at different energies. More recently, Vanderbilt12 proposed an ultrasoft pseudopotential, which relaxes the norm conservation condition used in the HSC pseudopotential by introducing weighted, orthonormal pseudo wave functions. The resulting KS equation has the form of a generalized eigenvalue problem. The HSC pseudopotential takes the form11,36

Vps ) Vloc(r) + ∑∆Vl(r)|lm〉〈lm|

(25)

lm

where |lm〉 is the l, mth spherical harmonic. Vloc(r) is a local potential equal to -Zval/r outside the core region. Here Zval is the magnitude of the valence electronic charge. The second term of Vps is local in the radial coordinate and nonlocal in the angular coordinates. The quantity ∆Vl(r) vanishes outside the core region. Different partial waves have different pseudopotentials due to the angular momentum projection operator. Kleinman and Bylander (KB)37 discovered a simple way to replace the semilocal HSC pseudopotential by an effective, separable nonlocal form, which leads to considerable savings in applying the Hamiltonian. The KB form of the HSC pseudopotential is

Vps ) Vloc + ∑

0 0 |∆Vlφlm 〉〈φlm ∆Vl|

lm

0 0 〈φlm |∆Vl|φlm 〉

(26)

(24)

k

are points in the finer and coarser grids, where xnj and xn-1 k respectively, and ∆xn-1 is the grid spacing of the coarser grid. The form of the DAF function IM(xnj - xn-1 k ) is given by eq 17. We note that this multigrid iteration procedure differs from the full multigrid algorithm (FMG or nested iteration)31 in that it does not require “N” cycles. Methods using other MG techniques in DFT have also been reported.32-35 D. Separable Nonlocal Pseudopotential. The approximation of the core orbitals (frozen in their atomic configurations) is used. In this approximation, the core eigenfunctions are assumed to be localized around the corresponding nuclei and the valence electrons are assumed to move in a pseudopotential produced by the Coulomb potential due to the nuclei and the core electrons. Thus, the pseudopotential method allows one to eliminate the need to include atomic core states and to replace

0 〉 are the wave functions of the atom including the where |φlm HSC pseudopotential. The KB form of the HSC potential can be improved by adding more radial basis functions.38 An accurate and efficient scheme for calculating nonlocal Hamiltonian matrix elements is absolutely essential to largescale calculations. We now describe our new scheme for this purpose, constructing a separable form of the HSC potential in a discrete variable representation (DVR).39-41 The DVR is closely related to a standard Lagrange interpolation for particular choices of the mesh. It provides an approximation to a wave function, defined in an interval (a, b) by its values at the mesh points, according to

N

Ψ(x) ) ∑Ψ(xi) fi(x) i)1

(27)

Direct Approach to DFT

J. Phys. Chem., Vol. 100, No. 19, 1996 7907 the interval (0, rc), thereby obtaining the expression for fk(r)

where the interpolation functions fi(x) satisfy

fi(xj) ) δij

(28)

∫abf*i (x) fj(x) dx ) wiδij

(29)

and

with wi being a Gaussian integration weight. These two conditions are sufficient to ensure a complete set of functions {f˜i(x) (≡fi(x)/xwi), i ) 1, ‚‚‚, N} for representing discretized wave functions according to eq 27. It follows that operators are of the form

O ˆ ) ∑ 〈f˜i|O|f˜j〉|f˜i〉〈f˜j|

(30)

ij

on the mesh. In the DVR, the scalar product of two functions is given by N

∫a Φ*(x) Ψ(x) dx ) ∑wiΦ*(xi) Ψ(xi) b

(31)

i)1

For a set of N basis functions, φk(x), orthonormal in the interval (a,b)

∫abφ*(x) φ(x) dx ) δkl

k, l ) 0, ‚‚‚, N - 1 (32)

one can show that N-1

fi(x) ) wi ∑ φ* k(xi) φk(x)

(33)

k)0

For a mesh with uniform spacing, φ(x) can be taken to be sines or cosines or Fourier basis functions which are dependent on the boundary conditions.42 The mesh can also be constructed using the roots of orthogonal polynomials of degree N, so that a Gaussian quadrature results. In the DVR, the HSC pseudopotential can thus be discretized and expressed in the separable nonlocal form DVR ) Vloc + ∑ Vps

[

(

∫0N∆rfk(r) fl(r)r2 dr ) wkδkl

DVR |Ψ〉 ) Vloc|Ψ〉 + ∑∆VlxwkΨ(rk|lm)|f˜klm〉 Vps klm

(37)

The integration weight wk is r2k ∆r. Besides the reduction in the required effort to set up matrix elements, the separable nonlocal pseudopotential given by eq 34 has the advantages that (1) unlike other basis functions, it has the general form of a nonlocal pseudopotential without requiring orthogonalization, (2) the pseudo wave functions are accurately represented on the given mesh, and (3) no spurious states are produced. It is expected that the DVR procedure can be carried out in other pseudopotential constructions, such as the Vanderbilt pseudopotential, where basis functions are used to construct nonlocal KB-type pseudopotential directly. E. A Fast Method for Poisson’s Equation in Nonperiodic Systems. Pseudopotentials facilitate the optimal use of uniform grids in 3-D Cartesian coordinates for polyatomic molecules. There are well-developed rapid methods for the solution of Poisson’s equation in the 3-D Cartesian coordinates. In the supercell technique,5 the FFT method is used in a plane wave approach to solve for the static potential for periodic systems and localized nonperiodic systems with periodic boundary conditions in an enlarged space. For a given boundary condition, a multigrid approach is another rapid method for solving Poisson’s equation and has been used in localized systems.14,35 In this paper, we present a fast method using FFTs for nonperiodic systems without having to enlarge the region. The procedure is to transform the static potential so that the boundary condition of Poisson’s equation is of the Dirichlet type. In particular, Poisson’s equation is

∇2φ ) -4πF r φ(b)| r boundary ) φ0(b)

where ∆Vl(rk) ) 〈f˜k|∆Vl|f˜k〉 by eq 30. It can be shown that the nonlocal form given by eq 34 is precisely equivalent to the original form of the HSC pseudopotential given in eq 25. Thus, for a wave function |Ψ〉, 〈f˜klm|∆Vl|Ψ〉 ) (wk)1/2∆Vl(rk) Ψ(rk|lm) and 〈lm|Ψ〉 ) ∑kΨ(rk|lm)|fk〉 by eqs 31 and 27 in the DVR, so DVR that Vps |Ψ〉 gives exactly the same result as the HSC pseudopotential. That is

]

fk(rj) ) δkj

(34) ∆Vl(rk)

)

rk sin π(rk - r)/2∆r π(N - 1)(rk - r) cos Nr sin π(rk - r)/2N∆r 2N∆r π(N - 1)(rk + r) sin π(rk + r)/2∆r (36) cos 2N∆r sin π(rk + r)/2N∆r

∆Vl|f˜klm〉〈f˜klm|∆Vl

klm

(38)

where φ is the static potential and F is the electronic density. To transform the boundary condition in eq 38 to a Dirichlet boundary condition, one simply defines

r δφ(b) r ) φ(b) r - φ0(b) δF(b) r ) F(b) r +

1 2 r ∇ φ0(b) 4π

(39)

so that eq 38 becomes

) Vloc|Ψ〉 + ∑∆Vl|lm〉〈lm|Ψ〉 lm

) Vps|Ψ〉

fk(r) )

(35)

The pseudo wave functions are zero at r ) 0, and the term ∆Vl in eq 34 vanishes outside the core region r > rc; the volume element of the radial coordinate is r2 dr. We choose φk(r) ) (1/r)(2/Nr)1/2 sin(πkr/N∆r) as the orthonormal basis set in eq 33 for an even spaced mesh and zero boundary condition on

∇2δφ ) -4πδF δφ(b)| r boundary ) 0

(40)

In general, in a DFT calculation, the grid size must be large enough to enclose all of the electron density, and consequently the density vanishes at the boundary of the grid. Outside the grid, the static potential can be obtained conveniently by using

7908 J. Phys. Chem., Vol. 100, No. 19, 1996

Kouri et al.

a multipole expansion,43 i.e.,

φ>(b) r )∫

F(b′) r |b r - b′| r

TABLE 1: MGI Procedure for the Ground State of the Ne Atom

db′ r

Qlm Ylm(θ,φ)

) 4π ∑ lm 2l + 1

(41)

rl+1

where the multipole moments Qlm are given by l r db r Qlm ) ∫Y* lm(θ,φ) r F(b)

(42)

which can be calculated easily. We now want to construct a b) which is smooth inside the grid and approaches φ>(r b) at φ0(r the boundary. One possible choice for φ0 is

()

1 r r ) 4πQ00Y00 erf + φ0(b) r r0

[ ( )]

Qlm Ylm(θ, φ) r 4π∑ erf r0 lm 2l + 1 rl+1

l+3

(43)

where erf is the error function. The value of r0 is chosen so that the error function is essentially one at the boundary. The b) is finite. power (l + 3) of the error function assures that ∇2φ0(r b) is simple and analytic in this treatment. The calculation of φ0(r Since the solution δφ(r b) and the density δF(r b) are zero at the boundary, they can be expanded in terms of sine waves on the 3-D grid, according to

δφ(b) r ) ∑δˆ φmnk sin nmk

δF(b) r ) ∑δˆ Fmnk sin nmk

nπx Nx∆x nπx Nx∆x

sin

sin

mπy Ny∆y mπy Ny∆y

sin

sin

kπz Nz∆z kπz Nz∆z

(44)

Substituting eq 44 into Poisson’s equation, eq 40, one obtains the solution for δˆ φmnk

δˆ Fmnk δˆ φmnk ) 4π 2 px (n) + p2y (m) + pz2(k)

(45)

∆r (au)

Nca

total iterations

energy (au)

1 2 3 4 5 No MGI

4.0 × 10 2.0 × 10-2 1.0 × 10-2 0.5 × 10-2 0.25 × 10-2 0.25 × 10-2

10(5) 20(5) 30(5) 40(5) 50(5) 50(40)

50 150 300 500 750 2000

-119.7 -124.7 -126.4 -127.0 -127.2 -127.2

-2

a N is the number of iterations in eq 6. The number in the c parentheses is the number of updates of the KS Hamiltonian using the latest electron density.

Table 1. Without the MGI procedure, 2000 iterations are required for convergence using a grid step size ∆r ) 2.5 × 10-3 as calculated previously.20 As shown in Table 1, with the MGI procedure, the number of iterations at the finest grid level, ∆r ) 2.5 × 10-3, is 250. Furthermore, although the total number of iterations is 750 in the MGI procedure, the coarser grids have many fewer mesh points, thereby requiring less computational effort. As a result, the 500 coarse grid iterations involve a computational effort which is equivalent to 150 iterations on the finest grid. Thus, the efficiency of the direct approach is improved by a factor of 5 by using the MGI procedure. For a polyatomic system, the fineness of the grid required is reduced by using a pseudopotential. To test further the methods presented in this paper, we consider the C2 and O2 diatoms. The first-row homonuclear diatoms are well studied (both experimentally and computationally)44-46 and have become standard examples for testing new methods (e.g., in ref 13) developed in DFT. For the exchange and correlation energy functional in the LSD approximation

Exc[F,ξ] ) ∫Fxc(F,ξ) db r

(47)

and we use the Dirac LSD approximation for exchange. The LSD approximation for correlation was obtained by Vosko, Wilk, and Nusair26 by interpolating accurately the Monte Carlo results of Ceperley and Alder.47 The energy of the valence electrons is then given by48

r Fσv [12φv + Ev[Fv] ) ∑σ∑Nn vσ,n - ∑σ ∫ db σ

The expansion in sine waves can be evaluated using FFTs. Finally, the physical static potential φ(r b) for the entire system is obtained by

r φ(b) r ) δφ(b) r + φ0(b)

level

(46)

b) is where δφ(r b) is given from eq 44 by FFTs and φ0(r analytically obtained by eq 43. This method is simple and can be used for general nonperiodic systems such as polyatomic molecules, clusters, liquids, and glasses. Computationally, the FFT method allows the use of highly developed, efficient computer codes which are available in most math library subroutine packages. III. Numerical Studies First, we have tested the multigrid iteration (MGI) procedure by calculating the ground state energy of the Ne atom.20 Spherical coordinates were used in the calculation. The angular coordinate dependence is expanded using spherical harmonics, and radial DAFs are used for the radial variable discretization.23,24 Five grid sizes are used in the MGI procedure. At the nth level, the number of grid points is 2n times as many as at the 1st level (the coarsest level). The results are given in

r Fσc [xc(F,ξ) Vxc,σ(F,ξ) - xc(F,ξ)] + ∑σ ∫ db xc(Fc,0)] +

ZciZcj 1 (48) 2∑i*j Rij

where Nσv is the total number of valence electrons with spin σ, σ,n are the eigenvalues of the KS equation, eq 3, for the valence electrons, Fσv and Fσc are the densities of valence and core electrons respectively, φv is the electrostatic potential generated by valence electrons, and Zci is the net charge on the ith atom ()Zi - the number of core electrons of the ith atom). Equation 48 is rigorous if the atomic core densities, Fc, are frozen and do not overlap. The other part of the total energy, which is contributed by the core, is independent of internuclear separation and can be omitted. Using the pseudopotential method, the valence energies, σ,n, are the eigenvalues of the Hamiltonian

1 r + φv(b) r + Vxc,σ(b) r (49) Hσ[Fv,ξ] ) - ∇2 + Vps(b) 2 where Vps is the pseudopotential. In our direct approach with

Direct Approach to DFT

J. Phys. Chem., Vol. 100, No. 19, 1996 7909

TABLE 2: Spectroscopic Constants for Diatomic Molecules C2 and O2 C2

O2

re (au) Experimenta current method other methodb experimenta current method other methodb experimenta current method other methodb a

2.35 2.43 2.35 De (ev) 6.3 8.2 7.3 ωe (cm-1) 1860 1923 1880

2.28 2.39 2.27 5.2 8.5 7.6 1580 1625 1620

From ref 44. b From ref 46.

the pseudopotential, the first term in eq 48 is given by σ

Nv

r σ[Fv,ξ] h(EF - Hσ[Fv,ξ])|b〉 r ∑σ ∑n σ,n ) ∑σ ∫ dbr 〈b|H

(50)

b) is the HSC pseudopotential in the separable In this study, Vps(r b) is calculated by the fast nonlocal form, eq 34, and φv(r algorithm discussed in the previous section. The calculations have been done using 3-D Cartesian coordinates in the following range: x (-6.0, 6.0); y (-6.0, 6.0); z (-7.0, 7.0), with the internuclear axis being aligned with the z axis. Two grid spacings are used in the MGI procedure. The grid spacings ∆ are 0.4 and 0.2 au., respectively, and are taken to be the same for all the components of the 3-D Cartesian coordinates. In the coarse grid calculations, we update the Hamiltonian using the electron density obtained after every 25 (Nc in eq 6) iterations. The Hamiltonian is updated a total of five times in the course of the calculation. We use the notation “25(5)” to indicate such a combination of iterations and updates. In the fine grid portion of the MGI calculations, the Hamiltonian is updated after every 50 iterations, and the total number of updates is 5, which is therefore denoted as 50(5). It is found that 5 updates are insufficient when one does not use an MGI procedure, but rather 10 updates of the electron density are required even though all of the iterations are done with the fine grid. The DAF parameters are M/2 ) 11 and σ/∆ ) 1.65. The bandwidth of the DAF kinetic energy matrix is 31 grid points. The spectroscopic constants are obtained by fitting calculated results of several interatomic distances and using a harmonic oscillator approximation. The results for C2 and O2 are given in Table 2. The values of re and we are within 5.0% agreement with the experiment results. However, the binding energies De (the energy difference between the dimers at the equilibrium re and their constituent atoms) are substantially larger than the experimental values. This overbinding (which is typical of DFT treatments of the first-row diatomics using pseudopotentials46) will be examined in future publications, but we believe it is due either to a need to employ more grid levels in the MGI procedure or to inadequacies in the LDA and the pseudopotential used in the calculations. IV. Conclusion We have presented a combination of methods developed for efficient density functional calculations of electronic structure and properties in polyatomic systems using 3-D Cartesian coordinates. The approach emphasizes a direct strategy for

calculating the Kohn-Sham ground state electron density using a Chebychev polynomial expansion of the Heaviside-Fermi level operator, instead of solving the KS equation. With appropriate scaling of the effective Hamiltonian, this approach yields an automatic minimization procedure which is comparable to the molecular dynamics and conjugate gradient approaches. The calculations involve banded, sparse matrix-vector multiplications and orthonormalization of N vectors. Due to this orthonormalization procedure, the computational scaling of the present method is order O(N2). We are working on modifications of our approach which should reduce this to linear in N. To improve this direct approach, a multigrid iteration procedure has been developed so that a set of optimal trial wave functions, close to the true ground state of the system, can be obtained efficiently and conveniently. A separable, nonlocal form of pseudopotential has been constructed using a DVR. This separable nonlocal pseudopotential is free of spurious states and reduces the effort involved in setting up and applying the Hamiltonian matrix elements. Finally, since the direct approach requires that the effective Hamiltonian be frequently updated by new electron densities, a fast method for obtaining the electrostatic potential in nonperiodic systems has been developed. An FFT algorithm without the supercell technique is used in the method. We have illustrated our methods by calculating the ground state energy of Ne and the spectroscopic constants of C2 and O2 diatoms. Except for the overbinding, our results agree with the experiment within 5.0%. In fact, this overbinding is a wellknown difficulty observed in the calculations carried out by other groups using SDFT, and techniques for handling it have been developed.49 We shall be carrying out additional studies to elucidate this effect further. In particular, we shall use Vanderbilt’s pseudopotential for the interaction of nuclei and valence electrons and the generalized gradient approximation49 for the exchange and correlation energy functional. In conclusion, we believe the methods in this paper hold promise for many problems and will be applying them in calculations for larger polyatomic systems in the future. For such large systems, the fast DAF/convolution30 algorithm can be used to partition the system into segments according to the bandwidth of the uniform grid DAF, with each segment treated separately. This strategy is especially convenient for parallel computing architectures. Further, the localized electronic orbitals obtained from the Heaviside function can lead to a linear computational scaling scheme.50 Acknowledgment. One of the authors, Y.H., acknowledges the partial support by the National Science Foundation through a grant to the Institute for Theoretical Atomic and Molecular Physics at Harvard University and Smithsonian Astrophysical Observatory. The authors are delighted to dedicate this paper to the celebration of Jim Kinsey’s 60th birthday. Note Added in Proof: We have learned that similar polynomial expansion techniques have been used recently for calculating the Fermi distribution for a system of 64 atoms by S. Goedecker and co-workers.51,52 References and Notes (1) Hohenberg, P.; Kohn, W. Phys. ReV. 1964, 136, B864. (2) Kohn, W.; Sham, L. J. Phys. ReV. 1965, 140, A1133. (3) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (4) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford: New York, 1989. (5) See: Payne, M. C.; Teter, M. P.; Allan, D. C.; Arias, T. A.; Joannopoulos, T. D. ReV. Mod. Phys. 1992, 64, 1045 and references therein. (6) Laasonen, K.; Sprik, M.; Parrinello, M.; Car, R. J. Chem. Phys. 1993, 99, 9080.

7910 J. Phys. Chem., Vol. 100, No. 19, 1996 (7) Teter, M. P.; Payne, M. C.; Allan, D. C. Phys. ReV. 1989, B40, 12255. (8) Arias, T. A.; Payne, M. C.; Joannopoulos, T. D. Phys. ReV. 1992, B45, 1538. (9) Phillips, J. C. Phys. ReV. 1958, 112, 685. (10) Cohen, M. L.; Heine, V. Solid State Phys. 1970, 24, 37. (11) Hamann, D. R.; Schlu¨ter, M.; Chiang, C. Phys. ReV. Lett. 1979, 43, 1494. (12) Vanderbilt, D. Phys. ReV. 1990, B41, 7892. (13) Chelikowsky, J. R.; Troullier, N.; Saad, Y. Phys. ReV. Lett. 1994, 72, 1240. (14) Gygi, F.; Galli, G. Phys. ReV. 1995, B52, R2229. (15) Seitsonen, A. P.; Puska, M. J.; Mieminen, R. M. Phys. ReV. 1995, B51, 14057. (16) Hoshi, T.; Arai, M.; Fujiwara, T. Phys. ReV. 1995, B52, R5459. (17) Jing, X.; Troullier, N.; Dean, D.; Binggeli, N.; Chelikowsky, J. R.; Wu, K.; Saad, Y. Phys. ReV. 1994, B50, 12234. (18) Li, X.-P.; Nunes, R. W.; Vanderbilt, D. Phys. ReV. 1993, B47, 10891. (19) Daw, M. S. Phys. ReV. 1993, B47, 10895. (20) Huang, Y.; Kouri, D. J.; Hoffman, D. K. Chem. Phys. Lett. 1995, 243, 367. (21) Hoffman, D. K.; Nayar, N.; Sharafeddin, O. A.; Kouri, D. J. J. Phys. Chem. 1991, 95, 8299. (22) Kouri, D. J.; Hoffman, D. K. J. Phys. Chem. 1992, 96, 9631. (23) Hoffman, D. K.; Kouri, D. J. J. Phys. Chem. 1992, 96, 1179. (24) Hoffman, D. K.; Kouri, D. J. J. Phys. Chem. 1993, 97, 4984. (25) von Barth, U.; Hedin, L. J. Phys. 1972, C5, 1629. (26) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (27) Perdew, J. P.; Zunger, A. Phys. ReV. 1981, 23, 5048. (28) Golden, S. ReV. Mod. Phys. 1960, 32, 322. (29) Huang, Y.; Kouri, D. J.; Hoffman, D. K. J. Chem. Phys. 1994, 101, 10493. (30) Huang, Y.; Kouri, D. J.; Arnold, M.; Marchioro, T. L., II; Hoffman, D. K. Comput. Phys. Commun. 1994, 80, 1.

Kouri et al. (31) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran, 2nd ed.; Cambridge University Press: Cambridge, 1992. (32) Bernholc, J.; Yi, J.-Y.; Sullivan, D. J. Faraday Discuss. 1991, 92, 217. (33) Davstad, K. J. Comput. Phys. 1992, 99, 33. (34) White, S. R.; Wilkins, J. W.; Teter, M. P. Phys. ReV. 1989, B39, 5819. (35) Iyer, K. A.; Merrick, M. P.; Beck, T. L. J. Chem. Phys. 1995, 103, 227. (36) Bachelet, G. B.; Hamann, D. R.; Schlu¨ter, M. Phys. ReV. 1982, B26, 4199. (37) Kleiman, L.; Bylander, D. M. Phys. ReV. Lett. 1982, 48, 1425. (38) Blo¨chl, P. E. Phys. ReV. 1990, B41, 5414. (39) Lill , J. V.; Parker, G. A.; Light, J. C. Chem. Phys. Lett. 1982, 89, 483. (40) Light, J. C.; Hamilton, I. P.; Lill, J. V. J. Chem. Phys. 1985, 82, 1400. (41) Baye, D.; Heenen, P.-H. J. Phys. 1986, A19, 2041. (42) Colbert, D. T.; Miller, W. H. J. Chem. Phys. 1992, 96, 1982. (43) Jackson, J. D. Classical Electrodynamics, 2nd ed.; John Wiley & Sons: New York, 1975. (44) Herzberg, G. Molecular Spectra and Molecular Structure, I. Spectra of Diatomic Molecules, 2nd ed.; Van Nostrand: New York, 1950. (45) Painter, G. S.; Averill, F. W. Phys. ReV. 1982, B26, 1781. (46) Becke, A. D. Phys. ReV. 1986, A33, 2786. (47) Ceperley, D. M.; Alder, B. J. Phys. ReV. Lett. 1980, 45, 566. (48) Gunnarsson, G.; Harris, J.; Jones, R. O. Phys. ReV. 1977, B15, 3027. (49) Kutzler, F. W.; Painter, G. S. Phys. ReV. 1992, B45, 3236. (50) Mauri, F.; Galli, G. Phys. ReV. 1994, B50, 4316. (51) Goedecker, S.; Colombo, L. Phys. ReV. Lett. 1994, 73, 122. (52) Goedecker, S.; Teter, M. Phys. ReV. B 1995, 51, 9455.

JP953339U