J. Phys. Chem. 1994,98, 6231-6233
6231
Theoretical Calculation of the Static Dielectric Constant of Water at High Temperatures and Pressures Saul Goldman' and Chris J o s h Department of Chemistry and Biochemistry and the Guelph- Waterloo Centre for Graduate Work in Chemistry, University of Guelph, Guelph, Ontario, NIG 2WI Canada Evgeny A. Wasserman Department of Geology, University of Bristol, Bristol. BS8 1 RJ, United Kingdom Received: January 21, 1994; In Final Form: April 1 1 , 1994'
We report a computationally efficient method for the quantitative prediction of the static dielectric constant (e) of water at high temperatures and pressures. Values of e in this regime are needed in geochemistry for the prediction of ionic solubilities in hydrothermal fluids. The method involves the application of a recently derived perturbation theory to the SPCE model of water. The experimental database for c in this regime is sparse, and the alternative computational procedure-a direct simulation- is 1-4 orders of magnitude slower than the method we present. We find that the calculation provides quantitatively accurate values of e for water provided y I 2.4, where y is the "dipolar strength function" (=(47r/9)pZp/kT).
Introduction
Theory
The so-called primitive electrolyte solution model' wherein discrete ions (charged hard or soft spheres) are dissolved in a uniform dielectric continuum has proved useful for the prediction of ionic solubilities in hydrothermal fluid^.^^^ The reliable prediction of these ionic solubilitiesis important in geochemistry.3 This application of the primitive model implicitly involves the use of a Born charging process. This, in turn, requires the static dielectricconstant (e) of water and its two logarithmicderivatives with respect to temperature (T) at constant pressure (P).3 Accurate, experimentally determined e's for water exist for limited ranges of P and T in the high-P, high-T regime.4-6 Where applicable, the experimentally determined e's should certainly be used. For state conditions not covered by experiment, simulations or theory can instead be used to predict e. Thus, one of us (E.A.W.) has recently carried out a series of molecular dynamics simulations to get e at high T and P,7 using the SPCE model for water.8 This model was chosen because it is computationally efficient-it is a rigid, nonpolarizable, three point-charge model of water, with the charges embedded in a Lennard-Jones core-and because it was recently shown to provide e's that are remarkably close to those of real liquid water along the liquid-vapor coexistencecurve. Specifically, e(SPCE) = 81.0 and 6 at 300 K and the critical point vs 78.0 and 5.3, respectively,for real water.g Unfortunately, simulations on water models are too time-consuming to serve as a routine method for generating values of E over a large P-T surface. Therefore, a computationally more efficient way to do this was sought. In this paper we do several things. We first show that when a recently derived perturbation theoryI0Jl is applied to the SPCE model, it quantitatively reproducesthe molecular dynamics (MD) values of E for this model, in the high-P, high-T regime. Also, we give an indication of the range of validity of the theory, and we show that along the 400 OC isotherm the theory, together with the SPCE model, quantitatively reproduces the e values of real water in the range 250 5 P I 3000 bar. Finally, we show that the SPCE model is much more accurate than its antecedent, the SPC model,l2 in this application.
We give only a brief outline of the theory; its derivation and the explicit form of all the underlying expressions have been given elsewhere.I0J1 The theory consists of two ingredients: a y expansion of c and a perturbation-theory-based estimate of the angular correlation function that arises in the third-order term of the y expansion. Thus1OJIJ3J4 c = 1 +3y+3yz+3y3
0022-365419412098-6231$04.50/0
l)+.*.
(1)
where y (=(4~/9)p2p/kT)is the dipolar strength function and G(z)is thesecond-order perturbation theoryestimate of theangular correlation function G. The latter is defined by G = PJd[lZ
(co~(Y1z)g(_r,zQlQz) )n,n,
(2)
In the above, p is the permanent dipole moment of the molecules of the fluid, p is the number density, Tis temperature, and k is Boltzmann's constant. Also, 7 1 2 is the angle between the two dipoles of molecules 1 and 2, Q,E (c$&,xJrepresents the Euler angles of particle i, ( )Q,Q, means the unweighed angle average over particles 1 and 2, and g(rlzQ1Qz) is the full distance and angle-dependent pair distribution function. G(2)is obtained from eq 2 by expanding g(rl2Q,Q,)as g(r,zQ1Q,) = go(r1z) + gl(_rl*QlQz)+ gz(LlzQ1Qz) + " *
(3)
and substituting this expansion, truncated at the gz term, into eq 2. In eq 3, go(r12)is the pair correlation function of the isotropic reference state, whose potential is obtained from (4)
where u is the full pair potential, gl is linear in the anisotropic potential u, E ( u - UO),g2 is quadratic in u,, etc. One finds that for eq 3 truncated at g2 only the g2 term contributes to G(2). Specifically,lo.l~.l3-~~ G(2)= GA+ G,
* Address for correspondence at the Department of Chemistry and
Biochemistry. e Abstract published in Advance ACS Abstracts, May 15, 1994.
;;(--
GA e
1 / 2 8 z ~ J d ~ go(rlz)(ua(!lzQIQ,)2 lz
0 1994 American Chemical Society
(5)
cos(~lz))n,n, ( 6 )
Goldman et al.
TABLE 1: Reduced Spherical Multipole Moments for the SPCE Water Model
GA and GB represent the two- and three-body contributions to G(2),and go(r12r13r~~) is the triplet distribution function of the reference system. The problem thus reduces to evaluating GA and GB. For a purely dipolar anisotropy only GBis n o n ~ e r o , ’but ~ for the SPCE model under consideration both G A and GB contribute to e.
The angular integrations in eq 6 and 7 are done analytically by expanding the multipolar u:s in generalized spherical harmonics and applying their orthogonality properties, and go(r12r13r23) is obtained by using a recent refinement16 of the Kirkwood superposition approximation. The radial integrals are done numerically. Equation 1 is a highly accurate route for the calculation of e provided y 6 2.5. For this range of y’s the y3 contribution to e is small, so that a n uncertainty in G(2) (caused by truncation of eq 3 at the g2 term) leads to only a small error in t.
Calculations The calculations were carried out as described elsewhere.I0-l1 The underlying reference-state pair correlation function go(r) was obtained by Monte Carlo Simulations on the Lennard-Jones reference fluid, as described in ref 11. The spherical multipole moments (defined in ref 10 and 11) that we used for the SPCE model are given in Table 1. They can be obtained by the method given in Appendix A of ref 10 or, more simply, by multiplying the previously determined spherical multipole moments of the SPC model by a factor of 0.4238/0.41. This is the ratio of the partial charge assigned to each proton (or the oxygen atom) in the SPCE model relative to the corresponding value in the SPC model. This trivial scaling of the multipole moments arises because the two models differ from one another, only with respect to the magnitudes of their partial charges. The geometries and Lennard-Jones core parameters are identical for these two models.8
Results and Discussion Comparisons of the theoretically determined e’s with values obtained from MD simulations in the highd, high-T regime at both low and high densities are given in Table 2. It is clear from this table that, for the pionts shown, theory and simulations agree to within the noise of the simulations. On the basis of this agreement, and the fact that the SPCE model was previously shown to accurately mimic t of real water along the liquid-vapor coexistence curve up to the critical point (see Introduction), we decided to check the theory’s usefulness for the prediction of e under conditions relevant to hydrothermal fluids. Thus, we set out to predict e of water along the slightly supercritical 400 O C isotherm over the pressure range 250 IP I 5000 bar. Both an equation of state” (needed to obtain densities) and experimental e’s for comparison4d are available for these conditions. Our detailed results are given in Table 3, and the e’s are graphed in Figure 1. It is clear from Figure 1 that the SPCE model together with the theory quantitatively reproduces the experimental e’s of real water along the 400 OC isotherm up to 3000 bar. Beyond this pressure the theory begins to overestimate t. This suggests that the model and theory can be used to quantitatively predict e for real water, provided y I2.4 (see Table 3). This is essentially the range of validity of the theory found previouslylo in its applications to the S P C model. The range y 5 2.4 permits an application to a large part of the P-T surface that is of geochemical interest. For example, along the 1000 OC isotherm, with P in the range 0 < P I 16 000 bar, y will always be less than 2, for which the model and theory are
I 1 2 2 3 3 4 4 4 5 5 5 6 6 6 6 7 7 7 7
n 0 0 2 0 2 0 2 4 0 2 4 0 2 4 6 0 2 4 6
Qa
0.196169 1+1 0 -0.5656255+0 -0.1992994+0 -0.2729019+0 -0.7211817-1 -0.6515934-1 0.4309885-1 -0.6231197-2 0 0.2606697-1 0.4941233-2 0.6329071-2 0.924420&2 -0.3129177-2 0.2419791-2 0.2402589-2 0.1733714-2 -0.2210061-2
I
n
Oina
8 8 8 8 8 9 9 9 9 9 10 10 10 10 10 10
0 2 4 6 8 0 2 4 6 8 0 2 4 6 8 10
0.3875601-3 0.2779196-3 -0.1457423-3 -0.9730021-3 0.2220563-3 -0.1003775-3 -0.1343313-3 -0,2265286-3 -0.279633G-3 0.1765097-3 -0,75240114 -0,77539874 0.758 16674 -0.36857214 0.90239894 -0.15548394
O,,,
a For this model = 0 for odd n and &, E Q I ~ / ( s u ~ ~ +where ~ ) ’ / ~a and c are the Lennard-Jones core parameters in the SPCE model.* a/k = 78.23 K, c = 3.166 X 10-8 cm. The notation +n(-n) means that the preceding number is to be multiplied by lP(lO+); e.g., 0.1961691+1 means 1.961691.
TABLE 2 Values of t for SPCE Water from Theory and Simulations at High Temperatures and Pressures T (K) P (kbar) density (g/cm3) @(theory) $(simulations) 809.8 1091.2 772.2 781.6 1040.3 1278.1
0.5535 5.617 5.670 9.796 15.097 19.950
0.2570 0.7016 0.8760 1.ooo 1.000 1.om
3.76 8.35 18.70 23.37 14.88 11.06
* * *
3.88 0.36 8.76 f 1.02 18.68 2.30 21.06 & 2.24 13.30 1.68 9.97 1.22
a Obtained using eqs 1 and 5-7. From the molecular dynamics simulationsin ref 7. The uncertaintiesrepresent two standard derivations. Complete details are given in ref 7.
TABLE 3 400 OC p (kbar) 0.25 0.35 0.50 1.00 1.50 2.00 2.40 2.80 3.00 3.50 4.00 4.50 5.00
Theoretical and Experimental E’S for Water at pa
9
0.17609 0.502 80 0.613 12 0.735 11 0.796 63 0.840 13 0.868 77 0.892 10 0.905 77 0.931 70 0.954 96 0.975 78 0.99506
0.460586 1.315 14 1.60369 1.92278 2.083 69 2.197 47 2.272 38 2.333 40 2.369 16 2.436 98 2.497 82 2.552 28 2.60271
- G A / Y ~ ‘ G E / Y ~ ‘e(theoryId 0.038 41 0.921 16 2.98 0.017 54 0.913 70 9.43 0.01607 0.914 11 12.26 0.015 26 0.91978 15.82 0.015 00 0.926 81 17.88 0.014 96 0.933 55 19.49 0.014 93 0.938 38 20.61 0.014 95 0.942 29 21.57 0.014 94 0.945 06 22.16 0.015 01 0.950 77 23.34 0.015 04 0.955 59 24.43 0.015 07 0.959 66 25.44 0.015 20 0.96280 26.36
c(experiment) 2.7e 9.3e 12.5’ 15.48f 17.66‘ 19.4Of 20.59 21.5Of 22.le 22.9e 233 24.lC 24.6e
a ;(the reduced density) E pd(SPCE) = 1.060764d, where d is the density in g/cm3. For P -< 2.8 kbar, densities were taken from ref 6 . For P>-2.8 kbar weused theHaarequationofstate(ref 17) forthecalculation of p. b y (4*/9)pZp/kT; p(SPCE) = 2.35 D. From the spherical multipolemoments in Table 1 and the computational proceduresdescribed in refs 10and 11. All the sphericalmultipolemoments up to and including 1 = 10 were included in the spherical harmonic expansion of the pair potential.lOJ1 From eqs 1 and 5 and the entries for G”/y2 and G8/y2 in this table. e From ref 4. f From ref 6. quantitatively accurate. It is also seen from Figure 1 that the S P C E model provides a much more accurate prediction of e than does its predecessor, the SPC model. The effective dipole moment of the SPC model was clearly too small for use in these applications. W e conclude with a comment on the computer time required for these calculations. The rate-determining step is the time required to obtain the underlying Lennard-Jones pair distribution
The Journal of Physical Chemistry,Vol. 98, No.24, 1994 6233
Dielectric Constants at High Temperatures and Pressures
t
I
I
25ti.
20
References and Notes
lo.
5t O 1
0
Acknowledgment. We are grateful to the Natural Sciences and Engineering Research Council of Canada and the Royal Society for financial assistance.
1
2
3
4
5
I 6
p
(Kbars) Figure 1. Experimental and calculated values of the static dielectric constant of water at 400 O C . Points are the experimental values for real water (listed in column 7, Table 3). The curves were calculated using eq 1 and 5-7. The solid curve is for the SPCE model; the dashed curve is for the SPC model. (The slight irregularity in the curves at 3 kbar is due to our switching from the densities listed in ref 6 to those obtained using the equation of state in ref 17. It is irrelevant to our conclusions.)
functions (go(r)). For state conditions in which a previously published fit to MD go(r)’sl*can be used (Le.,0.5 IkT/e I 5; 0.35 Ipa3 Il . l ) , the calculation of e required about 1 min of CPU time per statecondition (calculationsperformedon a Silicon Graphics Challenge XL workstation) and is thus some 4 orders of magnitude faster than a direct simulation for e of the SPCE model. If go(r)has to be simulated, there will still be a computer time saving of 1-2 orders of magnitude, depending on the temperature and the accuracy sought. This saving arises because a simulation for go(r) of a Lennard-Jones fluid is much faster than one for e of a model such as SPCE water. Another option for obtaining go(r) may be to use one of the atomic liquid-state perturbation theories’9*20for it.
(1) Friedman, H. L.; Dale, W. D. T. In Statistical Mechanics, Part A: Equilibrium Techniques; Berne, B. J., Ed.; Plenum: New York, 1977; pp 85-135. (2) Helgtson, H. C.; Kirkham, D. H.; Flowers, G. C. Theoretical Predictionofthe ThermodynamicBehaviourof Aqueous Electrolytes at High Pressures and Temperatures: 1. Summary of the Thermodynamic/ Electrostatic Properties of the Solvent. A.J.S. 1974, 274, 1089-1 198. (3) McKenzie,W.F.; Helgeson,H. C. Geochem. Cosmochim.Acta 1984, 48, 2167-2177 and references therein. (4) Heger, K.;Uematsu, M.; Franck, E. U. Ber. Bunsen-Ges.Phys. Chem. 1980,84,758-762. (5) Uematsu, M.; Franck, E. U. J . Phys. Chem. ReJ Data 1980,9,12911306. (6) Deul, R.; Franck, E. U. Ber. Bunsen-Ges. Phys. Chem. 1991, 95, 847-853. (7) Waasennan, E. A.; Wood, B. J. Molecular Dynamics Study of the Dielectric Constant of Water Under High Temperature and Pressure Conditions Ber. Bunsen-Ges. Phys. Chem., in press. (8) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J . Phys. Chem. t!7,91,626p-4271>
(9) Guiasani, Y.; Guillot, B. J . Chem. Phys. 1993, 98, 8221-8235. (10) Goldman, S.; Joslin, C. J . Chem. Phys. 1993, 99, 3021-3029. (11) Goldman, S.; Joslin, C. J . Phys. Chem. 1993, 97, 12349-12355. (12) Berendsen, H. J.C.;Postma, J. P. M.;vanGunstern, W. F. Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981; p 331. (13) Goldman, S . Mol. Phys. 1990, 71, 491-507. (14) Tani, A.; Henderson, D.; Barker, J. A.; Hecht, C. E. Mol. Phys. 1983,48,863-869. (15) Gray, C. G.; Gubbins, K. E. Mol. Phys. 1975, 30, 1481-1487. (16) Hernando, J. A,; Gamba, Z. J . Chem. Phys. 1992,97, 5142-5147. (17) Haar, L.;Gallagher, J. S.; Kell, G. S.NBSINRC Steam Tables; Hemisphere: Bristol, PA, 1984. (18) Goldman, S.J . Phys. Chem. 1979,83, 3033-3037. (19) Barker, J. A.; Henderson, D. J. Chem. Phys. 1967,47,4714-4721. (20) Weeks,J. D.; Chandler, D.; Andersen, H. C. J. Chem. Phys. 1971, 54, 5237-5247.
