Review of quantum Monte Carlo methods and their applications

THEO CHEM ELSEVIER

Journal of Molecular Structure (Theochem) 394 (1997) 75-85

Review of quantum Monte Carlo methods and their applications Paulo H. Acioli Depnrtamento

de Fkica,

Universidade

de Brasilia,

Brasilia,

DF, 70.910-900,

Bruzd

Received 22 November 1995; accepted 2 April 1996

Abstract Correlation energy makes a small but very important contribution to the total energy of an electronic system. Among the traditional methods used to study electronic correlation are coupled clusters (CC), configuration interaction (CI) and manybody perturbation theory (MBPT) in quantum chemistry, and density functional theory (DFT) in solid state physics. An alternative method, which has been applied successfully to systems ranging from the homogeneous electron gas, to atoms, molecules, solids and clusters is quantum Monte Carlo (QMC). In this method the Schrijdinger equation is transformed to a diffusion equation which is solved using stochastic methods. In this work we review some of the basic aspects of QMC in two of its variants, variational (VMC) and diffusion Monte Carlo (DMC). We also review some of its applications, such as the homogeneous electron gas, atoms and the inhomogeneous electron gas (jellium surface). The correlation energy obtained by Ceperley and Alder (D.M. Ceperley and B.J. Alder, Physical Review, 45 (1980) 566), as parameterized by Perdew and Zunger (J.P. Perdew and A. Zunger, Phys. Rev. 823 (1980) 5469), is one of the most used in DF’T calculations in the local density approximation (LDA). Unfortunately, the use of the LDA in inhomogeneous systems is questionable, and better approximations are desired or even necessary. We present results of the calculations performed on metallic surfaces in the jellium model which can be useful to obtain better approximations for the exchange and correlation funcrionals. We have computed the electronic density, work function, surface energy and pair correlation functions for a jellium slab at the average density of magnesium (u, = 2.66). Since there is an exact expression for the exchange and correlation functional in terms of the pair correlation functions, the knowledge of such functions near the edge of the surface may be useful to obtain exchange and correlation functionals valid for inhomogeneous systems. From the exchange and correlation functional we can conclude that the exchange-correlation hole is nearly spherical in the bulk region but elongated in the direction perpendicular to the surface as the electron approaches the edge of the surface, showing the anisotropic character of the electronic correlation near the surface. Q 1997 Elsevier Science B.V. Keywords:

Quantum Monte Carlo; Correlation

energy; Jellium model; Variational Monte Carlo; Diffusion Monte Carlo

interaction (CI) and multi-configuration

1. Introduction

tent field (MCSCF). An accurate

determination

of correlation

energies

is

a goal of both quantum chemistry and solid state physics. To go beyond the single-particle picture one often resorts to methods with more complicated wavefunctions such as linear combination of Slater determinants as used in methods like configuration

There

self-consis-

are also perturbation

tech-

such as many-body perturbation theory (MBPT) and coupled clusters (CC). In solid state physics the most popular way to treat correlation effects is through the density functional theory (Dfl) of Hohenberg and Kohn [l]. In this approach the single-particle picture is maintained while niques

0166.1280/97/$17.00 0 1997 Elsevier Science B.V. All rights reserved PZZ SO 166.1280(96)0482 1 -X

76

P.H. Acioli/Journal

of Molecular Structure (Theochem) 394 (1997) 75-85

correlation effects are introduced via a functional of the density. Although formally exact, in practice one has to make approximations for the exchangecorrelation functional such as the local density approximation (LDA) and beyond. An alternative approach, which has demonstrated potential in both quantum chemistry and solid state physics, is the use of quantum Monte Carlo (QMC) methods. From the earlier applications in nuclear physics problems to the more recent applications to atoms, molecules, solids, surfaces and clusters, QMC has demonstrated a potential for feasible calculations with chemical accuracy. In this paper we will review the treatment of electronic correlation via QMC. Our starting point is the variational Monte Carlo (VMC), the simplest method in this class. In this method the expectation values of operators with respect to a trial wavefunction are computed using the Metropolis algorithm [2]. We discuss some of the important aspects of a good trial wavefunction and its role in QMC methods. We then review a more sophisticated method, namely, diffusion Monte Carlo (DMC). In DMC one starts with the Schrodinger equation in imaginary time, which can be thought as a diffusion equation and as such can be simulated by stochastic methods. The trial wavefunction is projected onto the ground state via successive application of Green’s function of the system. We will also discuss the potential and limitations of the method and will discuss some applications in quantum chemistry and solid state physics.

2. Quantum

Monte Carlo

