Quantum Monte Carlo

DAVID CEPERLEY AND BERNI ALDER

In the new numerical techniques, called quantum Monte Carlo

An outline of a random walk computational method for methods, the Schrodinger equation, which exactly describes nonrel-
solving the Schrodinger equation for many interacting ativistic particles, is represented by a random walk in the many-
particles is given, together with a survey of results dimensional space in such a way that physical averages are exactly
achieved so far and of applications that remain to be calculated. Monte Carlo or statistical methods are in fact the only
explored. Monte Carlo simulations can be used to calcu- general methods known for exactly solving problems in many
late accurately the bulk properties of the light elements dimensions, provided only that the problem can be formulated in
hydrogen, helium, and lithium as well as the properties of terms of probabilities.
the isolated atoms and of molecules made up from these Such a numerical approach has only become possible since the
elements. It is now possible to make reliable predictions advent of high-speed computers. In fact, it is particularly adaptable
of the behavior of these substances under experimentally to these machines because algorithms are simple and highly repeti-

Downloaded from www.sciencemag.org on May 17, 2012

difficult conditions, such as high pressure, and ofproper- tive, characteristics that can take fulll advantage of the fast arithmetic
ties that are difficult to measure experimentally, such as capabilities of modem computers. Furthermore, these methods can
the momentum distribution in superfluid helium. For be easily adapted to a variety of computer architectures. They yield
chemical systems, the stochastic method has a number of exact results within statistically determined error bars that decrease
advantages over the widely used variational approach to with the length of the computer run. The principal goal in develop-
determine ground-state properties, namely fast conver- ing algorithms is thus to find ways of increasing the efficiency ofthe
gence to the exact result within objectively established calculations. This can be done in a straightforward way through a
error bounds. technique called importance sampling, which uses previous knowl-
edge to provide a good starting approximation.
The application of statistical methods to quantum mechanical
problems is not without difficulties of its own, the most serious
N THE EARLY DAYS OF QUANTUM MECHANICS, P. A. M. DIRAc being the calculation of systems that have a wave function that is not
observed that the physical laws necessary for the mathematical everywhere positive. Nevertheless, considerable progress over the
theory of a large part of physics and the whole of chemistry are past few years has enabled us to carry out realistic simulations of
completely known and that it is only necessary to find practical systems composed of the light elements. We intend to show here
methods for the solution of the equations for complex systems (1). that there is no practical impediment to realizing Dirac's program
One could have expected that the advent of modern, high-speed for many other many-body systems, although these applications will
computers would have by this time made it possible to perform such require considerably more efficient algorithms and faster computer
computations. However, the endeavor implied by Dirac's statement hardware.
of principle remains largely unfulfilled. Many of the existing numeri-
cal methods provide only a qualitative understanding of the proper-
ties of isolated atoms and molecules or of their collective behavior in Diffusion Monte Carlo
the condensed state. These methods are quantitatively inadequate
either because the approximations they embody cannot be further Around 1945, Fermi remarked that stochastic methods could be
refined or because the numerical scheme converges too slowly. used to solve the Schrodinger equation (2). The earliest recorded
A numerical method developed relatively recently to solve the implementation was carried out in 1949 by Donsker and Kac (3) for
Schrodinger equation for many interacting particles has the poten- the hydrogen atom, but the results were unimpressive because the
tial to realize Dirac's goal. This method is a departure from the lack of an importance function led to very low efficiency. For the
conventional approach to many-body problems in mathematical same reason, an unpublished calculation by Rosenbluth (4) a few
physics; namely, it does not reduce a system with very many degrees years later for the ground state of liquid 4He gave unsatisfactory
of freedom by an approximation to equations of much reduced results. Meanwhile, stochastic processes came increasingly to be
dimensionality. used in the study of neutron transport and classical condensed-
Alternate numerical methods that do not invoke such an approxi- matter systems. An important step for quantum Monte Carlo
mation, such as the configurational interaction method in quantum methods was maded by McMillan (5) in 1965 when he used a
chemistry, nevertheless expand the wave function in a complete set variational method to simulate helium. He showed that a one-to--
of one-body functions, so that once again one only has to deal with one correspondence to a classical simulation could be made if one
low-dimensional mathematical objects. Moreover, the functions assumed a pair-product wave function. In the same period, Kalos (6)
used for the expansion are generally restricted, so that the low- developed what is known as the Green's function Monte Carlo,
dimensional integrations that appear in the theory can be performed which in 1974 culminated in an exact algorithm for calculating the
analytically. The price paid for these restrictions is that, even with ground-state properties of the hard-sphere boson fluid (7). We will
a large number of terms (frequently running into the millions),
the calculations do not converge with the accuracy desired for The authors are at Lawrence Livermore National Laboratory, University of California,
chemistry. Livermore 94550.

