Mathematical Challenges from Theoretical/Computational Chemistry

Many-body problems in physics are often treated by a Monte Carlo (MC) approach (e.g., Hammersley and Handscomb, 1964; Kalos and Whitlock, 1986). The Monte Carlo method is statistical and draws its name from the famous gambling casinos of Monaco because of the role of random numbers or coin tosses in the method.

Problems handled by Monte Carlo are generally of two types, probabilistic or deterministic, depending on whether they are connected with random processes. In the probabilistic case, the simple Monte Carlo approach is to observe the occurrence of random numbers, chosen in a way that they directly simulate the physical random processes of the original problem, and to infer the desired solution from the behavior of these random numbers. In the deterministic case, the power of the Monte Carlo approach is t he capability of carrying out numerical calculations in cases where the equations that describe the essence of a problem and its underlying structure are not solvable by alternative means. The underlying structure or formal expression also describes some unrelated random process, and therefore the deterministic problem can be solved numerically by a Monte Carlo simulation of the corresponding probabilistic problem.

The essential feature common to all Monte Carlo computations is that at some point one will need to substitute for a random variable a corresponding set of values with the statistical properties of the random variable. The values that are substituted a re called random numbers. They are not really random, however, because if they were it would be impossible to repeat a particular run of a computer program. An absolute requirement in debugging a computer code is the ability to repeat a particular run o f the program. If real random numbers were used, no calculation could be repeated exactly, and attempts to check for errors would be extremely difficult. It is essential that one be able to repeat a calculation when program changes are made or when the program is moved to a new computer.

For electronic computation it is desirable to calculate easily by a completely specified rule a sequence of numbers as required that will satisfy reasonable statistical tests for randomness for the Monte Carlo problem of interest. Such a sequence is ca lled pseudorandom and clearly cannot pass every possible statistical test.

Most of the pseudorandom number generators now in use are special cases of the relation (Heermann, 1986; Kalos and Whitlock, 1986)

x a x+ a x+ ... + aj xj + b (mod P)

One initiates the generator by starting with a vector of j + 1 numbers x, x, . . . , xj. The generators are charact erized by a period that in the best case cannot exceed P . The len gth of and the statistical properties of the pseudorandom sequences depend on the values of a, b, and P.

With the choice of a = 0, j 1, and b = 0, one obtains the multiplicative congruential generator

x x (mod P)

Recent work has shown that a set of parameters and P can be chosen with confidence to give desired statistical properties in many-dimensional spaces (Kalos and Whitlock, 1986). There are many other ge nerators available; the interested reader should consult the references.

With parallel computers and supercomputers capable of very large calculations, very long pseudorandom sequences are necessary. In addition, there remains the desire to have reproducible runs. The question of independence of separate sequences to be us ed in parallel remains a research issue.

Quantum Monte Carlo (QMC), as used here, refers to a set of methods to solve the Schrödinger equation (exactly, in two of the three variants listed below) to within a statistical error by random walks in the many-dimensional space. These methods-- variational MC, diffusion MC, and Green's-function MC--are based on the formal similarity between the Schrödinger equation in imaginary time and a multidimensional classical diffusion equation.

Variational Monte Carlo (VMC)

For T(R) a known approximate (trial) wavefunction, where R is the 3N set of coordinates of the N-particle system, VMC (Kalos and Whitlock, 1986; Hammond et al., 1994) uses the Metropolis algorithm (Hammersley and Handscomb, 1964; Heermann, 1986; Kalos and Whitlock, 1986; Hammond et al., 1994) to sample |T|. Ther efore, any expectation value with respect to this trial function can be computed including the variational energy of T. To begin a calculation, an initial distribution of walkers is generated. In order to create a distribution of |T|, these walkers take a series of Metropolis steps to equilibrate, followed by another series of Metropolis steps at which the local energy

