In the Classroom
Quantum Chemistry without Wave Functions: Diffusion Monte Carlo Applied to H and H2+
W
Heather L. Cuthbert and Stuart M. Rothstein* Department of Chemistry, Brock University, St. Catharines, ON L2S 3A1, Canada; *
[email protected] In the traditional quantum chemistry course, some effort is devoted to the exactly solvable problems, such as the H atom, arriving at an analytic solution to the Schrödinger equation (expressed in atomic units): {1/2 ∇ 2ψ + V ψ = E0 ψ
(1)
Here V is the potential energy, and ψ and E0 are the exact wave function and energy, respectively. The treatment of molecules is less satisfactory. For the simplest case, H 2+, normally the analytic solution is not discussed in any detail. Instead, the system is considered within the LCAO–MO approximation, as a prelude to the study of more complex molecules. In atomic units, ψLCAO–MO = N [exp({ α ra) + exp({ α r b)]
(2)
where N is a normalization constant, and ra and r b are the distance from the electron to the protons, labeled a and b, respectively. Even allowing for the orbital exponent α to be treated as a variational parameter, the LCAO–MO approximation for H2+ gives rise to a substantial error for the energy at the equilibrium internuclear separation (R nn = 2a0) (see Table 1). In this paper we show how to solve the H atom and H2+ Schrödinger equation by using computer simulation methods. The technique is called “diffusion Monte Carlo” (1). The theory can be presented in an intuitive manner and, for simple systems, the Monte Carlo simulation itself is sufficiently straightforward that it makes an attractive exercise. Prior knowledge of the form of the wave function is not required. The instructor can assign these systems to students in a basic quantum chemistry course as a computer exercise; the results will be in agreement with the exact energies. Programming skills needed to perform the simulation are minimal. At this university there is a computer laboratory associated with a third-year-level quantum chemistry course. The students are taught basic Fortran 77 for half the term, including intrinsic functions, do-loops, conditional statements, vectors, and two-dimensional arrays. After midterm, as a project, they program the diffusion Monte Carlo algorithm for H and H2+. Extension to H2 is easily achieved, but not recommended for large classes because in this case the algorithm is relatively CPU intensive. Theory In an early review of quantum Monte Carlo, Ceperley and Alder (2) credit Fermi in 1945 with the idea of using stochastic methods to solve the Schrödinger equation. In his early simulations on simple systems, Anderson (1, 3) rekindled interest in using random walks to solve electronic structure problems in chemistry. 1378
Let’s consider the time-dependent Schrödinger equation in atomic units and (for the sake of simplicity) in one dimension: ∂ψ / ∂ t = 1⁄ 2 ∂ 2 ψ / ∂ x 2 – (V (x) – E 0 ) ψ
(3)
This equation can be viewed as a diffusion equation to which a first-order rate term is added: ∂ C/ ∂ t = D ∂ 2C/ ∂ x 2 – kC
(4)
Here D is the diffusion coefficient, C is the concentration of imaginary “walkers” (particles located at x), and k is the rate constant for their disappearance or multiplication. In the case of the Schrödinger equation, D is equal to 1⁄ 2 and k is the difference between V and the exact energy. The gist of the diffusion Monte Carlo method is that one can “simulate” the Schrödinger equation by a process of “diffusion”, whereby each walker is moved at random a distance determined by D, followed by a process of “branching”, where each is allowed to multiply or disappear with a probability (branching probability) depending on the value of k. The Schrödinger equation is “solved” in the sense that the exact energy E0 is estimated from the branching probabilities (see eq 9), and the walkers are distributed in space with probability equal to ψ, the exact wave function. Algorithm One can simulate the Schrödinger equation for any system in a manner analogous to that described in ref 4. For our purposes we are concerned with only a single electron. Choose a specific small value of the time-step τ. Also place M walkers at arbitrary locations in space. Say, for walkers randomly placed within a cube 1.0 au on a side, centered on the origin (the proton for the case of H, the midpoint of the bond for the case of H2+): x (m)o = (ζ – 0.5)
(5)
and the same for the y and z-coordinates. (The superscript m denotes the mth individual walker, and the subscript o the oth iteration.) In eq 5 and in the equations below, ζ is a uniform random number from the interval (0,1), obtained from so-called pseudo-random-number generators, which are almost universally available. These are designed such that a new value of ζ is generated on each individual call to the function. First we do “diffusion”. Advance each walker’s x-coordinate as follows: x (m)i+1 = x(m)i + τ 1/2 χ
(6)
where χ is a random number drawn from a normal distribution of zero mean and unit variance. Repeat this for the y- and zcoordinates. It is necessary to save the previous coordinate values from the ith iteration, as they are required in eq 8, below.
Journal of Chemical Education • Vol. 76 No. 10 October 1999 • JChemEd.chem.wisc.edu
In the Classroom
One may generate χ from two calls to the standard random-number generator as follows: χ = sqrt( { 2 ln( ζ ))cos( πζ )
(7)
Next we do “branching”. For each walker, compute the potential energy V (Ri) and V (Ri+1), using the values of its old and new coordinates, respectively. The potential energy is readily computed from the Cartesian coordinates: for H atom it is given by (in obvious notation) {1/r and for H 2+ by {1/ra – 1/rb + 1/Rnn. Also compute branching probabilities for each walker: Bm = exp{{ [V (Ri(m))+V (Ri+1(m))]τ/2}
(8)
and find its average value . The average branching probability provides an estimate of the exact energy : E0 = { ln {}/ τ
(9)
Now make Nm copies of the new coordinates {x (m)i +1, y (m)i +1,
z(m)i +1}, where
Nm = int {Bm / + ζ }
(10)
(The symbol “int” means take the integer part of the quantity in parentheses.) We now have M ′ walkers, the sum of all these Nm values. The subsequent “moves” (iterations) will thus consist of a variable (M ′) number of walkers. The energy estimate is always given by { ln {}/τ, where is now the average value of B m over the set of M ′ walkers. However, one must use the following modified version of eq 10 to prevent a random “extinction” or “explosion” of the “population” of walkers, which would otherwise be an inevitable consequence of the well-known laws of stochastic processes: Nm = int {(Bm/) ? (M/M ′) + ζ }
(11)
Repeat these moves a large number of times to reduce the statistical error of the overall estimate of E 0 obtained by simple averaging of the individual estimates (eq 9). (One must exclude the initial iterations that are required for “equilibration”, that is, the walkers reaching their stationary distribution, ψ.) To determine the standard error of this overall mean, divide all iterations into several (L) large blocks of the same size, compute the block averages, and then combine them using the ordinary statistical formula for independent observations (which the block averages practically are) to get the error bar of the E0 estimate: σ /L1/2. The grand-mean estimates of E0 done at a single small value of the time-step are normally of sufficient accuracy for
Table 1. Ground-State Energies for 1- Electron Systems System
LCAO–MO
H H 2+ (R nn = 2a 0)
b
{ 0.5865
Diffusion Monte Carloa
Exact Energy
{ 0.499 ± 0.007
{ 0.5
{ 0.601 ± 0.005
{ 0.6026...
c
NOTE: All energies are in atomic units. aThe time-step τ = 0.01 au, and initially M = 1000 walkers. The energy is the grand mean of 10 blocks, after one block is ignored to allow the process to equilibrate. Each block is the average energy estimate (eq 9) from 500 moves or “iterations” of the walkers. bRef 7 . cRef 8.
teaching purposes. (Some caution is required here, as too small a value ruins the efficiency of the algorithm: the statistical error of the energy is proportional to 1/τ1/2 [5].) Applications The algorithm was implemented for the hydrogen atom and H2+. (A copy of our codes written in Fortran 77 for Silicon Graphics computers appears on the World Wide Web [6 ] and on JCE Online.W ) The Monte Carlo energy estimates are reported in Table 1. Our algorithm is exact only in the τ → 0 limit, so the energy estimates have an underlying time-step bias. We selected a sufficiently small value of the time-step such that the bias is not visible; the Monte Carlo estimates are all within statistical error of the exact energy. Note W A copy of our codes written in Fortran 77 for Silicon Graphics computers is available on JCE Online at http://jchemed.chem.wisc.edu/ Journal/issues/1999/Oct/abs1378.html.
Literature Cited 1. 2. 3. 4. 5. 6.
Anderson, J. B. J. Chem. Phys. 1975, 63, 1499–1503. Ceperley, D.; Alder, B. Science 1986, 231, 555–560. Anderson, J. B. J. Chem. Phys. 1976, 65, 4121–4127. Vrbik, J.; Rothstein, S. M. J. Comput. Phys. 1986, 63, 130–139. Rothstein, S. M.; Vrbik, J. J. Comput. Phys. 1988, 75, 127–142. Department of Chemistry; Brock University: St. Catherines, ON; http://chemiris.labs.brocku.ca/~chemweb/faculty/rothstein/jce/ index.html (accessed Jul 1999). Also posted as supplementary material on JCE Online.W 7. Weinhold, F.; Chinen, A. B. J. Chem. Phys. 1972, 56, 3798–3801. 8. Wind, H. J. Chem. Phys. 1964, 42, 2371–2373.
JChemEd.chem.wisc.edu • Vol. 76 No. 10 October 1999 • Journal of Chemical Education
1379