2.1. Variational Monte Carlo The essence of VMC is to compute the expectation values of operators as a statistical average over configurations (R = r, ...rN} sampled according to the square of the trial wavefunction i+(R)]’ using the Metropolis algorithm. For example, the expectation value of the energy is computed as

where E,(R) =HQ(R)N(R) is called the local energy. The main advantage of this method is the fact that one can use explicitly correlated wavefunctions. The emphasis now is on building a good trial wavefunction. One can associate a variance with the estimate of the variational energy, so that, in addition to minimizing the energy, a good trial wavefunction will have the role of minimizing the variance and therefore limiting the length of the Monte Carlo simulation. Because of this dual role, Umrigar et al. [3] proposed a method of optimizing the trial wavefunction by minimizing the variance of the local energy. The basis of their argument is that in the limit the trial wavefunction approaches the ground state, the variance should be zero. Of course the inclusion of all relevant terms in the trial wavefunction can be a formidable task. Fortunately, there are more accurate methods such as diffusion Monte Carlo (DMC) which project the true ground state from the given trial wavefunction. 2.2. DifSusion Monte Carlo The basis of DMC is the Schrodinger imaginary time aqr, t) - p= at

N fi2 C - 2mV2 + V(R, t) 1=l

equation

in

(2)

One can think of this equation as a diffusion equation (with diffusion constant D = h2/2rn) with an additional branching term, given by the potential energy. This diffusion process can be simulated using Monte Carlo techniques. Take the eigenfunction expansion of the projection operator associated with this equation, exp( - t(H CY~)),and apply it to a given initial state of the system I\k,J, we get the state l+(t))=

5 eXp(-t((Yi-(Yg))I~,X~;I\kO) i=l

(3)

where (Y, are the eigenvalues of H and I+,,) are the eigenfunctions. It is clear from Eq. (3) that as t - ~1, all the terms in this expansion will decay, except for the ground state. That is,

(1)

P.H. A&Ii/Journal

of Molecular Structure(Theochem) 394 (1997) 75-85

Hence, as long as the overlap between the initial state and the ground state wavefunction is non-zero ((901*a) # 0), the asymptotic state will be proportional to the ground state of the Hamiltonian. Our problem is to obtain an explicit form for the projection operator, denoted by G(R’,t’;R,t), and then simulate the following time evolution equation: @(R’? t’) =

J

dRG(R’, t’; R, t)+(R, t)

Monte Carlo methods. An expression for G(R’,t’;R,t) valid for small time steps is given by

- $V(R’)

+ V(R)

(6)

The first term in the exponential is identified as a diffusion process in 3N dimensions, with diffusion constant D. The second term is identified as a branching term controlled by the potential energy of the system. Using this expression in Eq. (5) one can project the ground state wavefunction from any trial wavefunction using Monte Carlo methods to perform the integration, But there can be problems with the branching term, since V(R) will usually have singularities and this can lead to uncontrolled branching. To avoid this problem we use the idea of importance sampling. We take the evolution equation (Eq. (5)) and multiply it on the left by Q(R’) and insert q(R)/ ‘P(R) on the right-hand side to get f(R’,t’)=

s

dRG(R’,t’R,t)f(R,t)

(7)

= \k(R)+(R,t) and G(R’, t’R, t)= Like before, we want to obtain an explicit form for G. Multiplying the Schrodinger equation in imaginary time on both sides by \k(R) we get where

f $5~)

*(RI) %!&!j%).

d! CR,f) -------_ at

- DV2f (R, t) + DVj(R, +(EL(R)-ETV'(R,~)

where F(R,t) = 2V*(R)M(R) and 9(R). This equation resembles the we had before, with an additional on the right-hand side), as we can

G(R’, t’ R, t) =

for G.

1 (4aDAt)“N/2

(3

using

x exp

time step approximation

t)F(R, t)

(8)

G(R’,t’;R,t)

can be thought of as a transition probability to R’ given the initial configuration R. The diffusion and branching process will proceed as follows. First, we sample R’ from a Gaussian centered at the drifted position R + DAtF with variance 2DAt. We then weight (or branch) the final confgurations according to exp -At ‘I (R)fi;E1-(R’)ET . This pro>I cess is repeated untilf( [ d ,t) reaches its asymptotic state f(R, t -

~0)= ‘I’@)%(R)

(10)

where (P,(R) is the exact ground state for the original Hamiltonian H of the system. In such calculations the main quantity of interest is the ground state energy, which can be evaluated as the average value of the local energy with respect to the asymptotic distribution fi =fiR,t - x), i.e. tEL)

=

s