7 FEBRUARY I986 ARTICLES 555

describe a simplified version of Green's function Adonte Carlo, or distinguishable particles in the ground state, the wave function
known as diffusion Monte Carlo. can be directly interpreted as a probability density, so that Eq. 1 can
The basis of diffusion Monte Carlo is that the Schrodinger be interpreted as a diffusion and branching process in 3N dimen-
equation written in imaginary time, t, will converge to:o the ground sions. A useful analogy is to bacteria randomly diffusing in a puddle
state exponentially fast. That equation for the wave furiction 4)(R, t) with a diffusion constant of h212m, in which the growth conditions
is are uneven and depend on the position in the puddle (that is, on the
potential energy), which determines the rate ofgrowth or decline of
d~(t)
at = -H = Ej
j m
Vj2 _ [V(R) - ET]+ (1) the bacterial population there. The bacteria do not interact with
each other since the Schrodinger equation is linear. Then Eq. 1 gives
the evolution of the distribution in time. Alternative numerical
where mj is the mass of particle j, V(R) is the total pot.ential energy, methods that rely on tabulating the distribution everywhere (for
R refers to the 3N set of particle positions, and N is tihe number of example, in a grid) will consume an exponentially large amount of
particles (H is the Hamiltonian operator). A constant:, ET, the trial computer time and memory as the dimensionality of the space (the
energy, has for convenience been subtracted from the potential number of particles) increases. Direct simulation by raridom walks,
energy. From a formal expansion of the wave function in a complete which samples the distribution selectively, appears to be the only
set of eigenvectors and eigenvalues, it can be readily demonstrated general way of numerically solving the quantum many-body prob-
that all excited states decay exponentially fast with a dLecay constant lem.
given by the excitation energy from the ground state The rate of Such a computer calculation is set up in the following way (Fig.
convergence to the ground state is hence govemed 1by the lowest 1). An initial ensemble of systems is constructed, usually from a
excited state that has a component in the initially chosen wave classical Monte Carlo calculation with some trial wave function as
function. proposed by McMillan (5). An ensemble consists of a number of
If it is assumed that +(t) is non-negative, as is the caase for bosons "snapshots" of the coordinates of all the particles, let us say of all the

Downloaded from www.sciencemag.org on May 17, 2012

