4866
J . Phys. Chem. 1987, 91, 4866-4873
The Diffusion Monte Carlo Method for Quantum Systems at Nonzero Temperatures
D. F. Coked and R. 0. Watts* Research School of Physical Sciences, The Australian National University, Canberra, Australia (Received: August 28, 1986)
An analogy between the Bloch equation for the equilibrium density matrix and the diffusion equation with source and sink terms is exploited to develop a nonzero temperature quantum simulation algorithm. The method is derived in terms of simple one-dimensional problems and then applied to gaseous neon and the water dimer.
Introduction In recent years there have been rapid developments in quantum simulation algorithms.' The methods can be categorized in two classes, zero temperature algorithms, in which a single quantum state is generated, and methods for sampling a thermal distribution of quantum states. Typical examples of the first class are the Green's function Monte Carlo method first proposed by Kalos2 and the random walk Monte Carlo method suggested by Metropolis and Ulam3 and implemented by Grimm and Storer4 and A n d e r ~ o n . ~Applications include the simulation of the crystalline and liquid phases of helium: the electron gas,7 metallic hydrogen; as well as other system^.^ More recently, quantum simulation has been extended to molecular quantum mechanics10and shows some promise of being competitive with conventional CI methods.'' Most methods proposed for sampling thermal distributions of quantum states are based on finite temperature path integral techniques.I2 Such methods have been used to study the properties of liquid He4,I3neon gas,14electrons in fused salts,15the properties of argon clusters,I6 the solvation of H atoms and muonium in water,17 and the electronic and vibrational spectra of molecules.16 The methods used to study these problems are based on an integral equation for the many-body density matrixI2 and are either variants on the classical Monte Carlo a l g ~ r i t h mor ' ~ use Fourier expansions about linear paths.'6g20 The random walk Monte Carlo method for generating a single quantum state depends upon a correspondence between the time-dependent Schrodinger equation in imaginary time and the diffusion equation with source and sink term^.^-^ Such an analogy also exists between the Bloch equation for the equilibrium density matrix12s2' and the diffusion equation. In this case, the time variable is replaced by the inverse temperature /3 = l / k T . We explore the equivalence in this paper, modifying the quantum random walk algorithm to generate samples drawn from the quantum equilibrium many-body density function. The following section outlines the relevant statistical mechanics and describes the computational algorithm used. Next, the accuracy of the method is established by using the simple harmonic oscillator and Morse oscillator as examples. The method is then used to generate quantum distributions for neon gas and water dimer. Comparisons with other methods, such as the matrix squaring technique of Storer,22the direct evaluation of the Slater sum,23and the computation of Feynman paths'2-20 indicate the general applicability of the random walk algorithm. Random Walk Evaluation of the Density Matrix In the coordinate representation, the density matrix is defined as
Present address: Department of Chemistry, Columbia University, New York, N Y .
0022-3654/87/2091-4866$01.50/0
where $,(r) and E , are the the eigenfunctions and eigenvalues of the Hamiltonian operator
and the sum includes both discrete and continuum eigenstates. The diagonal components of this matrix represent the probability of finding a system, which is in thermal equilibrium at temperature l / @ , in a volume element dr centered at point r in configuration space
dr,r;/3) = P(r;B)
(3)
(1) Kalos, M. H. In Recent Progress in Many Body Theory, Zabolitzky, J. G., dellano, M., Fortes, M., Clark, J. W., Eds..;Springer: Berlin, 1984; Lecture Notes in Physics No. 142. Kalos, M. H. Monte Carlo Methods in Quantum Problems; Reidel: Dordrecht, 1984. (2) Kalos, M. H. Phys. Rev. 1962, 128, 1791. J . Compt. Phys. 1967, 2 . Phys. Rev. A 1970, 2, 250. (3) Metrowlis. N.: Ulam. S. J . A m . Stat. Assoc. 1949. 44. 335. (4) Grimm, R.; Storer, R. G . J . Compt. Phys. 1969, 4, 230. J . Compt. Phys. 1971, 7 , 134. (5) Anderson, J. B. J . Chem. Phys. 1975.63, 1499. J . Chem. Phys. 1976, 65,4121. Int. J . Quantum Chem. 1979,15, 109. J . Chem. Phys. 1980, 73, 3897. Anderson, J. B.; Freihaut, B. H. J . Compt. Phys. 1979, 3Z, 425. (6) Whitlock, P. A,; Ceperley, D. M.; Chester, G. V.; Kalos, M. H. Phys. Rev. B 1980,19, 5598. Phys. Rev. B 1980, 21,999. Lee, M. A,; Schmidt, K. E.; Kalos, M. H.; Chester, G. V. Phys. Rev. Letr. 1981, 46, 728. (7) Ceperley, D. M.; Alder, B. J. Phys. Rev. Lett. 1980, 45, 566. (8) Ceperley, D. M.; Alder, B. J. Physica 1981, 1088, 875. (9) Ceperley, D. M.; Kalos, M. H. In Monte Carlo Methods in Statistical Physics, Binder, K., Ed.; Springer: Berlin, 1979. (10) Ceperley, D. M.; Alder, B.bJ. J . Chem. Phys. 1984, 81, 5833. (11) Oksiiz, I. J . Chem. Phys. 1984, 81, 5005. (12) Feynman, R. P. Statistical Mechanics; Benjamin: New York, 1972. Feynman, R. P.; Hibbs, R. H. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. (13) Pollock, E. L.; Ceperley, D. M. Phys. Rev. B 1984, 30, 2555. (14) Thirumalai, D.; Hall, R. W.; Berne, B. J. J . Chem. Phys. 1984, 81, 2523. (15) Parrinelio, M.; Rahman, A . J . Chem. Phys. 1984, 80, 860. (16) Freeman, D. L.; Doll, J. D.J . Chem. Phys. 1984,80, 5079. J . Chem. Phys. 1985, 82, 462. (17) de Raedt, B.; Sprik, M.; Klein, M. L. J . Chem. Phys. 1984,80, 5719. (18) Thirumalai, D.; Berne, B. J. J . Chem. Phys. 1983, 79, 5029. J . Chem. Phys. 1984, 81, 2512. (19) Barker, J. A. J . Chem. Phys. 1979, 70, 2914. (20) Doll, J. D.; Freeman, D. L. J . Chem. Phys. 1984, 80, 2239. (21) Hill, T. L. Statistical Mechanics: Principles and Selected Applications; McGraw-Hill: New York, 1956. (22) Storer, R. G. J . Math. Phys. 1968, 9, 964. Klemm, A . D.; Storer, R. G. Aust. J . Phvs. 1973. 26. 43. (23) Fosdick, L. D. J . Math. Phys. 1962, 3, 1251. Fosdick, L. D.; Jordon, H . F. Phys. Reu. 1966, 143, 58.
0 1987 American Chemical Society
The Journal of Physical Chemistry, Vol. 91, No. 19, 1987 4867
Diffusion Monte Carlo Method for Quantum Systems The density matrix can also be written in the momentum representation in which case P(P,P’;P) =
n
+n*(p’) +,(PI
(4a)
That is, p(p,p’) is the double Fourier transform of the density matrix in the coordinate representation. The eigenfunctions &(p) are the energy eigenfunctions in the momentum representation, and the components of p are the eigenvalues of the momentum operators for the system. In the limit of high temperatures, that 0, the correspondence principle allows the momentum is as p eigenvalues to be identified with the classical momenta of the particles. As with the function p(r,r’), the diagonal components, P(p;p), of p(p,p’) represent the probability of finding the system with momenta p at inverse temperature 0. If the functions P(r) and P(p) are known, they can be used to calculate the structural and thermodynamic properties of the system. Both these functions can be obtained from the partial Fourier transform of the density matrix
-
F(p,r;P) = jei/hpr’p(r,r’;B) dr’ by a second transform, either from r to p’ or from p to r’. The function F(p,r;P) can be written
repeated for all members of the ensemble. Finally, the temperature is decreased, by putting (3 = Po Ab, and the whole process re~eated.~,~ In the case that the diffusion algorithm is being used to simulate solutions to the time-dependent Schrodinger equation, the initial configuration does not influence the final results. It can be shown in that case that the population evolves into one whose particle positions represent a sample distributed as the ground-state wave f u n c t i ~ n Also, . ~ ~ ~it is clear that if the simple algorithm outlined above is used for many increments A& the population of replicas will eventually decay to zero. Again, this problem arises when the algorithm is applied to the time-dependent Schrodinger equation. However, a straightforward population-stabilizing algorithm has been developed to overcome such problem^.^ Both the initial conditions and the population stabilization method have to be reconsidered when the diffusion algorithm is applied to the Bloch equation. At high temperatures, the function F(p,r;@)can be written in terms of the classical Boltzmann distribution
+
F(p,r;P) =
e-isHciassel/hpr
= e-fl~~,2/2m,e-6v(r)er/ hpr
(10) (1 1)
where the quantities pI are the classical momenta of the particles. In the present work, the initial spacial distribution of replicas is determined by using the classical Monte Carlo algorithm of Metropolis et al.,25giving a population that is a sample from the distribution exp(-pv). Every member of the ensemble can then be given a weight chosen from the multidimensional Gaussian distribution
as was shown by K i r k ~ o o d , ~and ’ , ~satisfies ~ the Bloch equation
aF/ap = -A F(p,r;@)
(7)
with the initial condition F(p,r;O) = ei/hpr It is the fact that F satisfies the Bloch equation that enables us to use the quantum random walk method to generate numerical results for the density matrix. Writing eq 7 in full gives
aF
h*
- = C-V?F ap 2mi
- V(r)F
(9)
an equation that is isomorphic with the diffusion equation for a system with source and sink terms. The inverse temperature, @, represents the time variable and F the multidimensional concentration. Note that the momentum p appears parametrically in this equation. The random walk method for solving equations of this type has been described in detail by Grimm and Storer4 and A n d e r ~ o n . ~ Suppose at temperature Po a set of replica systems is established, for example 500 copies of a harmonic oscillator or 300 copies of a typical configuration for a small cluster of water molecules. The algorithm consists of propagating this initial ensemble along a trajectory in “time”, or more correctly for this example, 8. A discrete increment in inverse temperature is chosen A@. Every particle in the first replica is then displaced by an amount, Ari, chosen at random from a Gaussian distribution with variance 4DiAp, where Di = h 2 / 2 m i is the “diffusion coefficient” for the particle. Once the new particle positions have been determined, the total potential energy of the new configuration, V(r), is calculated. Next, follows the “birth/death” step. If V(r) > 0, the replica is destroyed with probability exp(-ApV); if V(r) < 0, it is replicated with probability exp(-APV) - 1. This process is (24) Kirkwood, J. G . Phys. Reo. 1933,44, 31.
and a phase factor, exp(ip.r/h). If the ensemble is now propagated along the /3 trajectory, with this weight being retained by both the initial replica and all its descendents, the required sample from the distribution F(p,r;P) is generated. In practice, this weighting procedure is not needed. Momentum variables only appear parametrically in the Bloch equation and hence the weights do not change as the simulation progresses. Also, any properties of the system are calculated by averaging over all values of the classical momenta. As we shall show shortly, the simple Gaussian form for the initial momentum distribution allows these integrals to be done analytically. However, for the moment we consider every replica to be weighted as in the preceding paragraph. The second problem associated with the algorithm is the fact that the population size decreases steadily over the duration of the simulation. This decrease results in large fluctuations in calculated properties, particularly at low temperature (Le. large p). As quantum effects are particularly important in this region, some method for reducing these fluctuations has to be introduced. In the case of the time-dependent Schrodinger equatio?, the simulation is designed to generate an eigenfunction of H , and details of the transient behavior are unimportant. Hence, as discussed in detail by Ander~on,~ the zero of energy can be adjusted during the simulation by introducing a parameter VR. The birth/death decisions are made as above, but in terms of the quantity V(r) - VR. A steady-state population is obtained when VR is equal in value to the lowest eigen energy of the Hamiltonian. In the case of the Bloch equation, the reason for performing the simulation is to generate the transient behavior of F(p,r;P)-that is, to calculate the properties of the system as a function of temperature. Consequently, the simple V , method outlined above has to be modified. The method used here is to implement the VR scheme, but keep track of the values taken by this pa-
(25) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M . N.; Teller, E.; Teller, M. J . Chem. Phys. 1953,21, 1087.
4868
The Journal of Physical Chemistry, Vol. 91, No. 19, 1987
rameter as the ensemble evolves. It is clear upon considering the birth/death algorithm that this procedure is equivalent to weighting all members of the ensemble, at step j, by the quantity eXp(vR(j)A@), where VR(j) is the value of VR used at s t e p j to stabilize the population. Consequently, the correct transient behavior is obtained by weighting the population at temperature Dl by the quantity J
W(P) = Wj =
C exp(-APVR(k))
k= 1
(12)
+
where P = Po jAP. The use of this device stabilizes the population and ensures that, in general, the replicas are concentrated near those regions of phase space where F(p,r;P) is largest. Hence the fluctuations in the calculated properties are much reduced. It is also possible to use importance sampling techniques4s5to improve convergence. Suppose we have some approximation to the function F(p,r;P), FT, and use this to define a new function f(P,r,@) = F(p,r;P) FT(p,r;P) Then if this definition is used in the Bloch equation, an equation for f is obtained a! =
w
Coker and Watts We now reconsider step 5 above and show that the required weighting and average over initial momenta can be done analytically. Averages over Position and Momentum We have seen that every replica needs to be weighted by a Boltzmann distribution of momenta and a phase factor. As p is incremented, the replica diffuses away from its initial configuration, r, to reach a new configuration, I-‘. However, the momenta associated with the initial configuration, p, do not change during this process. From eq 5 , we see that the diagonal part of the density matrix for an N particle system, P(r’;p) can be obtained from F(p,r’;P) by an inverse Fourier transform
In the simulation, F(p,r’;P) is represented pointwise by a set of replicas each of which carries a weight determined by the Boltzmann momentum distribution, a phase factor, and the population-stabilizing term W(p),eq 12. Let Dk(r-r’;Po-P) represent the kth replica at any step. The arguments indicate that the kth replica is a descendant of one that started a t position r at the initial temperature Po and has migrated to position r’ by temperature P. With this representation of the replica, we can write eq 13 as p(r’;p) = C’Dk(r-r’;@o-fl)
x
k
As in methods for the ground-state wave function discussed by Grimm and Storer4 and A n d e r ~ o n ,the ~ second term in this equation represents a “drift velocity” effect which enhances diffusive motion toward more favorable regions of configuration space. The third term in the equation represents a modified birth/death process. If FTis chosen appropriately then the simulation algorithm is more efficient. There is one problem with this approach, namely the fact that the required distribution function has to be reconstructed from the ratiof(p,r;@)/F,(p,r;P). In regions where f is small, and F is also small, the division can result in any fluctuations being greatly enhanced. Although we have made some calculations using importance sampling, the errors were sufficiently large that the method was not pursued. We can summarize the algorithm used to generate F(p,r;P) as follows: 1. Use a classical Monte Carlo method to establish an initial spacial distribution of replicas at some sufficiently high temperature. 2. Propagate the ensemble along a trajectory in (3 by using the quantum random walk method. This includes Gaussian displacements of particles and birth/death processes for replica^.^,^ 3. Adjust VR after every increment AB to keep the population approximately constant in size. 4. At temperature 6, weight every ensemble member by the accumulated quantity Wl (eq 12). 5. Associate every replica with a weight consisting of a momentum distribution exp(-PCpr2/2m,) dp and a phase factor exp(ip.r/ h ) . 6. Once the required final temperature is reached, repeat the process with a new initial distribution. With this procedure, averages over a large ensemble can be constructed by using a sequence of smaller ensembles. However, these ensembles should not be too small or fluctuations will be large. The value of V, used to stabilize the population is obtained from the average potential energy of the subsystem at any step. Consequently, the parameter VR(j) will vary from cycle to cycle, and these variations will be large if the population in any cycle is small. Clearly, the larger the fluctuations in VR(j). the larger the fluctuations in W, and hence the bigger the errors in the calculated properties. It follows that the size of the ensemble used in any cycle should be as large as is possible given the limitations imposed by the computer equipment available.
where the prime on the summation constrains k to those replicas whose position at temperature P is r’. Here, we have recognized that neither the initial population, generated by using the classical Monte Carlo position, nor the quantum diffusion algorithm depends in any way upon the momentum variables. The integral in eq 14 can be done analytically, so that the final expression for the configuration space position distribution is
Dk(r-r’;Po-P)
W(P) (15)
Note that by summing over k, that is by averaging a large number of initial replicas, the algorithm automatically averages over the initial Boltzmann distribution of positions. Any average of an operator that is a function of position can be obtained from eq 15 by appropriate multiplication (e.g. by V(r’ ) for the average potential energy) and summing over all replicas. In the case where the operator in configuration space is a differential operator it is more convenient to work in momentum space.’’ Again starting from eq 5 , one can give the probability density in momentum space by
Substituting for the Boltzmann momentum distribution, the phase factor, and integrating as above gives an expression analogous to eq 15. In the case that the single particle momentum distribution is required, the momenta p2, ..., pN are integrated out, and the momentum p1 is written in spherical polar coordinates. Integrating over the angle variables gives P”’(P1;P) =
exp(
E)
sinc
- r,’l] Dk(r-r’;Po-P)
(17)
Note that in this case the summation includes all replicas.
-/‘j
The Journal of Physical Chemistry, Vol. 91, No. 19, 1987
Diffusion Monte Carlo Method for Quantum Systems The kinetic energy of the system can be obtained as
4869
I /,
/ / /
and on making the appropriate substitutions this becomes 1.o
where
i
0’
‘ ,’ Hence in practice, the kinetic energy at temperature @ is obtained by averaging the quantity Q over alj the replicas that survive. Before leaving this section, we remark that all these expressions depend upon the initial temperature being sufficiently high that the classical distribution function is an accurate approximation to the quantum distribution. If such an assumption does not hold, then either the initial temperature must be increased or a better initial distribution has to be used. We will discuss the second possibility in a subsequent section. Some Applications In this section we consider four systems, the simple harmonic oscillator, neon gas, the 0-H bond modeled as a Morse oscillator, and the water dimer. The harmonic oscillator density matrix is know analytically, and this system can be used to assess the proposed simulation algorithm. Several other methods have been used to study dilute neon gas, in particular the matrix squaring algorithm of Klemm and Storer:2 and hence this system can also be used to test the accuracy of the random walk method. The 0-H bond is important in a wide range of molecules whose liquid-state properties are of considerable interest. Consequently, we study this system, and the water dimer, with a view to extending the random walk technique to molecular liquids. Harmonic Oscillator. With the change of variables
the well-known harmonic oscillator Hamiltonian becomes
0.0 1.o
0.0
Reduced temperature
Figure 1. Energy of the harmonic oscillator as a function of reduced temperature. The upper solid curve is the analytical result for the total energy and the lower Qne represents the kinetic (= potential) energy. Dashed curves represent the corresponding classical results. Other symbols are as follows: @,A,+) T(init) = 10, T(fina1) = 0.3; (o,V,X) T(init) = 1, T(fina1) = 0.1. The symbols are ordered as (total energy, potential energy, kinetic energy).
-3.4
3.4
Y
Figure 2. Distribution of displacements for the harmonic oscillator for reduced temperature from 0.70 to 0.30. The solid points represent the analytic results, from eq 24. As this form is relatively simple, we work with reduced energy E* = E / h w and reduced temperature T* = k T / h w , so that @* = hw@. The quantum simulation method described in the previous sections was implemented for this problem. Ensembles with an average population of 5000 replicas were propagated through 100 steps of size A@*. The initial distributions were generated by using a classical Monte Carlo program at Po* = 1.O,0.2, and 0.1, and the trajectories were terminated at @*= 3.33 and 10.0. This range of parameters allowed the influences of starting temperature and step size to be investigated. Figure 1 presents the kinetic, potential, and total energy of the oscillator, as a function of P,for several of the calculations. Also shown in the figure are the predictions of the classical equipartition theory and the analytic result of FeynmanIz
E*pot= Eatin =
1 (1
+ e+’)
4 (1
- e+’)
-
(23)
The potential energy was calculated by averaging y 2 over the ~~
~
~
~~
( 2 6 ) Hirschfelder, J. 0.; Curtiss, C. F.; Bird,
357.
R. B. Mol. Phys. 1956, 52,
configuration space probability density, eq 15, and the kinetic energy was obtained from eq 19. All the results given in the figure were obtained by averaging over 10 independent trajectories, each with an average population of 5000 replicas. The figure shows how the kinetic and potential energies depend upon the starting configuration. When the initial distribution is generated at a temperature To* = 1.0, both the kinetic energy and potential energy differ significantly from the exact result. The kinetic energy is too small throughout the trajectory and the potential energy is too large. As the initial temperature is increased, the differences between the simulated and exact results decreases, and at To* = 10 the simulation is in excellent agreement with the analytic result throughout the trajectory. The figure also serves to demonstrate that the results are essentially independent of step size for the conditions studied here. Also, it is interesting to note that the calculated total energy is always in excellent agreement with the analytic result-that is, the errors in E,, and Eklnare equal in magnitude but opposite in sign. Figure 2 shows the diagonal part of the density matrix in configuration space as a function of temperature. The simulation results are compared with the analytic resultl2 P ( y ; @ * )= (T coth @ * / 2 ) - I f 2 exp(-y2 tanh p * / 2 )
(24)
4870
The Journal of Physical Chemistry, Vol. 91, No. 19, 1987
Coker and Watt\
.
0.1 16
0.058
d 0.70 13.7 27.4 q Figure 3. Momentum distribution for the harmonic oscillator as a function of reduced temperature. The solid points represent the corresponding classical distributions. Clearly, the agreement is excellent. Figure 3 compares the calculated momentum distributions with the classical Boltzmann distribution, showing the difference between classical and quantum mechanics results at low temperature. As a consistency check on the simulated momentum distribution we numerically averaged the kinetic energy, obtaining excellent agreement with the values calculated from eq 23. The problems associated with using too low an initial temperature in the simulation, noted above, can be overcome if a more accurate initial distribution function is used. In the case of the harmonic oscillator, the off-diagonal terms in the density matrix are known analytically'* and can be used to show that
F(q,y;P*) = exp(-)/2(q2
(
+ y z ) tanh p*) exp co:yP*)
4.0
I;"':-.
44. 30.
0.0 1.o
3.0
r/a
Figure 4. Pair distribution function for neon gas calculated at N a 3 / V = 0.0093, for temperatures from 50 to 10 K. Solid curves are from quantum simulation, the dashed lines represent a classical calculation, and the solid points are from Klemm and Storer.22
n
.-> a, c
Y
Y
/
%
P
(25)
a,
This expression can be used in place of that given by eq 11 to generate the initial distribution of replicas, providing the appropriate changes are made to eq 15, 17, and 20. If this procedure is implemented, excellent agreement with the analytic result is obtained by using an initial temperature T* = 1 .O. Low Density Neon Gas. Klemm and Storer22 have used a matrix squaring technique to calculate the pair distribution function for low temperature neon gas. Their method begins with the high temperature density matrix, written in terms of spherical harmonics and the radial wave function, and generates lower temperature results by successive application of a matrix squaring routine. Every time the density matrix is squared, the temperature is halved. Although Klemm and Storer's method is very good for two- and three-body problems, it is not readily extended to more complex systems. We have used the simulation algorithm described in this paper to generate properties of low temperature neon gas and compared them with Klemm and Storer's calculations.22 The gas was simulated by placing two neon atoms in a periodic, cubic box whose dimensions were such that p* = pu3 = 0.0093. A Lennard-Jones (12,6) potential with c / k = 35.7 K and u = 2.789 was used to model the neon pair potential.26 The initial distribution in configuration space was generated by using a classical Monte Carlo algorithm at a temperature of 100 K. Random walks were generated for ensembles of 2500 replicas, with 100 steps Ap between T = 100 K and T = 10 K. A total of 500 separate trajectories were computed, every one starting from a different high temperature configuration. Figure 4 compares the simulated radial distribution functions, at several temperatures, with the corresponding classical distribution functions and the results of Klemm and Storer.z2 The agreement between the two sets of quantum calculations is very good, confirming the accuracy of the simulation algorithm. As expected, there are significant differences between the quantum and classical distribution functions, particularly at lower tem-
a
i
I
/
-
I 0 1 I
i
0
I I
I I
- 10.0 0.0
20.
40.
Temperature (Kelvin) Figure 5. Potential energy of the neon dimer as a function of temperature. The circles, joined by a dashed line, represent classical results; the squares are from quantum simulation. peratures. The most important differences between the two distributions occur in the region of the repulsive wall of the potential. Quantum tunneling enables the neon atoms to interpenetrate, creating classically forbidden configurations. This effect also shows up in the average potential energy of the gas, shown in Figure 5 . Here, the quantum calculation is always more positive than the classical result. The 0 - H Vibmtion. If the quantum simulation method is to be applied to molecular fluids, it is important that methods for modeling intramolecular vibrations be established. The Morse potential V ( r ) De( 1 - e-a(r-ro))2 (26) is frequently used to describe chemical bonds and has been shown elsewhere*' to give a good description of the 0-H vibration in the water molecule. For that reason, we study the isolated Morse oscillator with parameters De = 66,049 K, ro = 0.9572 A, and a = 2.14125 ~ k 'and ,~ then ' the water dimer. The bound state eigenfunctions of the Morse oscillator are known analytically28 and hence the diagonal elements of the (27) Reimers, J. R.; Watts, R. 0. Mol. Phys. 1984, 52, 356. Coker, D. F.; Watts, R. 0. J . Phys. Chem., to be published. (28) Morse, P. M. Phys. Rev. 1929, 34, 57. ter Haar, D. Phys. Rec. 1946,
70, 2 2 2 .
Diffusion Monte Carlo Method for Quantum Systems
2
.-->
4000.
0,
Y Y
t
The Journal of Physical Chemistry, Vol. 91, No. 19, 1987 4871
n
.-c>
4000.
-
s
Y
h
P $
w
8
2000.
2000.
/ /
0.
2000.
4000.
0.
Temperature (Kelvin)
Figure 6. Energy of the Morse oscillator as a function of temperature. The results are represented as follows: (-) total energy by direct summation of Slater sum; (-.-.) kinetic energy from Slater sum; ( - - - - ) potential energy form Slater sum; (---ktc-) corresponding classical results. Other symbols represent quantum simulations, ordered as (total energy, kinetic, potential), (0,+, A) T(init) = 10000 K, (0, X, V) T(init) = 5000 K. Results for T(init) = 7500 K, not shown in the figure, lie between the other two quantum simulations. density matrix and other properties can be calculated directly from the Slater sum, providing the temperature is such that kT