dRwL(R)

(11)

SdRfi The last equation can be rewritten lE

>=

L

(@o(R)IH~WOE @o(R)lWR)) ’

as

(12)

where we used the fact that H is hermitian and therefore &(R)IH = (90(R)IEo. Any operator that commutes with the Hamiltonian can be evaluated in the same way. For operators that do not commute with the Hamiltonian we define a mixed estimate (13) The error associated with this estimate is linear in the difference between the ground state and the trial wavefunctions. A better estimator is a simple linear extrapolation

EL(R) = (H+(R))/

diffusion equation drift term (second see from the short

(@,,, = 2(&i,

- (@,,,

(14)

where (b),,, is the variational estimator. One can easily show that the error from this estimator is

P.H. AcioWJournal

78

of Molecular Structure (Theochem) 394 (1997) 75-85

quadratic in the difference between the trial and ground state wavefunctions. For some quantities like densities and pair correlation functions, one can use a geometric extrapolation:

(15) In the limit of very accurate wavefunctions both estimators yield the same answer. In this method there is a systematic error due to the approximation of Green’s function for a short time step, but Reynolds et al. [4], showed that if A? is adjusted so that the acceptance ratio is larger than 99%, this error will be small. This systematic error can be eliminated by sampling the exact Green’s function of the system G(R, R, At) =

I 1 + At(E - EL)

(16)

This method is called Green’s function Monte Carlo (GFMC) and it was developed by Kalos et al. [5], but will not be discussed here, as it will not be used in this work. DMC produces the exact ground state for a system composed of bosons, but one must be careful in treating fermions, because of the antisymmetrization of the wavefunction, in which case the mixed distribution f(R,t) is not always definite positive. In the next section we will discuss the problems with fermions, and the so-called fixed-node approximation. 2.3. Fermi statistics and thefied-node

approximation

For the diffusion Monte Carlo algorithm to work, the mixed distribution f(R,t) must be always positive definite (aO(R)\k(R) 2 0). For bosons this is true, as the ground state wavefunction is nodeless, and a good trial wavefunction should be constructed in this way. For fermions, because of the antisymmetric character of the wavefunction, both aO(R) and ‘P(R) will have nodes, and usually they will not be the same. That means that there will be regions in space where f(R,t) will be negative, making direct sampling out of the question. One way to avoid this problem is to restrict the random walks to regions in space where the trial and the true ground state have the same sign. This is equivalent to imposing the boundary condition that

@0(R) = 0 whenever \k(R) = 0. Thus, the asymptotic distribution (as t - a) will be the product of the trial wavefunction and the bosonic ground state for each region bounded by the nodes. This procedure, known as the fixed-node approximation, will generate an antisymmetric ground state wavefunction and the energy obtained will be an upper bound on the ground state energy. A simple proof of the variational character of this approximation is given in [4]. Thus, in the limit that the trial wavefunction nodes are the same as those of the ground state wavefunction, the procedure will project the correct ground state wavefunction and energy. It is clear now that the role of the trial wavefunction is very important. First, because a good trial wavefunction will reduce the fluctuations in energy, meaning faster convergence and smaller variance. Second, because if the nodes are good then the energy obtained will be very close to the true ground state energy. The practical implementation of this approximation is very simple: we just include the extra step to kill a configuration whenever it reaches or crosses a node (i.e. *o(R)*(R) = 0). The typical error in this approximation is about 10% of the difference between the VMC and DMC energies 161. Beyond the fixed-node approximation there are procedures known as release-node [7] or transientestimate [8] methods, which allow one to obtain a more accurate estimate of the ground state properties. The release-node method allows the walkers to cross the nodes, not killing them instantly, and letting them continue for some time. But their contributions to averages will change sign each time the walkers cross a node. In this method the antisymmetry is taking into account exactly, but the process is numerically unstable and in practice is restricted to systems in which the nodes are very accurate or to systems where the Pauli principle is relatively unimportant. 2.4. The trial wavefunction Above we briefly discussed the role of a good trial wavefunction in QMC. As in any other variational methods, an accurate wavefunction means lower energy. In QMC a good trial wavefunction has the additional role of reducing the variance of the energy, meaning shorter simulations to achieve a desired statistical error and the very important role of

P.H. AcioWJournal

of Molecular Structure (Themhem)

providing accurate nodes, in which case the error of the fixed-node approximation is greatly reduced. In this section we would like to discuss some aspects of trial wavefunctions and the specific form we chose for our calculations. We start from the simplest antisymmetric wavefunction of a system of N fermions (in our case electrons). This function is a Slater determinant ]9] constructed from single-particle orbitals

(17)

The Hartree-Fock (HF) approximation [lo] assumes that the N-electron wavefunction consists of a single Slater determinant. After the single-particle orbitals have been optimized, the typical variational energy of an atom in this approximation is about 99% of its experimental value. Although this sounds remarkably good, and in fact it is, the missing 1% is a good fraction of the relevant energies in chemistry, like binding energies and electron affinities. Very sophisticated methods have been devised to try to recover this missing energy, known as the correlation energy. To improve upon the HF approximation one should start by including conditions on the many-body trial wavefunction which are known exactly. Among them is the so-called “cusp condition”. When the distance between two particles, two electrons or an electron and the nucleus approaches zero, the kinetic energy must diverge to cancel the divergence of the potential energy. It can be shown that to cancel such a singularity the trial wavefunction \k must satisfy the following conditions:

dru

=

I

i for like spins

79

single determinant wavefunction and how the cusp condition is described in them. Although a single determinant may not be a good description of some electronic systems, one can obtain a complete set of Slater determinants that spans the space of antisymmetric wavefunctions. So, in principle, one can write the exact ground state of any electronic Hamiltonian as a linear combination of Slater determinants. Examples of methods that use such an expansion are the configuration interaction (CI) and multi-configuration self-consistent field (MCSCF) [11,12]. In principle, the exact ground state wavefunction is a linear combination of an infinite number of Slater determinants. A large expansion is needed because the expansion is trying to reproduce a function with a cusp with a basis set of functions with no cusp. Despite the fact that an accurate description of the cusp may be hard, a multi-determinant wavefunction may be necessary to describe the electronic configuration of some systems. In the language of QMC, the role of the first few determinants on a multi-determinant wavefunction is to improve the nodes of the wavefunction while the role of the rest of the expansion is to reproduce the cusp conditions and the correlation hole. In the next approach, we discuss the opposite case, where the cusp condition is exactly reproduced and no attempt is made to correct the nodes of the single determinant trial wavefunction. The second approach is to use a correlated wavefunction [ 13,141. In this approach the trial wavefunction is a product of a Slater determinant times a function describing the correlation of pair of particles:

where

i for unlike spins

dln*(r,, =0)

394 (1997) 75-85

(IS)

and (19) where rv is the distance between two electrons and r,, is the distance between an electron and the nucleus. A good trial wavefunction for an electronic system will include these cusp conditions. We would like to discuss two different approaches to improve upon the

U(r) =

sr

(21)

and a and b are adjustable parameters and det( I ). det( ] ) are Slater determinants for spin-up and spin-down electrons, respectively. The parameter b will be a measure of the correlation hole and a is adjusted to satisfy the electron-electron cusp conditions (Eq. (18)). This kind of trial wavefunction satisfies exactly the cusp condition without the necessity of a large expansion of Slater determinants. Including the correlation function significantly

80

Table 1 Comparison

P.H. AcioWJournal

of calculated

of Moleculur Structure (Theochem) 394 (1997) 78-85

and measured energies for Be

Method VMC [I91 Fixed-node QMC [ 191 CI: 650000 determinants [20] CI: Extrapolation to large number of determinants Experiment [21] Experiment: non-relativistic correction

Energy (hartrees)

[20]

improves the energy, but there is still room for improvement. In terms of QMC, the nodes of the wavefunction are not improved, although it should reduce the fluctuations in energy, thus reducing the computational time. The functional form of the correlation function (Eq. (21)) exactly describes the cusp condition, but may not exactly describe the shape of the correlation hole. One can improve on that by using a more flexible form with more adjustable parameters and can also include higher order correlation like correlation between three electrons and two electrons and the nucleus. Such forms have been suggested and greatly improve the ground state energy [3,15]. But there is still the question of how to improve the nodes. Umrigar et al. [3] have combined the advantages of both approaches for the study of few-electron systems. They considered a trial wavefunction which is a product of a linear combination of a few Slater determinants (typically four) times a correlation function which included correlations between two electrons and two electrons and the nucleus which had 44 adjustable parameters. They recovered 99% of the correlation energy for the Be atom. This indicates that one must include the information about the electronic configuration of the system as well as a good description of its correlation holes. An alternative approach to improve the nodes of a single determinant trial function is the use of backflow correlations. In this case there is a change in variables of the single orbitals, which changes the nodes of the wavefunction. By optimizing the parameters on this change of variables and the parameters describing the two and higher order correlations one can effectively improve the VMC and DMC energies of the system. The inclusion of back-flow correlations has been shown to significantly improve the ground state energy of a two-dimensional electron gas 1161.

-

14.66648(l) 14.66718(3) 14.66557 14.66737(3) 14.669331 14.667375(25)

3. Applications

of QMC

3. I. All-electron

calculations

In this section we review some of the important results of realistic electronic structure calculations using QMC. In particular, we would like to emphasize that one can obtain accurate results within the fixednode approximation. This brief review is not intended to substitute for more thorough reviews such as those by Anderson [ 171 and Hammond et al. [ 181. We would like to start with all-electron atomic calculations, and emphasize that one can obtain very accurate wavefunctions such as the ones obtained by Umringar et al. 131 for Be and Ne. In Table 1 and Table 2 we present a comparison of QMC calculations with stateof-the-art CI calculations and with experiment. One can see that VMC calculations recover 99% of the correlation energy for Be, in good comparison with the extrapolation of the CI result using 650000 determinants and 9 1% of the correlation energy of Ne. These results were further improved with the use of fixed-node (F’N) DMC, recovering 100% of the correlation energy of Be and 98% for Ne. Another good example of all-electron calculations is the calculation of the LiH molecule, where there have been several calculations using DMC in both the fixed-node approximation and using the released-node method. In Table 3, extracted from [17], we show

Table 2 Comparison Method VMC VMC DMC Exact

[I51 [3] [22] [23]

of calculated

and measured energies for Ne Energy (hartrees) -

128.8771(5) 128.884(4) 128.9274(65) 128.937

P.H. AcioWJournal

Table 3 Comparison

of Molecular Structure (Theochem) 394 (1997) 75-RS

81

where of calculated

Method

and measured energies for LiH”.’ Energy (hartrees)

FN-DMC [24] RNDMC [24] RN-GFMC [25] Exact

-

8.0680(6) 8.0700(2) 8.07021(5) 8.07023

* Using internuclear distance R = 3.015 bohr. h Extracted from [ 161.

some of these results. One can see that the two most accurate released-node calculations are those of Caffarel and Ceperley [24] and Chen and Anderson [25], which are in agreement with the estimated true ground state energy of - 8.07023 hartrees. 3.2. Pseudopotential

calculations

One of the limitations of the method, not exclusive to QMC, is the treatment of heavy atoms, as the computation time scales as .Z5.5-6.5[26,27], where Z is the atomic number. Fortunately, one can use pseudopotentials and perform valence-only calculations. The advantages of using pseudopotentials is twofold. First, the number of electrons is reduced and, second, the computation time now scales as Z:i, where Z,rr is the number of valence electrons. The simulation of a valence-only Hamiltonian with VMC is straightforward. As in usual VMC the configurations are sampled from the square of the trial wavefunction, and the energy is computed as the average of the local energy. The only complication is the need to perform an additional angular integration of the non-local part of the pseudopotential. The details of this integration are discussed in [28,29]. Although this integration slows down the calculation, simulations are still feasible and have been performed [283 1] for solids and molecules. In DMC there is a fundamental complication related to the use of non-local pseudopotentials. To understand it better, let us rewrite the diffusion equation for the valence Hamiltonian as

4f (R t> - DV2f (R, t) + DV.F(R, t) -= at

- 6 f(R, t) + 4R) (22)

>

t(R) =

V,I,,WC 0 +(R,t)

Vn,oc’W,t) -

WR,t)

(23)

and V,I,, is the non-local pseudopotential. This equation is still a diffusion equation with drifting and branching, with an added term that represents nonlocal branching. To keep the random walk local, one is forced to make the approximation of neglecting the term E(R). This is called the locality approximation [29] and it is justified whenever an accurate wavefunction is obtainable, in which case E(R) - 0. In this approximation the non-local part of the Hamiltonian is evaluated in the same fashion as in VMC. The main consequence of this approximation is the loss of the variational character of the method. The impact of such a loss is reduced if the trial wavefunction is very accurate, as the error in the approximation is quadratic in the difference between the trial and exact wavefunctions. Simulations using this approach have been successfully carried out for small molecules [29], solid N [32], the Fe atom [33] and carbon and silicon clusters [3 1,341. In all cases the wavefunction was very well optimized, so the error due to the locality approximation was smaller than the statistical error of the simulation. The third approach for valence-only simulations tries to overcome the difficulties involving the use of non-local pseudopotentials by recasting the problem in terms of a local pseudohamiltonian with a modified kinetic energy tensor [35]. The main advantage of this approach is that DMC can be applied with only minor modifications of the algorithm and the variational character of the method is preserved. The basic idea of the method is to rewrite the valence Hamiltonian as

where the al, b, and vI are local functions adjusted to reproduce the pseudo-orbitals generated from the nonlocal pseudopotential. In order to have a well-behaved Hamiltonian the functions al and b, must satisfy a(r) > -1

a(r)+b(r)

> -1

(25)

82 Table 4 Companson

P.H. Acioli/Journat

of calculated

oj Molecular Structure (Theochem) 394 (1997) 75-S

and measured cohesive energies of diamond and solid silicon and aluminum

System

Method

Carbon (diamond)

VMC, pseudopotential Experimenta

Silicon

VMC, pseudopotential [30] VMC, pseudohamiltonian [36] FN-DMC, pseudohamiltonian [36] Experiment”

4.81(7) 4.38(4) 4.51(3) 4.63(S)

Aluminum

LDA, pseudopotential [39] LDA, pseudopotential [37] FN-DMC, pseudohamiltonian Experiment [40]

3.983 4.242 3.53(3) 3.392(44)

[ 161

Cohesive energylev [30]

atoms’

7.4X7) 7.37

[37]

a See lists in [30,36].

The question is, can we always obtain a(r) and b(r) such that Eq. (25) is satisfied? Unfortunately, the answer is no for many atoms and common norm-conserving pseudopotentials. However, for many atoms that form s-p bonded systems, an exact reproduction of s and p orbitals is possible for most atoms, with the exception of rare earths and transition metals. This pseudohamiltonian formulation has been applied in conjunction with QMC to small atoms and dimers [35] and for the study of solid Si [36] and solid aluminum [37]. There have also been applications of the formulation in the framework of DFT [38]. In Table 4 we present the results of QMC calculations using non-local pseudopotentials and the local pseudohamiltonian formulation for solid Si [30,36], C (diamond) [30] and Al [37]. One can see that both approaches lead to accurate determination of electronic properties such as the cohesive energy. As one can see, QMC simulation of heavy atoms is possible if used in conjunction with a non-local pseudopotential or a local pseudohamiltonian. The question now is whether the pseudopotentials generated within the framework of single-particle theories are accurate enough to be used in QMC simulations. An attempt to generate pseudopotentials from correlated wavefunctions based on the properties of the one-body density matrix was proposed by the author of [22]. The method was applied succesfully to Li and further work is necessary to extend it to atoms with more than one valence electron.

3.3. The homogeneous

electron gas

We would now like to address the simulation of the electron gas, both homogeneous and inhomogeneous. Most quantum chemistry methods which treat electronic correlation cannot be applied to extended systems. For these systems the study of electronic correlation is done with density functional theory (DFT) [l]. In DFT the exchange and correlation are treated through a universal functional of the density. Although the theory is formally exact, in practice approximations for the exchange-correlation are needed. The most popular is the local density approximation (LDA) [41] which is to consider that locally the exchange-correlation energy is that of the homogeneous electron gas at the same density. It is then clear that an accurate determination of the correlation energies for the electron gas is needed. There has been a large amount of work, using different methods, to determine such energies and the most accurate are considered to be those of Ceperley and Alder using GFMC [42,43]. LDA is expected to be more accurate for homogeneous systems; for inhomogeneous systems one must obtain more accurate approximations such as generalized gradient approximations (CGA) [4447]. Although the starting point for such approximations is the homogeneous electron gas, it is expected that the knowledge of correlation energies and pair correlation functions for the inhomogeneous case can provide an insight into the problem and make further improvements. For this purpose jellium surfaces have been widely studied as a test case for

83

P.H. AcioWJournal of Molecular Structure (Theochem) 394 (1997) 75-85 Table 5 Comparison

of the jellium surface energies LDA;

rx

LDA:

- 730 110

2.07 2.66

LM’

- 602 _

a From [50]. b From [St]. ’ The DFT calculation of [5 11 using the Langreth-Mehl d From [52]. ’ Fixed-node DMC calculations from [48,49]. ’ This work.

- 484 _

non-local

In this work we present the results of a fixed-node diffusion Monte Carlo (FN-DMC) simulation for a jellium slab at the average density of magnesium (r,y = 2.66). In this model the positive background density is given by zso

n+(z)=

(26) { 0 z>o This positive background density will create an external potential v(z) in which the electrons move. The Hamiltonian of the system (in atomic units) is then: f$f-

C ‘+ rii

Cv(zJ+const.

i
We

used

Table 6 Comparison

periodic

(27)

i

boundary

conditions

and

QMC - 464’ 320’

orthorhombic supercell, square in the xy-plane, with dimensions L, = Ly = 5.57r, A, Lz = 7.63r,? A and s = LJ4, with periodic boundary conditions in all three directions. The trial wavefunction is of the SlaterJastrow type (Eq. (21)), where the two-body is derived from the RPA approximation for the homogeneous electron gas, modified to account for the anisotropy of the system as proposed by Li et al. [49] in their calculation for r,? = 2.07. In Table 5 and Table 6 we present the surface energies and work functions for r,r II 2.07 and rs = 2.66 and compare them with results obtained using different methods. One can see from Table 5 that for rs = 2.07 the DMC result is in agreement with the DFT [51] using the Langreth-Mehl [44] non-local functional and lower than the Fermihypernetted-chain (FHNC) calculation of Krotscheck et al. [52]. For rs = 2.66 the DMC result approaches that of the FHNC which is expected to be higher than the result using the Langreth-Mehl functional. These results suggest that at lower densities, where electronic correlation should be a more important factor, the current approximations for DFT are less accurate and should be improved.

3.4. Jellium sut$aces

H=-

- 222 383

functional

new approximations. As different approaches have given different results, it is our hope that the investigation of such system using QMC can conclusively estabilish the correlation energies of this system for a range of electron densities [48].

no

FHNCd

an

of the jellium work functions

rr

LDA;

LDA;

LM‘

FHNCd

QMC

2.07 2.66

3.87 3.66

3.79 _

4.12 _

3.39 _

_

‘From [53]. b From [51]. ’ The DFT calculation of [51) using the Langreth-Mehl d From [52]. ’ This work, from fixed-node DMC calculations.

non-local

functional.

2.6Se

84

P.H. Acioli/Joumal of Molecular Structure (Theochem) 394 (1997) 75-85

4. Conclusion We have reviewed briefly VMC and DMC methods and how the electronic correlation is treated in these methods. We showed that DMC in the fixed-node approximation yields energies with chemical accuracy and has been successfully applied to atoms, molecules, clusters, solids and surfaces, and should be considered an important tool to study the electronic structure of systems where correlation effects are important. We have also presented results from a fixed-node DMC simulation of the inhomogeneous electron gas (jellium surface) at T,~= 2.66. The surface energy from our calculations is higher than the results using the Langreth-Mehl non-local functional, but still lower than the FHNC values. The work function is lower than the LDA value obtained by Lang and Kohn [53].

Acknowledgements I would like to acknowledge the collaboration with D. Ceperley without whom this work would not have been possible. This work has been supported by NSF grant DMR-91-17822 and CNPq (Brazil). The computational work used the workstations and supercomputers of the National Center for Supercomputing Applications and the Cray-C90 of the Pittsburgh Supercomputer Center.

References [l] P. Hohenberg and W. Kohn, Phys. Rev., I36 (1964) B864. [2] N. Metropolis, A.W. Rosembluth, M.N. Rosembluth, A.H. Teller and E. Teller, J. Chem. Phys., 21 (1953) 1087. [3] C.J. Umrigar, K.G. Wilson and J.W. Wilkins, Phys. Rev. Lett., 60 (1988) 1719. [4] P.J. Reynolds, D.M. Ceperley, B.J. Alder and W.A. Lester, Jr., J. Chem. Phys., 77 (1982) 5593. [5] M.H. Kalos, Phys. Rev., 128 (1962) 1791. [6] D. Ceperley, in J.G. Zabolitsky (Ed.), Recent Progress in Many-Body Theories, Springer-Verlag, Berlin, 1981, p. 262 [7] D.M. Ceperley and B.J. Alder, J. Chem. Phys., 81 (1984) 5833. [8] D.M. Ceperley and B.J. Alder, Phys. Rev. Lett, 45 (1980) 566. [9] J.C. Slater, Phys. Rev., 34 (1929) 1293. [lo] C.C.J. Roothan, Phys. Rev., 97 (1951) 1474. [I l] I. Shavitt, in H.F. Schaeffer III (Ed.), Methods of Electronic Structure Theory, Plenum Press, New York, 1977.

[I21 P.E.M. Siegbahn, in G.H.F. Diercksen and S. Wilson (Eds.), Methods in Computational Molecular Physics, D. Riedel. Dordrecht, 1983. [ 131 S.F. Boys and NC. Handy, Proc. R. Sot. London, Ser. A, 3 10 (1969) 43. [I41 NC. Handy, J. Chem. Phys., 58 (1973) 279. [15] K. Schmidt and J.W. Moskowita, J. Chem. Phys., 93 (1990) 4172. [ 161 Y.K. Kwon, D.M. Ceperley and R.M. Martin, Phys. Rev. B, 48 (1993) 12037. [17] J.B. Anderson, Int. Rev. Phys. Chem., 14 (1995) 85. [18] B.L. Hammond, W.A. Lester, Jr., and P.J. Reynolds, Monte Carlo Methods in Ab Initio Quantum Chemistry, World Scientific, Singapore, 1994. [19] C.J. Umrigar, M.P. Nightingale and K.J. Runge, J. Chem. Phys., 94 (1993) 3657. [20] J. Olsen and D. Sundholm, reported in A.-M. MArtersonPendril, Phys. Rev. A, 43 (1991) 3355. [21] J.E. Holmstrom and L. Johansson, Ark. Fys., 40 (1969) 133. [22] P.H. Acioli and D.M. Ceperley, J. Chem. Phys., 100 (1994) 8169. [23] A. Veillard and E. Clementi, J. Chem. Phys., 49 (1968) 2415. [24] M. Caffarel and D.M. Ceperley, J. Chem. Phys., 97 (1992) 8415. [25] B. Chen and J.B. Anderson, in press. [26] L. Hammond, P.J. Reynolds and W.A. Lester, Jr., J. Chem. Phys., 87 (1987) 1130. [27] D.M. Ceperley, J. Stat. Phys., 43 (1986) 815. [28] S. Fahy, X.W. Wang and S.G. Louie, Phys. Rev. Lett., 61 (1988) 1631. [29] L. MitaS, E.L. Shirley and D.M. Ceperley, J. Chem. Phys., 95 (1991) 3467. [30] S. Fahy, X.W. Wang and S.G. Louie, Phys. Rev. B, 42 (1990) 3503. [31] J. Grossman and L. MitB, Phys. Rev. Lett, 74 (1995) 1323. [32] L. MidS and R.M. Martin, Phys. Rev. Lett., 72 (1994) 2438. [33] L. MitaS, Phys. Rev. A, 49 (1994) 441 I [34] J. Grossman and L. MitaS, Phys. Rev. B, 52 (1995) 16735. [35] G.B. Bachelet, D.M. Ceperley and M.G.B. Chiocchetti, Phys. Rev. Lett., 62 (1989) 2088. [36] X.-P. Li, D.M. Ceperley and R.M. Martin, Phys. Rev. B, 44 (1991) 10929. [37] P.H. Acioli and D.M. Ceperley, unpublished work. [38] A. Bosin, V. Fiorentini, A. Lastri and G.B. Bachelet, unpublished work. [39] R.J. Needs and M.J. Godfrey, Phys. Rev. B, 42 (1990) 10933. [40] JANAF Thermodynamic Tables, M.W. Chase, Jr., C.A. Davies, J.R. Downey, Jr., D.J. Frurip. R.A. MacDonald and A.N. Syverud (Eds.), J. Phys. Chem. Ref. Data, I4 (Suppl. I) (1985) 535. [41] W. Kohn and L.J. Sham, Phys. Rev., 140 (1965) Al 133. [42] D. Ceperley, Phys. Rev. B, 18 (1978) 3 126. [43] D.M. Ceperley and B.J. Alder, Phys. Rev. Len, 45 (1980) 566. [44] D.C. Langreth and M.J. Mehl, Phys. Rev. B, 28 (1983) 1809. CD. Hu and D.C. Langreth, Phys. Ser., 32 (1985) 391. [45] J.P. Perdew and Y. Wang, Phys. Rev. B, 33 (1986) 8800; 40 (1989) 3399.

P.H. AcioliIJournal

of Molecular Structure (Theochem) 394 (1997) 75-85

[46] J.P. Perdew. J.A. Chevary, S.H. Vosko, K.A. Jackson, M.R. Pederson, D.J. Singh and C. Fiolhais, Phys. Rev. B, 46 (1992) 6671; 48 (1993) 4978. (471 A.D. Becke, Phys. Rev. A, 38 (1988) 3098. [48] P.H. Acioli and D.M. Ceperley, Phys. Rev. B, 54 (1996) 17199.. [49] X.-P. Li, R.J. Needs, R.M. Martin and D.M. Ceperley, Phys. Rev. B, 45 (1992) 6124.

85

[50] N.D. Lang and W. Kohn, Phys. Rev. B, 1 (1970) 4555. [51] Z.Y. Zhang, D.C. Langreth and J.P. Perdew, Phys. Rev. B, 41 (1990) 5674. [52] E. Krotscheck, W. Kohn and G.-X. Qian, Phys. Rev. B, 32 (1985) 5693. [53] N.D. Lang and W. Kohn. Phys. Rev. B, 3 (1971) 1215.