Fluid Phase Equilibria 193 (2002) 277–287
Theory and computer simulation of the coordination number of square-well fluids of variable width J. Largo, J.R. Solana∗ Departamento de F´ısica Aplicada, Universidad de Cantabria, 39005 Santander, Spain Received 3 September 2001; accepted 31 October 2001
Abstract Computer simulations have been performed to determine the coordination number of square-well (SW) fluids for wide ranges of densities, temperatures, and potential widths. These data are used to test the performance of a new coordination number model derived on the basis of the quasi-chemical approximation. The model presents a considerable improvement over other previously reported theoretical models based on similar grounds. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Theory; Computer simulation; Coordination number model; Quasi-chemical approximation; Statistical mechanics
1. Introduction In a square-well (SW) fluid, particles interact by means of a potential of the form ∞, if r < σ u(r) = −ε, if σ ≤ r ≤ λσ 0, if r > λσ
(1)
where σ is the diameter of the particles, −ε the potential depth, and λσ the range or potential width. This is perhaps the simplest model of a fluid that includes attractive as well as repulsive contributions to the potential. Therefore, SW potential is the simplest potential model whose shape resembles that of the potential of a real fluid with spherically symmetrical interactions. In spite of the fact of being excessively simple, SW fluid has been profusely studied from a theoretical viewpoint as well as from computer simulations. This is because it can be considered as the first stage in the search of a theory for the thermodynamic properties of real fluids, whereas its theoretical treatment is generally easier. Among the most successful theories for this kind of fluids is the Barker–Henderson perturbation theory [1] which provides fairly close agreement with the simulation data for the equation of state, although the predicted values for the internal energy are less accurate. Another class of theories frequently used ∗
Corresponding author. Tel.: +34-942-20-1447; fax: +34-942-20-1402. E-mail address:
[email protected] (J.R. Solana). 0378-3812/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 3 8 1 2 ( 0 1 ) 0 0 7 4 4 - 0
278
J. Largo, J.R. Solana / Fluid Phase Equilibria 193 (2002) 277–287
for SW fluids is that of integral equation theories [2], although they have the drawbacks of their greater complexity and, in most cases, their non-analytical character, apart from the fact that their performance is not always enough satisfactory. Finally, we will mention here the coordination number theories or generalized van der Waals theories [3–9]. This group of theories have in common their considerable simplicity, both conceptually and operationally. However, the theoretical models existing for the coordination number do not provide, in general, accurate results so that parametrizations, often lacking of a theoretical basis, of the simulation data are frequently used instead. Regarding the study of SW fluids by means of simulation, the amount of data available is considerable for the range λ = 1.5, whereas for other ranges the data are more scarce. The paper is organized as follows. Section 2 summarizes the existing models for the coordination number of SW fluids based on the quasi-chemical approximation. The new model is derived in Section 3. Finally, in Section 4 all these models are compared with the new simulation data, which are also reported in this section, and results are discussed. 2. Coordination number models for SW fluids Several theoretically-based coordination number models for SW fluids have been derived on the basis of the quasi-chemical approximation. Thus, Lee et al. [5] obtain Nc =
zV0 eε/2kT , V + V0 (eε/2kT − 1)
(2)
where Nc is the coordination number, that is, the average number of particles within the potential well, z is the lattice coordination number, that is, the maximum number of particles that can lie within the potential well, and V and V0 the actual en close-packed volumes, respectively, of the system. In terms of the reduced number density ρ ∗ = ρσ 3 , the preceding expression reads zρ ∗ eε/2kT Nc = √ , 2 + ρ ∗ (eε/2kT − 1)
(3)
zρ ∗ eε/αkT Nc = √ , 2 + ρ ∗ (eε/αkT − 1)
(4)
zρ ∗ eε/kT . ∗ ρmax + ρ ∗ (eε/kT − 1)
(5)
√ √ √ √ For a face centered cubic (fcc) lattice [10], z = 12 for λ < 2, 18 for 2 ≤ λ < 3, 42 for 3 ≤ λ < 2, and so on. However, for consistency a slightly modified expression has been proposed [6]
√ ∗ ∗ ∗ /(ρmax − ρ ∗ ), with ρmax = 2 being the density corresponding to V0 . Eq. (4) behaves where α = ρmax more satisfactorily than Eq. (3) both in the low density limit, where it is compatible with the right second ∗ virial coefficient of the SW fluid, and in the high density limit ρ ∗ = ρmax , in which Nc reaches its maximum value. A slightly different expression has been derived by Heyes [9] Nc =
J. Largo, J.R. Solana / Fluid Phase Equilibria 193 (2002) 277–287
279
The main differences between this expression and Eq. (3), apart from the factor 2 in the denominator ∗ of the exponent, arise in the values of ρmax and z. For the first of these quantities, the value 1 was taken [9], very close to the value 0.9428 which corresponds to the density of the hard-sphere (HS) fluid in equilibrium with the corresponding solid. To determine z use is made of the relationship between the coordination number and the radial distribution function (r.d.f.) g(r) λ ∗ Nc = 4πρ g(r ∗ )r ∗2 dr ∗ , (6) 1
where r ∗ = r/σ . In the limit ρ → 0 g(r) = e−βu(r) ,
(7)
and Eq. (6) gives Nc = 43 πρ ∗ eε/kT (λ3 − 1).
(8)
On the other hand, in the limit ρ ∗ → 0, expanding Eq. (5) in power series of the reduced density ρ ∗ and retaining only up to the first-order term Nc = zρ ∗ eε/kT .
(9)
Equating both results we obtain z = 43 π(λ3 − 1).
(10)
Assuming that this result holds at any density and introducing it into Eq. (5) one finally obtains [9] Nc =
(4/3)πρ ∗ (λ3 −1)eε/kT . 1+ρ ∗ (eε/kT −1)
(11)
This expression provides the correct low density limit of the coordination number. 3. A new model for the coordination number of SW fluids In the high temperature limit T → ∞, Eq. (5) reduces to Nc =
ρ∗ z. ∗ ρmax
(12)
On the other hand, at infinity temperature, the SW fluid behaves as a HS fluid, for which the average number NcHS of particles up to a distance λ of a given particle is given by Eq. (6) by taking for the radial distribution function g(r) that corresponding to a HS fluid g HS (r). Therefore, we can write z=
∗ ρmax N HS . ρ∗ c
(13)
For the HS fluid r.d.f., we can take the Percus–Yevick solution which has been obtained in analytical form by a number of authors. We will consider here the formulation of Chang and Sandler [11] from which we only need the contribution corresponding to the first coordination shell, corresponding to the range 1 ≤ r ∗ ≤ 2, since we are interested in potential ranges 1 ≤ λ ≤ 2. In the same paper the authors
280
J. Largo, J.R. Solana / Fluid Phase Equilibria 193 (2002) 277–287
Table 1 Coordination number Nc for square-well fluids of variable width from MC calculations λ
ρ∗
T ∗ = kT/ 0.5
0.7
1.0
1.5
2.0
3.0
5.0
1.97 2.65 3.46 4.05 4.66 5.29 5.96
1.13 1.69 2.27 2.89 3.55 4.26 5.09
0.81 1.26 1.76 2.33 2.96 3.70 4.55
0.62 1.00 1.44 1.96 2.57 3.29 4.16
0.54 0.89 1.30 1.79 2.38 3.10 3.97
0.47 0.79 1.17 1.63 2.21 2.91 3.78
0.42 0.71 1.07 1.51 2.07 2.76 3.63
2.37 3.31 4.09 4.90 5.59 6.39 7.34
1.58 2.35 3.12 3.92 4.79 5.75 6.85
1.23 1.90 2.63 3.42 4.32 5.34 6.51
1.09 1.72 2.43 3.22 4.12 5.16 6.35
0.97 1.56 2.24 3.02 3.93 4.98 6.19
0.88 1.44 2.10 2.87 3.78 4.84 6.06
1.1
0.20 0.30 0.40 0.50 0.60 0.70 0.80
1.2
0.20 0.30 0.40 0.50 0.60 0.70 0.80
1.3
0.20 0.30 0.40 0.50 0.60 0.70 0.80
2.47 3.50 4.48 5.42 6.41 7.55 8.84
1.89 2.82 3.78 4.80 5.92 7.15 8.49
1.68 2.57 3.53 4.56 5.71 6.97 8.33
1.50 2.36 3.30 4.35 5.51 6.79 8.16
1.38 2.21 3.14 4.19 5.37 6.65 8.03
1.4
0.20 0.30 0.40 0.50 0.60 0.70 0.80
3.75 4.99 6.04 6.96 8.00 9.27 10.64
2.62 3.80 4.96 6.18 7.48 8.88 10.27
2.33 3.47 4.66 5.92 7.27 8.68 10.09
2.09 3.21 4.41 5.69 7.07 8.50 9.90
1.93 3.03 4.22 5.52 6.91 8.35 9.75
1.5
0.20 0.30 0.40 0.50 0.60 0.70 0.80
3.44 4.89 6.22 7.59 9.06 10.56 11.89
3.04 4.45 5.85 7.32 8.85 10.36 11.71
2.73 4.12 5.56 7.08 8.64 10.15 11.53
2.53 3.91 5.37 6.90 8.47 9.99 11.38
1.6
0.20 0.30 0.40 0.50 0.60 0.70 0.80
4.50 6.16 7.57 9.08 10.71 12.29 13.67
3.84 5.51 7.13 8.80 10.49 12.08 13.47
3.45 5.12 6.81 8.55 10.27 11.87 13.28
3.20 4.87 6.59 8.35 10.08 11.70 13.11
J. Largo, J.R. Solana / Fluid Phase Equilibria 193 (2002) 277–287
281
Table 1 (Continued ) λ
ρ∗
T ∗ = kT/ 0.5
0.7
1.0
1.5
2.0
3.0
5.0
6.38 7.87 9.22 10.71 12.54 14.30 15.89
4.79 6.70 8.52 10.40 12.30 14.06 15.64
4.23 6.20 8.16 10.14 12.05 13.81 15.39
3.94 5.92 7.93 9.93 11.85 13.61 15.18
1.7
0.20 0.30 0.40 0.50 0.60 0.70 0.80
1.8
0.20 0.30 0.40 0.50 0.60 0.70 0.80
5.98 8.12 10.08 12.22 14.38 16.44 18.38
5.13 7.41 9.67 11.92 14.10 16.14 18.06
4.76 7.08 9.40 11.69 13.87 15.90 17.80
1.9
0.20 0.30 0.40 0.50 0.60 0.70 0.80
8.38 10.20 12.05 14.33 16.88 19.43 22.04
6.13 8.78 11.37 13.96 16.53 19.04 21.60
5.67 8.37 11.06 13.70 16.25 18.74 21.26
2.0
0.20 0.30 0.40 0.50 0.60 0.70 0.80
7.37 10.42 13.32 16.40 19.51 22.71 26.08
6.69 9.83 12.95 16.06 19.16 22.31 25.63
reported also an analytical expression for the coordination number NcHS of the HS fluid up to a distance λ, which is reproduced in the Appendix A for convenience. √ ∗ Regarding the behavior of ρmax , the maximum density a HS system can reach is ρ ∗ = 2, which corresponds to a volume V0 , the regular close packing volume. This upper limit holds also for SW systems. However, for more irregularly packed systems, other limiting values of the density of the HS system are possible [12], and the same must be true for a SW system. Let us consider the infinity range limit λ → ∞ of Eq. (5), with z given by Eq. (13). In this limit, particles move in an uniform background and the potential well has no effect on the coordination number, that is, Nc (λ = ∞) = NcHS = N. From ∗ Eqs. (5) and (13) this requires that ρmax = ρ ∗ in this limit. In the opposite limit λ = 1 the width of the potential well is zero, but particles can stick to other particles so that Nc√ = NcHS except in the limit ∗ ∗ . The last situation requires that ρmax = 2 for λ = 1. Perhaps the T → ∞ or in the limit ρ ∗ = ρmax
282
J. Largo, J.R. Solana / Fluid Phase Equilibria 193 (2002) 277–287
Fig. 1. Coordination number for a square-well fluid with range λ = 1.1 as a function of the reduced density ρ ∗ = N σ 3 /V for several values of the reduced temperature T ∗ = kT/. Points: simulation data of the present paper. Crosses: simulation data from [9]. Continuous line: Eq. (15); dash-dotted line: Eq. (3); dashed line: Eq. (4); dotted line: Eq. (11).
J. Largo, J.R. Solana / Fluid Phase Equilibria 193 (2002) 277–287 ∗ simplest way of interpolating the values of ρmax between the two extreme values of λ is to take 1 √ ∗ ρmax = ρ ∗ + 3 ( 2 − ρ ∗ ). λ
283
(14)
From Eqs. (5) and (13), the expression of the coordination number we propose is Nc =
∗ ρmax NcHS eε/kT , ∗ ρmax + ρ ∗ (eε/kT − 1)
(15)
∗ given by Eq. (14). with ρmax
4. Results and discussion Computer simulations of the coordination number of SW fluids have been reported by several authors [5,8,9] for values of λ = 1.1, 1.2, and 1.5. Here we present more extensive simulations for values of λ ranging from 1.1 to 2.0. Monte Carlo simulation have been performed in the NVT ensemble for a system of 512 SW particles initially placed in a regular configuration. System was allowed to equilibrate
Fig. 2. As in Fig. 1 for λ = 1.5, except for the fact that here the crosses (for T ∗ = 2.0) represent the simulation data from [5]. The crosses are hardly distinguishable from our own data at the scale of the figure.
284
J. Largo, J.R. Solana / Fluid Phase Equilibria 193 (2002) 277–287
Fig. 3. As in Fig. 1 for λ = 1.7, except for the lack of simulation data from other authors.
for 2 × 104 cycles, the first 20% of them at an extremely high temperature, after which averages were taken over 5 × 104 additional cycles, each cycle consisting of an attempt move per particle. Results are listed in Table 1. They are in agreement with those previously reported in the cases where comparison is possible. In Figs. 1–3 we compare the results from Eqs. (3), (4), (11) and (15) with the simulation data for λ = 1.1, 1.5, and 1.7. One can see that Eqs. (3) and (4) are in reasonable agreement with simulation data for λ = 1.5, whereas the values they predict are too high for λ = 1.1 and too low for λ = 1.7. Similar results are obtained for other values of λ. The reason for this behavior is the discontinuous variations of z with λ in these models, as mentioned in Section 2. Eq. (11) is exact at low densities for any value of λ, but at high densities agreement with simulation is unsatisfactory for ranges λ ≤ 1.7. The expression (15) derived in this paper gives excellent results in all cases, except for temperatures T ∗ ≤ 0.7 and λ ≤ 1.2. The averaged relative deviations of Eqs. (3), (4), (11) and (15), with respect to the whole set n of simulation data listed in Table 1, calculated as (1/n)Σi=1 |Nccalc − Ncsim |/Ncsim , where n is the number of simulation data, are 45, 52, 25, and 4%, respectively. This shows clearly the excellent performance of Eq. (15) developed in this paper From the preceding considerations, we can conclude that the expression derived in this paper for the coordination number of SW fluids provides a remarkable improvement over the previously reported theoretically-based models. However, further improvements are needed before we can obtain accurate
J. Largo, J.R. Solana / Fluid Phase Equilibria 193 (2002) 277–287
285
results for the thermodynamic properties of SW fluids on the basis of this model, because these properties are very sensitive to small deviations of the behavior of the coordination number model with respect to simulation data. One possible improvement would consist in using for the HS r.d.f. a more accurate expression than provides the Percus–Yevick theory, such as those derived in [13–15]. However, this results in more complicate expressions for the coordination number without a significant increase in the accuracy of the results. This is because the calculated values of the coordination number of the HS fluid are relatively insensitive to accuracy of the r.d.f. used. List of symbols A, B, C functions implicit in Eq. (A.1) N number of particles in the system Nc coordination number T temperature T ∗ = kT/ reduced temperature V volume V0 regular close packing volume a1 , a2 , a3 functions implicit in Eq. (A.1) f function implicit in Eq. (A.1) k Boltzman’s constant g(r) radial distribution function r radial distance r ∗ = r/σ reduced distance u potential energy function u∗ u/ y+ , y− functions implicit in Eq. (A.1) z coordination number of the lattice z d , zs functions implicit in (A.1) Greek symbols β 1/kT depth of the potential well γ λ−1 η packing fraction λ radial distance or range of the potential in units of σ ρ = N/V number density ∗ 3 ρ = ρσ reduced density ∗ ρmax maximum (reduced) density of the system σ diameter of a hard-sphere
Acknowledgements We are grateful to the Spanish ‘Dirección General de Investigación’ (DGI) for the financial support under Grant No. BFM2000-0014.
286
J. Largo, J.R. Solana / Fluid Phase Equilibria 193 (2002) 277–287
Appendix A The coordination number of a HS fluid up to a distance λ in the Percus–Yevick approximation is given [11] NcHS = 4πρ ∗ (a1 t1 + a2 t2 + a3 t3 ),
(A.1)
where t1 =
exp(Aγ )(λA − 1) − A + 1 , A2
t2 =
λexp(Bγ )[B cos (Cγ ) + C sin (Cγ )] − B B2 + C2 −
t3 =
exp(Bγ )[(B 2 − C 2 ) cos (Cγ ) + 2BC sin (Cγ )] − (B 2 − C 2 ) , (B 2 + C 2 )2
(A.2)
(A.3)
λexp(Bγ )[B sin (Cγ ) − C cos (Cγ )] + C B2 + C2 exp(Bγ )[(B 2 − C 2 ) sin (Cγ ) − 2BC cos (Cγ )] + 2BC (B 2 + C 2 )2
(A.4)
a1 =
−2η(1 − η − 3η2 ) + (1 − 3η − 4η2 )zd + (1 + η/2)zd2 , 3(2η2 + zd2 )(1 − η)2
(A.5)
a2 =
η(2 + 4η − 3η2 ) − (1 − 3η − 4η2 )zd + 2(1 + η/2)zd2 , 3(2η2 + zd2 )(1 − η)2
(A.6)
a3 =
(1 − 3η − 4η2 )(4η2 + zd2 ) + η(2 − 5η2 )zd . √ 3zs (2η2 + zd2 )(1 − η)2
(A.7)
−
and
with γ = λ − 1, A=
−2η + zd , 1−η
−2η − (1/2)zd 1−η √ 3zs C= , 2(1 − η) B=
(A.8) (A.9) (A.10) (A.11)
zd = y+ − y− ,
(A.12)
zs = y+ + y− ,
(A.13)
J. Largo, J.R. Solana / Fluid Phase Equilibria 193 (2002) 277–287
y± = (2ηf )
1/3
1/2
2η4 +1 f2
287
1/3 ±1
f = 3 + 3η − η2 ,
,
(A.14) (A.15)
and η = (π/6)ρ ∗ , the packing fraction. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
J.A. Barker, D. Henderson, J. Chem. Phys. 47 (1967) 2856. C. Caccamo, Phys. Rep. 274 (1996) 1, and references therein. J.H. Vera, J.M. Prausnitz, Chem. Eng. J. 3 (1972) 1. S.I. Sandler, Fluid Phase Equilib. 19 (1985) 233. K.-H. Lee, M. Lombardo, S.I. Sandler, Fluid Phase Equilib. 21 (1985) 177. K.-H. Lee, S.I. Sandler, N.C. Patel, Fluid Phase Equilib. 25 (1986) 31. R.J. Lee, J.C. Chao, Mol. Phys. 65 (1988) 1253. M. Guo, W. Wang, H. Lu, Fluid Phase Equilib. 60 (1990) 37. D.M. Heyes, J. Chem. Soc., Faraday Trans. 87 (1991) 3373. J.O. Hirschfelder, C.F. Curtiss, R.B. Bird, Molecular Theory of Gases and Liquids, Wiley, New York, 1964, p. 1037. J. Chang, S.I. Sandler, Mol. Phys. 81 (1994) 735. L.V. Woodcock, Ann. New York Acad. Sci. 371 (1981) 274. Y. Tang, B.C.-Y. Lu, J. Chem. Phys. 103 (1995) 7463. S. Bravo Yuste, A. Santos, Phys. Rev. A 43 (1991) 5418. S. Bravo Yuste, M. López de Haro, A. Santos, Phys. Rev. E 53 (1996) 4820.