electrons and nuclei. In actual calculations, the ensemble is made up
of about 1000 such snapshots. The evolution is accomplished by
Possible new
Generation configurations considering each snapshot in turn, displacing each of the particles by
Old of new and their New a random amount with a mean square displacement given by
configurations configurations multiplicity, m configurations Th2I2m, where T is the time step. Then branching is done; a number
U of copies of the snapshot equal to the integer part of {exp[-'T/
;;O;" ~~000 2(Vold + Vnew)] + u} is made, where u is a uniformly distributed
s
random number in [0,1]. Thus a new ensemble is generated with a
m=0.8
different number of snapshots. As the ensemble is evolving and its
population is varying, the trial energy ET must be adjusted with a
feedback mechanism so that the population remains stable. If the
population becomes too large, ET is made smaller, and if the
. population diminishes, ET is increased. The value of ET necessary to
stabilize the population is then the ground-state energy. The
_. snapshots generated once steady state is reached (that is, when
a4/at = 0) are then samples of the ground-state wave function.
/
a
0m ~Importance Sampling
For most problems the above algorithm is not satisfactory because
the branching process is uncontrolled. Whenever the potential
energy becomes large and negative (as it will, for example, when an
electron approaches a nucleus), the branching process blows up, and
a huge number ofcopies of that snapshot is created. Luckily, there is
~~~~0 a very simple and elegant way of solving this problem: importance
0
sampling.
0
Importance sampling means changing the underlying probability
m=0.3 distribution in a known way so that the calculation will spend more
time in the important regions. For this purpose the trial function,
Fig. 1. Schematic of the Green's function Monte Carlo calculation with
importance sampling, demonstrating the evolution of the "snapshots." X)T, is introduced as an approximation to the ground-state wave
Illustrated is a three-electron system in a box; the squares represent old function (derived, for example, from a Hartree calculation), and
positions, the circles new ones. The old ensemble consists of four snapshots, f(R, t) = 4T(R)4(R, t) is defined. The Schrodinger equation can be
and the new one consists of three. The process leading to a new ensemble written in terms off by some algebraic manipulations (8), resulting
involves (i) adding to each electron's coordinates a drift term equal to in the importance-sampled equation
'r/2Vln I4TI/M (shown by arrows), where is the time step; (ii) adding a
r

random displacement (shown by a wiggly line) representing the diffusive

step, whose mean square displacement in each dimension is Th2/2m; and (iii) aJ(R,dfw,~
t)t)=
V
Nh
~h Vf
{V - fj in X)T II H4T/4)T -ET()
branching, which creates from zero to several snapshots in the new ensemble,
with the number determined by the integer part of m = [exp(-'rEL) + u],
at j m
where u is a uniformly distributed random number in [0,1] and EL = HOT/ which has a structure very similar to the original Eq. 1 but with
T ET is the local energy. In the fixed-node approximation for dealing
-
with fermions, a snapshot is deleted if the sign of 4T has changed during this some important differences. The first term on the right-hand side is
step. With these three new snapshots, the process is repeated. the gradient of something; therefore it conserves probability and
556 SCIENCE, VOL. 23I
does not cause branching. It represents diffusion with a superim- necessitated by the Pauli exclusion principle. Antisymmetry is built
posed drift. In our analogy, the water in the puddle is not still but in into the trial function by multiplying the pair-product function by a
steady motion, carrying the bacteria around the puddle while they Slater determinant. The Slater determinant is made from single
are diffusing. The second term again gives rise to branching, but electron orbitals obtained from the Hartree-Fock or local density
now the rate is determined by the local energy EL(R) = H40T/'T functional method. These orbitals determine the nodes of the trial
- ET. function, where it changes sign. The procedure for fermions has
The process of simulating this equation proceeds as follows (Fig. evolved in two steps: the fixed-node approximation, and the subse-
1). An initial population is generated according to '4T12. This quent exact algorithm that releases those nodes.
ensemble is evolved by examining each snapshot and diffusively In the fixed-node approximation (9, 12) only one additional rule
displacing each coordinate, as before, but in addition displacing each needs to be added to the previously discussed algorithm: if the
position by a drift term equal to ri2VIIn X~I/m. The effect of the drift random walk crosses a node of the trial function, that is, when
is to push the random walk away from unimportant regions (that is, (kT(old)4T(new) < 0, that snapshot is deleted. This will occur relative-
where the trial function is small), since there the drift velocity is ly rarely, since the drift term will push the walks away from the
large. The number of copies is now calculated from exp{-T/2[EL(O1d) places where the wave function vanishes. In our analogy of the
+ EL(new)]} mapped onto an integer with a random number as bacteria, we must add the condition that the bacteria in the puddle
before. The crucial improvement is that if IT has been chosen close die if they reach a boundary. Thus, this method solves the Schro-
to the ground-state wave function, the local energy will be small, so dinger equation in each nodal region separately. It can be shown
that branching is much less. Importance sampling makes it practical- that the energy so obtained is the best upper bound consistent with
ly possible to solve the Schrodinger equation for several hundred these conditions (13). It has been found that this approximation is
particles. Also, for good trial functions the asymptotic distribution often numerically very accurate because the node locations are not
of snapshots will be equal to the square of the wave function, which crucial for determining the energy. Of course, if the correct node
is just what one would expect physically. locations are known, as in one-dimensional problems, the exact

Downloaded from www.sciencemag.org on May 17, 2012

The statistical error of the ground-state energy is approximately result is obtained. One can apply the fixed-node approximation to
given by [2(Ev - EO)/(TnP*)] /2, which has the familiar dependence calculate any excited state for which a variational principle applies.
on the inverse square root of the number of steps, with a propor- The releasing of the nodes (8, 14) to get their correct locations
tionality constant given by the difference between the variational leads to the exact fermion ground state; however, the computer time
energy of the trial function Ev and the correct ground-state energy required may become exponentially large because of a numerical
EO. Here P* is an effective population of the ensemble (the average instability. Snapshots are not deleted when they hop across a node
number not counting duplicates), and n is the total number of time but now carry a plus or minus sign corresponding to the sign of the
steps. Thus the computational efficiency is determined by the trial function when the walk was begun. The estimate of the wave
accuracy of the trial function (Ev - ET) times the computer time function is the difference in the number of positive and negative
necessary to evaluate the trial function, the drift, and the local energy snapshots that arrive at a given point, and the trial energy is correctly
for a single snapshot. There is a trade-off between trial functions that adjusted when this difference at any given point is constant in time.
are accurate and those that are fast to evaluate. Although the procedure is mathematically correct, the signal-to-
noise ratio for a given amount of computer time decreases exponen-
tially as the positive and negative walks become mixed (Fig. 2); thus
the computer budget may run out before satisfactory results are
Green's Function Monte Carlo obtained. Because the method is unstable, it has been called a
Green's function Monte Carlo (9) is a reformulation of the transient estimate (15). In principle, there are ways of canceling
diffusion process such that no systematic errors due to the finite time positive and negative snapshots to prevent the exponential growth
step T arise. The method is so named because the differential in the population (16). In practice, this is difficult to carry out in
equation is converted into an integral equation, the kernel or many dimensions because the probability that a positive and nega-
Green's function of which is sampled exactly. The procedure is a tive snapshot will have the positions of all the particles identical
generalization of a Monte Carlo method suggested by von Neu- within a possible relabeling is too small. In spite of these difficulties,
mann and Ulan (10) to solve systems of linear equations. In Green's satisfactory results have often been obtained (8). A rigorous and
function Monte Carlo there are several additional elements in the stable method to simulate fermion systems by a stochastic process
algorithm. For one, the time step itself must be sampled for each remains a most challenging problem.
move; the walk does not advance by a predetermined time. Also,
intermediate snapshots that are not legitimate members of the
ensemble are generated. They serve only to sample the correct
Green's function and hence contribute only to the propagation of Other Quantum Monte Carlo Techniques
the walk. In some cases this more complicated procedure is also The variational quantum Monte Carlo method, an adaptation of
more efficient, since a larger average time step results (11). Further- the classical Metropolis algorithm (17), was previously mentioned in
more, the procedure does not require user adjustment of the time connection with finding good importance finctions and beginning
step: the exact result will be automatically obtained. the ensemble. It can also be used to determine the ground-state wave
function in the same sense that the traditional variational methods
(Hartree-Fock or configuration interaction) are used. In the config-
Fermi Statistics uration interaction procedure, the wave function is expanded in a
complete set of functions, each of which is an antisymmetric product
It would appear from the discussion so far that the method is of single-particle orbitals. For variational quantum Monte Carlo
limited to the calculation of systems in which the wave function is methods, there is no such restriction on the basis set, since the
non-negative. A few years ago the method was extended to ground- required integrals are obtained with Monte Carlo rather than
state fermion systems, where the wave function is real but equally analytically. An intriguing possibility is that of a self-learning
positive and negative because of the requirement of antisymmetry mechanism, in which the output of the Green's function Monte
7 FEBRUARY I986 ARTICLES 557
Carlo simulation is used to improve the form of the trial funiction. there have been very few simulations of fermion systems in continu-
For finite-temperature quantum mechanical calculation s, the ous space at finite temperatures, although the zero-temperature
path-integral Monte Carlo is employed. The origin of the metthod is techniques discussed earlier are applicable. Most applications, partic-
based on the observation that the density matrix (at high ternpera- ularly those for the lattice gauge theories of high-energy physics,
tures, the Maxwell-Boltzmann distribution; at low temperatur,es, the have been restricted to lattice and spin systems (20).
square of the wave function) can be factored into a prodLuct of An entirely different method that has also been tried primarily on
density matrices, each at a higher temperature. lattice problems (21) and on one-dimensional systems (22) should
be explored further for three-dimensional Fermi systems in continu-
AMR
(ROle-S-H IRM) _= (RoIeIm/IRI
\
R-IM'IRm)
)(R1 ... IRM-I)(RM-IleC ous space, since Fermi statistics pose no special difficulty with this
method. In this method, pair interactions in the Hamiltonian are
(3) replaced by interactions with an external, random, time-dependent
where P is l/kT (T is the temperature). The chosen numiber of field. By this so-called Hubbard-Stratonovitch (23) transformation,
products, M, is sufficiently large that at the effective temp :rature the many-body calculation is exactly transformed to one'for a system
MT an accurate expression for the density matrix exists, usual] ly such of one-body noninteracting fermions. The one-body problem must,
that the Boltzmann distribution becomes valid. That transformed however, be solved at each step, and this requires considerable
density matrix can then be evaluated by an analogy to a cllassical computer time. A further disadvantage of this scheme appears to be
system of N closed polymer rings ofM links (18). The simtalation that there is no way to allow for the introduction of pair-product
problem at finite temperature is then reduced to finding an elfficient importance functions.
procedure for sampling all energetically contributing polymc-r con-
figurations representing all the intermediate positions Ri of the
links, or, in other words, all contributing paths. For thoat, the
classical Metropolis (17) method is used. This path-integral iMonte Condensed-Matter Applications

Downloaded from www.sciencemag.org on May 17, 2012

Carlo scheme differs from the Green's function Monte Carlo nrnethod The calculation of the properties of liquid 'He at zero tempera-
in that the paths must close on themselves because thermodyn amical ture was the first large-scale application of Green's function Monte
properties are obtained from the trace or the diagonal part of the Carlo (7, 24). It assumed a pair-interaction potential between
density matrix. For the polymer system, one samples a sp)ace of helium atoms deduced from theoretical considerations and experi-
3N x M instead of the 3N x P dimensions (P is the populaltion of mental data. Such an interaction potential can now be accurately
the ensemble) in the Green's function method. Also, there is no calculated directly by Monte Carlo methods (25). Equilibrium
explicit importance sampling in path-integral Monte Carlo me thods. properties, such as energy versus density, pressure versus density,
Boson statistics are introduced by allowing neighboring polynners to crystallization pressure at zero temperature, and the structure factor
cross-link. The efficient sampling of this polymerlike system has obtained by x-ray or neutron diffraction, come out very close to
made it possible to perform accurate simulations (19) of liqui id4He experimental values (24); the differences can be ascribed to the
both above and below the superfluid transition point. In cointrast, inadequacy of the assumed pair potential.
The most spectacular properties of liquid helium, the dynamical
ones resulting from its superfluidity, are difficult to simulate with
Fig. 2. Node release this method. However, many of these unusual transport properties
for the first excited are believed to result from the fact that, in superfluid helium, a finite
state of a particle in a *** fraction of the atoms have condensed into a zero-momentum state.
box. (A) In the initial The difficult neutron-scattering measurements of the momentum
guess, the node was
not properly located at distribution needed (26) to confirm this theory are given in Fig. 3 at
the middle of the box, B 00o000o 1 K. As mentioned above, similar calculations at finite temperltures
but the wave fimction are now being completed.
was made orthogonal a~~~~~~~
0 *n The electrons in the conduction band of a simple metal are often
0
to the ground state. on
The circles represent modeled by replacing the ions by a uniform positive background.
the distribution of the U2 - c: Although the model dates back more than 50 years, there have been
positive population, no convincing calculations of its properties, in spite of many
the squares the nega- attempts, except in the asymptotic limit of high and low densities.
tive population, and C Even the order of magnitude of the melting density is not generally
the solid curve repre-
sents the wave func- 0
0 m
a agreed upon. Wigner (27) predicted that, contrary to the usual
tion (the difference be- 0
*
* 0o
0
situation, a crystal phase occurs in the low-density regime. The
tween the two distri- 0 a properties of the electron gas have been calculated (8, 28) by the
butions). (B) Both the 0 0

variational, fixed-node, and release-node Monte Carlo methods, and

positive and negative 0

populations evolve to- the melting transition has been located. In contrast to the simula-
ward the symmetric, 00 tions of liquid helium, there are no direct experimental results to
nodeless ground state, Om0
compare with these electron gas calculations. Hence, in Fig. 3, the
while the difference OM 0 momentum distribution of an electron gas at a density approximate-
has an improved node 02
or ly equal to that in the conduction band of potassium is compared
location. (C) The 04
process continues, viv- 00 ° with that of an ideal, noninteracting fermi gas. Recent calculations
idly demonstrating QI
Q, similar to those on the electron gas have been performed on 3He, a
on
that the difference is cm CP fermi liquid; the results are in good agreement with experiment
increasingly hard to CP (29). The electron simulation results are now often taken as standard
determine numerically >
as the positive and
input to the approximate solid-state calculations that use local
negative populations density functional theory.
grow. The electron gas has a fairly rich phase diagram. At zero tempera-
558 SCIENCE, VOL. 23I
ture and normal metallic densities, the gas is a regular, spin-paired, to achieve a given statistical error per atom increases only as the
diamagnetic fermi liquid. As the electron density is reduced by a number of particles squared. For very large systems, one can use
factor of 106, the electrons spontaneously spin-polarize. The polar- sparse matrix techniques to lower this power further.
ization increases as the density is lowered by another factor of 106, As a first example of a few-body system, the energy of three
when the electrons undergo Wigner crystallization. Such ferromag- hydrogen atoms for several positions of the three protons has been
netism nearly occurs in liquid He at low temperature and partly calculated (33). The potential surface for this molecule needs to be
accounts for the unusual magnetic behavior of that liquid. These known to establish the barrier for the simplest chemical reaction,
Monte Carlo calculations, involving several hundred particles, have namely the exchange of a hydrogen atom with one in a hydrogen
not yet attained the precision to explore such subtle effects as the
superfluidity of 3He or the superconductivity of metallic hydrogen.
The effects energetically are too small and the relevant length scales
for the phenomena too large compared to the size ofthe system that B
2 ---
can be simulated. However, it is possible to calculate the response of 0

the electron gas to various types of extemal fields, test charges, or 6

impurities. Such calculations are underway. c
1-
Hydrogen, as the simplest of the elements, provides a natural
extension of the previous work on the electron gas in which the
uniform background is replaced by the actual protons. It is much 0
IIl
10 2fl 3XD
more difficult to use a molecular pair potential in Monte Carlo k
calculations of the properties of hydrogen than it is with helium.
The pair interaction for hydrogen is more complex and less certain, 0OX Fig. 3. (A) The momentum
and at high pressures there are nonpairwise additive effects. Thus the 0.05 - distribution, n(k), as a func-

Downloaded from www.sciencemag.org on May 17, 2012

simulations were performed directly with protons and electrons, tion of momentum, k, for liq-
interacting only through their coulomb potential. Both the protons uid 4He at 1 K. The momen-
tum distribution as obtained
and the electrons have a sizeable quantum motion. This is taken into
3 ~~~~~~by
0 Monte Carlo calculations
account in the simulation by letting the protons drift and diffuse, 0 1 2 at zero temperature consists of
but at a rate 1836 times slower than the electrons, which is the ratio a zero momentum compo-
of their masses. nent, the condensate, represented by a delta function, 8(k) (heavy line at the
origin), comprising about 10 percent of the particles, and a normal
These calculations carn be compared with experiments (30) at component, n*(k) (solid curve). Ignoring the temperature-dependence, the
pressures above 106 atmospheres (1 Mbar) reached with a diamond normal component can be compared with the results of neutron scattering
anvil apparatus. However, the molecular-atomic transition to a experiments at 2.27 K (above the superfluid transition; open circles). The
metal, first predicted by Wigner in 1935 (31), has not been observed solid circles represent the measured momentum distribution at 1 K (below
the superfluid transition). The condensate fraction, found by integrating the
experimentally. Simulation of hydrogen (32) in both the molecular difference between the experimental distribution at 1 K and at 2.27 K, is
and atomic phase established that this phase transition occurs at found to be 14 percent, which is comparable to the theoretical value at 0 K.
about 2.8 Mbar. It has also been established that the protons (B) The momentum distribution of an interacting electron gas calculated by
undergo a melting transition at low temperature but under astro- the Monte Carlo method at zero temperature and at a density approximately
physical conditions (108 Mbar). The excellent agreement between equal to that of the valence electrons in potassium under standard conditions
(dots and heavy line) compared with that of an ideal fermi gas (dashed line).
the theoretical equation of state in the molecular phase and the
experimental one is shown in Fig. 4. The energy resolution achieved
with these calculations of 0.001 Rydberg per atom is also enough to
Pressure (Mbar)
determine crudely the pressure (roughly 1 Mbar) at which hydrogen
molecules stop rotating, that is, when they become aligned in the
crystal.
These simulations of hydrogen have only scratched the surface of
the interesting properties that could be reliably calculated. It would
be relatively straightforward to obtain band gaps, bond frequencies, 0-
and dielectric properties of hydrogen even though the amount of 0
computer time required is large (10 hours of CRAY-1 time per
computation). Calculations of the properties of mixtures of hydro-
gen and helium, as they occur, for example, in the core of Jupiter, -

and simulation of metallic lithium are further possible extensions of 1.05

CD
CK
LU
this work.

Few-Body Problems 0 2 4 6 8 10
The applications described so far have all been to bulk systems, for Volume (cm3 molr1)
which periodic boundary conditions were used to minimize finite Fig. 4. The energy of hydrogen at zero temperature plotted against volume
system effects. In calculations of the properties of few-particle and corresponding pressure. The results of the diamond-anvil experiment are
systems, such as molecules or clusters, the correct boundary condi- given by the solid curve and an extrapolation of the data is indicated by the
tions are those of an isolated system. In such studies the advantage dashed curve. The Monte Carlo results for the molecular phase are given by
of the Green's function Monte Carlo method lies not only in its open squares and for the atomic phase by closed squares joined by a dotted
line. The double-tangent construction indicated by the straight line deter-
rigor but also in its ability to deal with many more electrons than mines the molecular-to-atomic transition region to be between the two
alternative ab initio methods. The amount of computer time needed horizontal markers at 2.8 Mbar.
7 FEBRUARY I986 ARTICLES S59
molecule. The energy for the barrier as calculated by Monte Carlo energy is insensitive to the value of the wave function where it is
was 9.60 ± 0.05 kcal mol'-, which is 0.3 kcal mol' lower than the needed, at the coalescence point. However, the wave function
best variational configuration interaction upper bound (34) but is can be calculated by Monte Carlo in a very simple way (36). Many
exactly the same as the variational result ifa correction that estimates random walks are started from the desired configuration of the
the rate of convergence of the expansion is applied. particles. When the walks have reached their steady-state distribu-
Another example is a cluster oflithium atoms. Calculations on the tion, the ratio ofthe exact wave function to the trial wave function is
lithium dimer gave more accurate total energies (8, 33) than the proportional to the average population of the ensemble. The
configuration interaction calculation. Similarly, the ground-state method is rigorous, and even at this highly singular point impor-
energy of the trimer in various geometries has been determined, tance sampling permits calculation of the wave function with a
encouraging the study ofmuch larger clusters. If accurate results can relative error less than 0.5 percent.
be obtained for sufficiently large clusters with reasonable computer
time, it will be possible to simulate metal surfaces realistically and
subsequently to study the properties of molecules absorbed on
them. In such cluster calculations the lithium nuclei should not be Conclusion
held fixed, but the geometry of the cluster should come out of the The Green's function Monte Carlo technique is still very much in
calculation. A good approximation would be to assume that these the developmental phase, and applications have been limited to a
nuclei behave classically, which calls for an algorithm that can deal few long-standing problems. Nevertheless, it is already apparent that
with both classical and quantum mechanical degrees of freedom at the simplicity, rigor, and adaptability of this approach to different
the same time. Such an algorithm would be very useful in treating computer architectures make it likely that it will become a standard
any liquid system, such as water, starting with classically behaving computational tool in physics and chemistry. Future algorithmic
oxygen and hydrogen nuclei (the latter requires quantum correc- improvements include finding better numerical methods for dealing
tions but not a full quantum treatment) and electrons that behave with fermions, finding ways of improving the trial wave function to

Downloaded from www.sciencemag.org on May 17, 2012

quantum mechanically. Such calculations would be considerably make importance sampling more effective, computing excited states,
more realistic than the present classical simulations ofwater that are combining quantum and classical Monte Carlo calculations, finding
based on rigid molecules interacting with pair forces. However, an efficient procedure for calculating energy differences, and dealing
with current computers and methods, such a simulation of water with fermions at finite temperature. The systems that could be
would be too costly. One needed improvement is a way of directly simulated range from nuclear matter to plasmas in outer space,
calculating the change in the electronic energy as classical degrees of whenever more than two particles interact.
freedom change. Such a need arises in many applications. Work is
underway on a differential Monte Carlo scheme that computes such REFERENCES AND NOTES
energy differences. I. P. A. M. Dirac, Proc. R. Soc. London Ser. A 123, 714 (I929).
The reason that only relative energies are required is that fre- 2. N. Metropolis and S. Ulam,J. Am. Stat. Assoc. 44, 247 (1949).
quently one is interested only in the difference in the energies of 3. M. D. Donsker and M. Kac,J. Res. NatI. Bur. Stan. 44, SSi (i950).
4. M. Rosenbluth, private communication.
various arrangements of the atom in a molecule to find the 5. W. L. McMillan, Phys. Rev. A 138, 442 (I965).
arrangement of highest stability. Another application would be in 6. M. H. Kalos, ibid. 128, 179I (I962); J. Comp. Phys. 2, 257 (I967).
7. , D. Levesque, L. Verlet, Phys. Rev. A 9, 2178 (1974).
the simulation of atoms larger than neon. The two main problems 8. D. M. Ceperley and B. J. Alder, Pbys. Rev. Lett. 4S, 566 (I980); P. J. Reynolds, D.
for heavier elements are (i) that the time step needed to follow the M. Ceperley, B. J. Alder, W. A. Lester, J. Chem. Phys. 77, 5593 (I982).
9. D. M. Ceperley and M. H. Kalos, in Monte Carlo Methods in Statstical Physics, K.
inner electrons becomes smaller as the core electrons are bound Binder, Ed. (Springer, Berlin, I979), pp. i45-I97.
more tightly, leading to slow convergence, and (ii) that the differ- io. G. E. Forsythe and R. A. Leibler, MTAC 4, I27 (i950).
ii. D. M. Ceperley,J. ComP. Pbys. Sr, 404 (I983).
ence between the fermion ground state and the distribution to I2. J. B. Anderson, J. Chem. Pbys. 63, 1499 (I975); ibid. 65, 4121 (I976).
which the system relaxes upon nodal release becomes larger for I3. D. M. Ceperley, in Recent Progress in Many-Body Theories, J. G. Zabolitzky, Ed.
atoms with more electrons, so that there is increasing difficulty in (Springer, Berlin, I98I), pp. 262-269.
14. M. A. Lee, K. E. Schmidt, M. H. Kalos, G. V. Chester, Phys. Rev. Lett. 46, 728
reliably extracting the difference between the positive and negative (1981).
populations. If the inner-shell electrons could be accurately repre- i5. K. E. Schmidt and M. H. Kalos, in Monte Carlo Methods in Statitil Physuics, K.
Binder, Ed. (Springer, Berlin, 1984), vol. 2.
sented by a pseudopotential, these difficulties would be ameliorated, I6. D. Amow, M. H. Kalos, M. A. Lee, K. E. Schmidt,J. Chem. Phys. 77, 5562 (1982).
since then only the part of the random walk concerned with the 17. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. M. Teller, E. Teller, ibid.
21, I087 (I953).
electrons in the valence region would contribute. The accurate M8. R. P. Fynmann, Statistical Mechanis (Benjamin, Reading, MA, i97).
replacement ofthe inner electrons by a pseudopotential requires that I9. E. L. Pollock and D. M. Ceperly, Phys. Rev. B. 30, 2555 (1984).
20. See M. H. Kalos, Ed., Monte Carlo Methods in Quantum Problems (Reidel,
the core electrons be insensitive to the valence electrons. If for that Dordrecht, 1984).
purpose a nonlocal pseudopotential must be introduced, the Monte 21. J. E. Hirsch, Phys. Rev. B 28, 4059 (I983).
22. G. Sugiyama and S. E. Koonin, in preparation.
Carlo calculation would be much more involved. 23. R. L. Stratonovitch, Sov. Phys. Dokl. 2, 4I6 (I958); J. Hubbard, Phys. Rev. Lett. 3, 77
As a final example of few-body problems, consider the sticking (I9S9) -
24. P. A. Whiitock, D. M. Ceperley, G. V. Chester, M. H. Kalos, Phys. Rev. B I9, 5598
probability of a muon to an alpha particle. In this situation the (1979).
Monte Carlo method can be used not only to obtain energies but 25. R. E. Lowther and R. L. Coldwell, Phys. Rev. A 22, I4 (I980).
26. V. F. Sears, E. C. Svenssen, P. Martel, A. D. B. Woods, Phys. Rev. Lett. 49, 279
also to calculate the value of the wave function in an extremely (1982).
improbable arrangement of the particles, namely when two nuclei 27.
28.
E. P. Wigner, Phys. Rev. 46, I002 (I934).
D. M. Ceperley, Phys. Rev. B 18, 3126 (1978).
are in the process offusing. The total number offusions catalyzed by 29. J. Carlson, R. M. Panoff, K. E. Schmidt, M. H. Kalos, in preparation.
a single muon placed in a deuterium-tritium liquid is experimentally 30. J. van Straaten, R. J. Wijngaarden, I. F. Silvera, Phys. Rev. Lett. 48, 97 (I982).
found (35) to be more than 100. This value is limited by the 31. E. Wigner and H. B. Huntington, J. Che. Phys. 3, 764 (193S).
32. D. M. Ceperley and B. J. Adler, Physica Io8B, 875 (I98I); unpublished information.
probability that an alpha particle will capture a muon immediately 33. , Chem. Pys. 8I, 5833 (I984); F. Mentch and J. B. Anderon, ibid. 80, 2675
after a fusion occurs. Calculation of this capture probability requires (I84).
34.. B. Liu, ibid., 8o, 2675 (1984).
knowing the value of the wave function of a molecule composed of a 35. S. Jones et al., Phys. Rev. Lett. SI, 1757 (i985).
muon bound simultaneously to a deuteron and a triton. Such a 36. D. M. Ccperley and B. J. Alder, Phys. Rev. A 31, 1999 (I985).
37. Supported by Lawrence Livermore National Laboratory under contract to the
three-body calculation is difficult by traditional methods, since the Department of Energy (W-7405-ENG-48).
560 SCIENCE, VOL. 231