3304
J, Phys. Chem. 1980, 8 4 , 3304-3306
Thermodynamics of Liquid and Supercritical Water to 900 O C by a Monte-Carlo Method Seamus
F. O'Shea"
Department of Chemistry, University of Lethbrldge, Lethbridge, Alberta T7K 3M4, Canada
and Peter R. Tremalne Research Chemistry Branch, Atomic Energy of Canada Limited, Whiteshell Nuclear Research Establishment, Pinawa, Manitoba ROE 7L0, Canada (Received: Ju/y 75, 7980)
The Matsuoka-Clementi-Yoshimine water-water configuration interaction potential was used in a Monte-Carlo simulation of liquid and supercritical water at temperatures up to 900 "C and several densities. Isochoric heat capacities derived from the simulation agree with experimental values to within 5 J mol-l K-I. Discrepancies ~, between the calculated and experimental values for the internal energy, -5 kJ mol-l at p = 1.0 g ~ m -are largely independent of temperature but decrease with decreasing density, consistent with the neglect of long-range polarization effects and three-body interactions in the model.
Introduction Ab-initio theoretical calculations of the properties of liquid water have been carried out by a number of workers using both the Monte-Carlo (MC)l-12and molecular dynamics (MD)13-16approaches. With the exception of one early study,14these calculations have been restricted to 25 "C and densities near 1.00 g/cm3; no simulations of liquid water above 100 "C or of supercritical water have yet been reported. Accurate experimental values for the thermodynamic properties of water have been determined at temperatures up to 1000 "C and pressures as high as 1000 MPa.17J8The high-temperature data, which include a wide range of densities, provide a unique opportunity to study the success of simulation techniques in describing bulk properties under a variety of conditions. Such calculations are also of practical interest in that they can, potentially, be used to provide thermodynamic data for the experimentally inaccessible, extreme conditions which can occur in industrial and geological systems. In principle both the MD and MC methods are exact for classical systems. For molecular fluids, such as water, which contain at least one heavy atom, quantum effects are believedlg to be sufficiently small that they can be either neglected or, at worst, treated by an expansion in powers of h2. In practice, the quantum corrections are usually neglected. The technical difficulties which arise in a finite calculation for these near-classical fluids are reasonably well understood.2&22The major difficulty has been the lack of an accurate formulation for the waterwater pair interaction potential. Recently, Lie and Clementi6 and Swaminathan and Beveridges have obtained superior results at 25 "C and 1 g/cm3 by using the analytical pair potential developed by Matsuoka, Clementi, and Y o ~ h i m i n etheir , ~ ~ so-called configuration interaction (MCY-CI) model. Here we report the results of similar Monte-Carlo calculations on high-temperature liquid water and on supercritical steam. Methods and Results The calculations were done by using the conventional MC A cubic unit cell containing 108 molecules was used with periodic boundary conditions, and the dimension of the cube was chosen to give the required density. For the initial run at each density, performed at 600 "C, the cell was divided into 33 = 27 subcubes of equal volume and four molecules were placed at random in each 0022-3654/80/2084-3304$0 1.OO/O
subcube. In this way a roughly uniform density was generated for the starting configuration. The system was then equilibrated by allowing at least lo5 attempted moves according to the Metropolis algorithm20before beginning to collect the configurational energies. Subsequent calculations at other temperatures and the same density used the final configuration from the 600 "C simulation to initiate the equilibration procedure. Values for the configurational internal energy (eq 1)and E(config) = (l/N)CEj J
(1)
the configurational heat capacity (eq 2) were obtained from C,(config) = ( l / N R P ) C ( E j E(config))2 J
(2)
the molecular energy, E., of each of the configurations produced in a preequilibrated run of N 2 8 X lo5 attempted moves. In extending the Metropolis algorithm to include rotational degrees of freedom we took particular care to ensure that the sampling in angle space is random and uniform. The rotations and translations were tested separately during the equilibration runs in order to obtain rotational and translational amplitudes which yielded an acceptance ratio of 0.5.20v21In the final runs, composite moves were used, and, in scaling the amplitudes to maintain this acceptance ratio, the ratio of rotational to translational amplitudes was held fixed at the value determined during equilibration. The intermolecular interactions were truncated by using the spherical cutoff technique,20and in every case the cutoff distance was taken to be half the cell side. The total thermodynamic internal energy, E , and the isochoric heat capacity, C,, also include nonconfigurational contributions from the translational and internal kinetic energy of the molecules. These were assumed to be identical with those listed by McBride and Gordonz4for the real gas in its standard state, so that E = E(config) + Eo(gas) (3) C, = C,(config) + C,O(gas) Configurational energies and heat capacities from the MC simulation, and the resulting values for E and C,, are listed in Table I. Experimental values of comparison, also listed in the table, are from the semiempirical equation of state reported by Haar et al.18 for the range 0 It I900 "C and 0 S P I1000 MPa. The differences between the exper0 1980 American Chemical Society
The Journal of Physical Chemistry, Vol. 84, No. 24, 1980 3305
Thermodynamics of Liquid and Supercritical Water
TABLE I: 13xperimentala -and Calculated Parameters for Water internal energy, kJ mol-' -. tP C! p l g cm+ P/MPa E( config) Eb 2 51 3001 600 900
1.00
200
0.865 0.86 5 0.374
600
600
-36.6 -26.6 -1 9.4 -14.0 -29.0 -19.3 -12.0
6.4 505 1051 1483 1.9 616 100
1.00 1.00 1.00
-29.1, -11.9, 4.1, 17.7, -17.0, 4.3, 11.5,
-34.1 -1 7.2
heat capacity, J mol-' K-l
C"b
CJconfig)
70.8, 55.6, 49.9, 49.0 62.3, 51.1, 43.9,
45.5 27.9 18.6 17.3 35.6 19.8 12.6
-1.0 (16.9)
-20.7 0.4 10.8
74.2 54.7 (53.9) 59.8 51.1 45.1
a Experimental values, in italics, are from the equation of state of Haar et a1.l' 'Values in parentheses are extrapolations, subject to ~ n c e r t a i n t y . ~ ~Calculated parameters are the sum of the configurational term and the real gas kinetic termsz4
0
J-lot,
60 04
08
06
p/g
10
, ,
200
I
I
400
1
I
600
i
0
cm-3
t/"C
Flgure 1. Deviations from the p = 1.00 g c w 3 experimental loschor. 6~ = €(config) .t Eo(gas)- E(exptl), SC, = C,(config) Cvo(gas)C,(exptl). Error bars in SC, correspond to tho fluctuations in C,(conflg) during the last 2 X lo5 moves in the simulation.
+
imental values and those from our simulation, 6E and 6C,, are plotted in Figures 1and 2, which show values for the p = 1.00 g cm-3 isochor and the 600 " C isotherm, respectively. We estimate the precision of the calculated internal energies to be -0.5 kJ mol-l, based on the following method. Each rim was segmented into units of lo5 moves; local averages over these units were then used to form the global average, and the standard deviations provide the quoted error limit. It is important to distinguish these uncertainties from the systematic errors which may arise because of the use of small repeat units, truncation of the potential, and the possible existence of undetected longlived states in the calculations, all of which are discussed below. Reliablle estimates of the uncertainty in the calculated heat capacities are more difficult to provide because C, is known1° to converge slowly. Consequently, we simply based our estimated precision on the maximum variation observed in CUduring the final 2 X lo5 configurations of each run. The results are shown as error bars in Figures 1 and 2. In common with all of the models used to date, the MCY-CI potential energy function contains a number of deficiencies which result in calculated pressures in very poor agreement23 with experimental values. Pressure calculations are so time consuming that they would have sharply reduced the number of temperatures and densities we could study and, as a result, were inot included in the simulations. Also, since both MC simulations14J6and experimental s t u d i e P show that the radial distribution functions for water retain little detailed structure as the
Flgure 2. Deviationsfrom the 600 "C isotherm. &€and GC,are defined in the Figure 1 caption.
temperature is raised above 100 "C, we chose not to calculate them.
Discussion there has been considerable discussion in the literature concerning technical details of the simulation of polar fluids, particularly water. These center on two topics: the method of handling the long-range interactions, and the number of attempted moves needed to ensure a representative statistical sample for calculating bulk thermodynamic properties. While there is every reason to believe that these are real difficulties in calculations for room-temperature systems, the elevated temperatures studied here have the effect of reducing the long-range correlaitions and increasing the mobility of the sampling point in phase space. For this reason we feel that these considerations are less important here. Our results for a trial calculation at 25 "C gave a slightly higher energy and a lower heat capacity but are entirely consistent with the detailed analysis given by Mezei et al.1° With increasing temperature and decreasing density, the convergence became increasingly rapid, and no evidence of the multistate behavior found at room temperature was observed in any of the high-temperature runs. On the basis of the results reported here, it appears that the discrepancy between the calculated and experimental values of the internal energy is largely independent of temperature but deceases rapidly with decreasing density, This is consistent with the expected deficiencies of the model. Short-range polarization effects are considered since they are included in the Hartree-Fock calculations which are fitted to construct the potential. However, the model potential does not include terms which reduce to
3306
The Journal of Physical Chemistry, Vol. 84, No. 24, 7980
the correct form for the long-range polarization terms, Furthermore, no attempt is made to account for three-body interactions, many of which are due to polarization in this case. The error induced by these omissions should decrease with the decreasing density but be less strongly affected by temperature changes, in agreement with the accuracy of the simulated internal energies in Table I and in Figures 1 and 2. We note that, at the high temperatures used here, some configurations may include 0-0 separations shorter than those used in the fitting of the potential energy function. One danger in such situations is that the model may have “holes” or regions of low energy in physically unreasonable configurations; the form of the MCY potential makes this unlikely in the present case. The effect of systematic errors in the potential at short separations is probably small. The heat capacity is harder to assess. Any residual convergence difficulties will be more evident in the heat capacity than in the internal energy; furthermore, the experimental data are significantly less reliable. The accuracy of the simulated C, data at the lower densities is particularly significant in this regard, since the kinetic contribution is well defined under these conditions as is the experimental data for comparison. Systematic errors in the simulated internal energies at high densities, of the type discussed above, should not affect C, since they are nearly temperature independent (C, = (dE dT),). In Table I, we have used the data for the nonrigid2 rather than the rigid molecule gas to correct the internal energy and the heat capacity for kinetic effects. This makes relatively little difference for the internal energy but is much more important for the heat capacity at high temperatures, where the real gas heat capacity deviates substantially from 3R. The real gas approximation is undoubtedly valid at low densities and, in fact, yielded excellent agreement for supercritical steam at p = 0.374 and 0.865 g cm-3 (Figure 2). Increasing the density perturbs the internal vibrational and, at densities approaching 1.0, this may cause the kinetic contributions to decrease toward the rigidmolecule ideal gas values. More accurate experimental C, data at p = 1.00 and temperatures above 500 “C would be required to resolve the point.29 Our results support the contention that numerical simulation is potentially useful for estimating thermodynamic parameters of systems under extreme conditions. They also suggest that, whatever its faults,30the purely theoretical MCY-CI model for the pair interaction potential of water molecules is remarkably good for this purpose. Finally, we stress that using the data for water at elevated temperatures available from the engineering and geochemical literature broadens the range of material available
d
O’Shea and Tremaine
for testing models and that this important source of information should not be neglected.
Acknowledgment. We are indebted to Dr. L. Haar for supplying data from his equation of state and to the staff of the WNRE computer center for their assistance with the MC simulations. Helpful discussions with Dr. M. L. Klein are also gratefully acknowledged. The work was funded by Atomic Energy of Canada, Limited.
References and Notes (1) J. A. Barker and R. 0. Watts, Chem. Phys. Lett., 3, 144 (1969). (2) J. A. Barker and R. 0. Watts, Mol, Phys., 28, 792 (1973). (3) R. 0. Watts, Mol. Phys., 28, 1069 (1974). (4) H.Popkie, H. Kistenmacher, and E. Ciementi, J . Chem. Phys., 59, 1325 (1973). (5) H. Klstenmacher, G. Lie, H. Popke, and E. Clementi, J. Chem. Phys., 61, 546 (1974). (6) 0. Lle and E. Clementl, J . Chem. Phys., 62, 2195 (1975). (7) G. Lie, E. Clementl, and M. Yoshimine, J. Chem. Phys., 84, 2314 (1976). (8) S. Swaminathan and D. L. Beveridge, J. Am. Chem. Soc., 99,8392 (19771. (9) M.-%zei, S. Swamlnathan,and D. L. Beveridge, J. Am. Chem. Soc., 100, 3255 (1978). (10) M. Mezei, S. Swaminathan, and D. L. Beverldae, J. Chem. Phys., 71, 3366 (1979). (11) C. Pangali, M. Rao, and B. J. Berne, Chem. Phys. Lett., 55, 413 (1978). (12) M. Rao, C. Pangali, and B. J. Berne, Mol. Phys., 37, 1733 (1979). (13) A. Rahman and F. H. Stilllnger, J. Chem. Phys., 55, 3336 (1971). (14) F. H. Stilllnger and A. Rahman, J. Chem. Phys., 57, 1281 (1972). (15) F. H. Stillingr and A. Rahman, J. Chem. fhys., 80, 1545 (1974). (16) R. W. Impey, M. L. Klein, and I. R. McDonald, J. Chem. Phys., in press, 1981. (17) See, for example: (a) H. C. Helgeson and D. H. Kirkham, Am. J. ScL, 274, 1089 (1974); (b) K. Tijdheide In “Water, A Comprehensive Treatise”, Vol. I, F. Franks, Ed., Plenum Press, New York, 1972, Chapter 13. (18) L. Haar, J. Gallagher, and G. S. Kell, Proc. Int. Conf. Prop. Water Steam, 9th, 7979, in press. (19) J. G. Powies and G. Rickayzen, Mol. Phys., 38, 1875 (1979). (20) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J . Chem. Phys., 21, 1087 (1953). (21) J. P. Valleau and S. G. Whittington, “Statistlcal Mechanics Part A: Equilibrium Techniques”, B.J. Berne, Ed., Plenum Press, New York, 1977. (22) J. A. Barker and D. L. Henderson, Rev. Mod. Phys., 48, 587 (1976). (23) 0. Matsuoka, E. Clementi, and M. Yoshimine, J. Chem. Phys., 84, 1351 (1976). 1241 B. J. McBrlde and S. Gordon. J . Chem. fhvs.. 35. 2198 (1961). i25j A. H. Narten and H. A. Levy, J . Chem. Phis., 55, 2263 (i97ij. (26) A. J. C. Ladd, Mol. Phys., 33, 1039 (1977). (27) D. J. Adam, E. M. Adam, and G. J. Hills, Mo.l phys., 38, 387 (1979). (26) . . G. S. Kell In “Water, A Comprehensive Treatise”, Vol. I, F. Franks, Ed., Plenum Press, New York, 1972, Chapter 10. (29) The p = 1.OO lsochor for C, from Haar’s thermodynamic surface” shows a minimum near 500 OC rather than contlnuing to decrease monotonically toward the real gas values with Increasing temperature as mlght be expected from ref 17b and the calculated results In Table I. No experimental C, data have been reported at these high temperatures and densities, and the extrapolated values for C, are partlcularly sensitive to small errors in the surfa~e.”~‘~ (30) I. R. McDonald and M. L. Klein, J. Chem. Phys., 68, 4875 (1978).
-
