THEO
CHEM
ELSEVIER zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
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 zyxwvutsr
I\k,J, we get the state
l+(t))=
5
i=l
e Xp (-t((Yi-(Yg ))I~,X~;I\kO )
(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 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
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’) =
dRG(R’, t’; R, t)+(R, t)
J
time step approximation
for G.
1
G(R’, t’ R, t) =
(4aDAt)“N/2
(3
G(R’,t’;R,t) can be thought of as a transition probabilMonte Carlo methods. An expression
for
ity to R’ given the initial configuration R. The diffuG(R’,t’;R,t) valid for small time steps is given by zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIH
sion 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 prox exp
- $V(R’)
+ V(R)
(6)
>I
cess is repeated untilf(
[
d ,t) reaches its asymptotic state
using
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’)=
dRG(R’,t’R,t)f(R,t)
s
where
(7)
f $5~)
= \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
*(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
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 = zyxwvutsrqponmlkjihgfedcbaZYXWVU
(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
t)F(R, t) zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
(8)
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
78
P.H. AcioWJournal
of Molecular Structure (Theochem) 394 (1997) 75-85
@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
(15)
energy obtained will be an upper bound on the ground
state energy. A simple proof of the variational charIn the limit of very accurate wavefunctions both estiacter
of this approximation is given in [4]. Thus, in the
mators yield the same answer.
limit
that
the trial wavefunction nodes are the same as
In this method there is a systematic error due to the
those
of
the
ground state wavefunction, the procedure
approximation of Green’s function for a short time
will
project
the
correct ground state wavefunction and
step, but Reynolds et al. [4], showed that if A? is
energy.
It
is
clear
now that the role of the trial waveadjusted so that the acceptance ratio is larger than
function
is
very
important.
First, because a good trial
99%, this error will be small. This systematic error
wavefunction
will
reduce
the
fluctuations in energy,
can be eliminated by sampling the exact Green’s funcmeaning
faster
convergence
and smaller variance.
tion of the system
Second, because if the nodes are good then the energy
I zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
obtained will be very close to the true ground state
G(R, R, At) =
(16)
1 + At(E - EL)
energy.
This method is called Green’s function Monte Carlo
The practical implementation of this approximation
(GFMC) and it was developed by Kalos et al. [5], but
is very simple: we just include the extra step to kill a
will not be discussed here, as it will not be used in this
configuration whenever it reaches or crosses a node
work.
(i.e. *o(R)*(R) = 0). The typical error in this approxDMC produces the exact ground state for a system
imation is about 10% of the difference between the
composed of bosons, but one must be careful in treatVMC and DMC energies 161.
ing fermions, because of the antisymmetrization
of the
Beyond the fixed-node approximation
there are
wavefunction,
in which case the mixed distribution
procedures known as release-node [7] or transientf(R,t) is not always definite positive. In the next secestimate [8] methods, which allow one to obtain a
tion we will discuss the problems with fermions, and
more accurate estimate of the ground state properties.
the so-called fixed-node approximation.
The release-node method allows the walkers to cross
the nodes, not killing them instantly, and letting them
2.3. Fermi statistics and thefied-node
approximation
continue for some time. But their contributions
to
averages will change sign each time the walkers
For the diffusion Monte Carlo algorithm to work,
cross a node. In this method the antisymmetry is takthe mixed distribution f(R,t) must be always positive
ing into account exactly, but the process is numeridefinite (aO(R)\k(R) 2 0). For bosons this is true, as
cally unstable and in practice is restricted to systems
the ground state wavefunction is nodeless, and a good
in which the nodes are very accurate or to systems
trial wavefunction should be constructed in this way.
where the Pauli principle is relatively unimportant.
For fermions, because of the antisymmetric character
of the wavefunction, both aO(R) and ‘P(R) will have
2.4. The trial wavefunction
nodes, and usually they will not be the same. That
means that there will be regions in space where
Above we briefly discussed the role of a good trial
f(R,t) will be negative, making direct sampling out
wavefunction
in QMC. As in any other variational
of the question.
methods, an accurate wavefunction
means lower
One way to avoid this problem is to restrict the
energy. In QMC a good trial wavefunction
has the
random walks to regions in space where the trial
additional role of reducing the variance of the energy,
and the true ground state have the same sign. This is
meaning shorter simulations
to achieve a desired
equivalent to imposing the boundary condition that
statistical error and the very important
role of
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:
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:
dln*(r,, =0)
dru
=
I
394 (1997) 75-85
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 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
(IS)
i for like spins
(21)
U(r) =
sr
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
and a and zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIH
b are adjustable parameters and det( I ).
det( ] ) are Slater determinants
for spin-up and
spin-down electrons, respectively. The parameter zyxwvutsrq
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
81
of Molecular Structure (Theochem) 394 (1997) 75-RS
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
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 non4f (R t> zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
-=
- DV2f (R, t) + DV.F(R, t)
local pseudopotential. In order to have a well-behaved
Hamiltonian the functions al and b, must satisfy
at
- 6 f(R, t) + 4R) (22)
>
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. AcioW Journal of M olecular Structure (Theochem) 394 (1997) 75- 85
Table 5 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
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
sut$aces
{
no
zso
0
z>o
(26)
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:
H=-
f$f-
Cv(zJ+const.
C ‘+
rii
i<j
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.
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
n+(z)=
- 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].
3.4. Jellium
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
_
2.6Se
‘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.
_
84
P.H. Acioli/Joumal of M olecular 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.