is calculated for each walker; here H is the Hamiltonian. Averaging the local energy over the walkers at the sampling points yields the variational energy, which is an upper bound to the exact energy of the ground state. An alternative strateg y is to minimize the variance of the local energies. A strength of the VMC method is the capability of treating trial functions that depend explicitly on interelectronic coordinates--there is no integral problem associated with the use of trial functions containing such coordinate dependence for many-electron systems.

Diffusion Monte Carlo (DMC)

If one multiplies the time-dependent Schrödinger equation in imaginary time by T and rewrites it in terms of a new probability distribution f(R,t) (R,t)T(R), one obtains

where D / 2m, EL is the local energy, and ET (the trial energy) represents a constant shift in the zero of energy. At large simulation (imaginary) time, the function (R ,t) tends to the ground state wavefunction.

The algorithm (Hammond et al., 1994) is initiated with a distribution (ensemble) of several hundred walkers taken from f(R,0) = |T(R)|, which is then evolved forward in time after a sequence of equilibration steps. The three terms on the right-hand side then correspond to diffusion with diffusion constant D, a drift term assoc iated with the trial function, and a branching term that derives this designation from the DMC equation being, in the absence of the first two terms on the right-hand side,a first-order kinetic equation. Because f can, in general, assume both pos itive and negative values, which would preclude the interpretation of f as a probability for fermion systems, one alternative is to impose the nodes of the ground state wavefunction T on so that f is always positive. This is the fixed-node approximation, which can be shown to give an upper bound to the ground state energy.

After sufficiently long simulation time in which the steady-state solution of the DMC equation has been attained, the ground state energy is obtained as the average value of the local energy now averaged over the mixed distribution (R)T(R). This is an improvement over VMC and associated sampling from |T| because DMC contains information on the exact ground state (in the fixed-node approximation). The ground state energy has the zero-variance property of QMC; that is, in the limit that T approaches , the variance of the MC estimate of the energy approaches zero.

Green's Function Monte Carlo (GFMC)

The integral equation form of the Schrödinger equation can also be modeled by a stochastic process, which leads immediately to the consideration of Green's functions. GFMC approaches (Kalos and Whitlock, 1986) have been investigated for the time-i ndependent as well as the time-dependent Schrödinger equations. In practice, however, there are only convergent when the Fermi energy is close to the Bose energy and the trial function has sufficiently accurate nodes. New directions with GFMC that overcome these limitations still encounter limitations that restrict their application, at present, to first-row atoms. Nevertheless, GFMC remains an area for continued scrutiny.

Research Opportunities

The various forms of Monte Carlo for solving the Schrödinger equation (VMC, DMC, and GFMC) could each be improved by better sampling methods. Wavefunctions typically used for importance sampling often recover at best 80 to 95% of the correlation en ergy, the energy difference between the Hartree-Fock approximation and the nonrelativistic limit, for molecules consisting of first-row atoms. Beyond the first row, computer time dependence on atomic number, which has been estimated as z^ where = 5.5 to 6.5, is a major limiting factor that has led to the introduction of analytical functions to describe inner-shell electrons (i.e., pseudopotentials or effective c ore potentials). Another area in which improvements are eagerly sought is in methods that go beyond the fixed-node approximation. Currently, these methods are limited to atoms and molecules containing no more than 6 to 10 electrons.

References

Hammersley, J.M., and D.C. Handscomb, 1964, Monte Carlo Methods, Methuen, London.

Hammond, B.L, W.A. Lester, Jr., and P.J. Reynolds, 1994, Monte Carlo Methods in Ab Initio Quantum Chemistry, World Scientific, Singapore.

Heermann, D.C., 1986, Computer Simulation Methods, Springer-Verlag, Berlin.

Kalos, M.H., and P.A. Whitlock, 1986, Monte Carlo Methods, Volume 1: Basics, John Wiley & Sons, New York.

NAS Home Page | NAP Home Page | Reading Room | Report Home Page