- Email: [email protected]
Calculation of the chemical potentials of hard-sphere solutes solvated in hard-sphere solvents using the radial free-space distribution function
5 March 1999
Chemical Physics Letters 301 Ž1999. 481–486
Calculation of the chemical potentials of hard-sphere solutes solvated in hard-sphere solvents using the radial free-space distribution function Byoung Jip Yoon ) , Kyoung K. Baeck, Sang Il Jeon Department of Chemistry, Kangnung National UniÕersity, Kangnung, Kangwondo 210-702, South Korea Received 5 October 1998; in final form 24 November 1998
Abstract The chemical potentials of hard-sphere solutes solvated in hard-sphere solvents were calculated by a Monte Carlo method using the radial free-space distribution function. This method is based on calculating the entropy by comparing the free volume of a molecule with that of an ideal gas, and is applicable even when the diameter of solute is very large and the density of system is high. The results are in good agreement with previous simulations using the insertion method Žfor small solutes. and analytical expressions derived from Boublic–Mansoori–Carnahan–Starling–Leland equation of state. q 1999 Elsevier Science B.V. All rights reserved.
1. Introduction The chemical potential is one of the most important thermodynamic properties since it determines the directions of spontaneous processes. However, its theoretical estimation cannot simply be made as an ensemble average w1x. The insertion method w2–4x, i.e., calculating the probability of inserting a particle into an equilibrated ensemble has been considered as a practical and exact method that has been used in Monte Carlo ŽMC. simulations. The chemical potentials of hard-sphere solutes in hard-sphere solvents have been calculated w5,6x by the insertion method and compared with the analytical expression derived from the equation of state of Boublic–Mansoori–
) Corresponding author. Fax: q82 391 647 1183; e-mail: [email protected].
Carnahan–Starling–Leland ŽBMCSL. w7,8x. The insertion method for calculating the chemical potential has a limit when the density of a system is high and the sizes of solutes are very large. This limit is reached when the chemical potentials are greater than about 14 kT, as this corresponds to insertion probabilities that are less than 10y6 , i.e., the lowest probability detected when the usual number of simulation trial steps are of the order of millions. In this work, we use the different method employing the radial free-space distribution function ŽRFSDF. which gives information regarding the probability that a molecule can move freely Žin the cavity formed by surrounding nearest-neighbor molecules.. The method employing RFSDF overcomes the difficulties that arise when the density is high and the size of solute is very large. This method has provided good results for estimating the free energies of non-solvated single component systems
0009-2614r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 Ž 9 9 . 0 0 0 7 8 - 0
B.J. Yoon et al.r Chemical Physics Letters 301 (1999) 481–486
482
with hard-sphere and Lennard-Jones potentials w9,10x. In this work, the method is applied to solvated systems with a hard-sphere potential, and the results have been compared with the data of insertion method w6x and BMCSL predictions w7,8x.
as half of the average nearest-neighbor distance. The definition of the cut-off distance means that the free volume for an ideal gas is simply Vidf s 4r3p Ž r f . 3 because RFSDF for an ideal gas is unity at all values of r. The cut-off distance for a hard-sphere fluid solvent, r 0f , has been defined w9,11x as follows:
2. Method
r 0f s
The RFSDF has been defined in order to calculate the entropy andror chemical potential of fluids and solids by MC simulation method w9x. The entropy is easily calculated from the RFSDF without any of the difficulties which might arise in the usual method 1. The RFSDF, z Ž r ., is given in the MC procedure by the following ratio:
z Ž r. s
Acceptances of displacement of r Trials of displacement of r
.
Ž 1.
z Ž r . starts from unity at r s 0 and decreases at large r. The RFSDF is the distribution function of the particle-cavity Žincluding the cavity formed by removing the particle. structure. In the case of hardsphere fluids and solids, the particle and cavity are interchangeable and thus the function is related to the usual two-cavity distribution function w11x, y Ž r ., i.e., yŽ r .
z Ž r. s
y Ž 0.
.
Ž 2.
R
s ln
Vf
ž / Vidf
Ž 3.
in which R is the gas constant, V f is the free volume of a molecule in the system, and Vidf is that of an ideal gas. The free volume V f is calculated from the RFSDF. f
V s4p
rf 2
H0
r z Ž r. dr ,
Ž 4.
where r f is a cut-off distance to specify the molecular volume, which has been formally defined w9,11x
1
See the references cited in Ref. w9x.
2 r 1r3s
,
Ž 5.
where r is the number density and s is the diameter of the hard sphere. This definition is equal to half of the nearest-neighbor distance assuming the structure to have a simple cubic lattice. In this work, it is not easy to define the cut-off distance for the solute molecule surrounded by solvent molecules, but we used the following mean value to evaluate the pressure in Eq. Ž8. below, rfs
Ž1qd. 2
r 0f ,
Ž 6.
where d is the ratio of the diameter of solute to that of solvent. To determine the excess chemical potential over an ideal gas, mex , then the entropy term is added to the pressure term,
mex
p ex V s
kT
The excess entropy over an ideal gas is calculated from the free volume as follows w9x, S ex
1
S ex
Ž 7.
y RT
R
in which k is the Boltzmann constant and T is temperature. The excess pressure, p ex , of solute over an ideal gas is calculated from the radial distribution function ŽRDF. as follows w12x: p ex V
2 1qd s
RT
3
ž
2
3
/
prs 3 g 10
ž
1qd 2
/
Ž 8.
In Eq. Ž8., g 10 is the radial distribution function of solvent molecules around a solute molecule, and g 10 ŽŽ1 q d .r2. is the value when the surfaces of solute and solvent molecule are in contact. For the entropy of solvent molecule, the RFSDF of solvent has been used and the cut-off distance given in Eq. Ž5. has been used. For the entropy of the solute molecule, the RFSDF of the solute has been calculated separately and the cut-off distance
B.J. Yoon et al.r Chemical Physics Letters 301 (1999) 481–486
483
Table 1 The excess pressures, entropies and chemical potentials of solvent and solute of hard-sphere system at different solute–solvent ratios. Small numbers are the standard deviationa d
Solvent p 0ex VrRT
Solute ex
yS rR
ex
m rkT
p ex VrRT
yS exrR
Ž10y3 .
Ž10y3 .
mexrkT
mexrkT b
mexrkT c
0.059 0.104 0.169 0.255 0.368 0.502 0.678 0.886 1.146 1.334 1.753 2.163 2.606 3.056 3.589 4.230 4.845 5.563
0.054 0.095 0.153 0.232 0.334 0.464 0.624 0.817 1.047 1.317 1.630 1.989 2.397 2.858 3.378 3.963 4.610 5.315
0.054 0.095 0.153 0.232 0.334 0.464 0.623 0.815 1.044 1.311 1.621 1.977 2.381 2.837 3.348 3.916 4.546 5.240
0.232 0.453 0.776 1.210 1.816 2.639 3.540 4.720 6.151 7.913 9.748 11.839 14.524 16.579 20.894 24.309 27.610
0.235 0.443 0.764 1.224 1.850 2.669 3.711 5.004 6.586 8.457 10.621 13.122 15.4 – – – –
0.235 0.442 0.762 1.219 1.839 2.648 3.669 4.929 6.452 8.264 10.389 12.853 15.680 18.897 22.527 26.597 31.130
0.580 1.348 2.617 4.397
0.543 1.175 2.309 4.096
0.543 1.171 2.293 4.062
3
rs s 0.1 Ž10y4 .
Ž10y4 .
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4
0.356 3 Žy. d 0.266 2 Žy. 0.266 1Žy. 0.2652 Žy. 0.268 3 Žy. 0.268 2 Žy. 0.2674 Žy. 0.270 2 Žy. 0.271 2 Žy. 0.271 3 Žy. 0.274 3 Žy. 0.274 2 Žy. 0.276 3 Žy. 0.278 3 Žy. 0.280 1Žy. 0.282 3 Žy. 0.2844 Žy. 0.2852 Žy.
0.239 8 0.239 1 0.23510 0.241 4 0.2396 0.239 14 0.2338 0.242 9 0.2438 0.244 7 0.246 3 0.25114 0.25112 0.2559 0.256 9 0.261 9 0.26710 0.2695
rs 3 s 0.4 Ž10y3 .
Ž10y4 .
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2
1.134 3 Žy. 1.072 7 Žy. 1.071 4 Žy. 1.074 7 Žy. 1.0794 Žy. 1.083 4 Žy. 1.092 4 Žy. 1.100 7 Žy. 1.110 8 Žy. 1.120 5 Žy. 1.134 14 Žy. 1.1535 Žy. 1.171 6 Žy. 1.194 11Žy. 1.219 7 Žy. 1.2455 Žy. 1.290 6 Žy.
1.516 2 1.5236 1.499 3 1.534 2 1.520 4 1.512 3 1.541 4 1.561 4 1.578 5 1.560 2 1.613 2 1.6555 1.658 3 1.7254 1.811 5 1.842 10 1.920 7
rs 3 s 0.8 Ž10y3 .
Ž10y4 .
0.0 0.2 0.4 0.6
3.316 9 Žy. 3.299 21Žy. 3.313 26 Žy. 3.32314 Žy.
6.6079 6.604 9 6.650 9 6.658 10
0.595 0.505 0.501 0.506 0.507 0.507 0.500 0.512 0.514 0.515 0.520 0.525 0.527 0.533 0.535 0.543 0.550 0.554
2.650 2.596 2.569 2.608 2.599 2.595 2.633 2.671 2.687 2.680 2.747 2.808 2.855 2.920 3.030 3.087 3.210
9.923 9.903 9.963 9.981
0.028 4 0.049 2 0.0794 0.120 4 0.1755 0.236 4 0.321 9 0.42117 0.54948 0.686 6 0.834 17 1.052 28 1.271 56 1.47711 1.740 91 2.083133 2.40583 2.758 97
0.032 3 Žy. d 0.054 2 Žy. 0.090 3 Žy. 0.1354 Žy. 0.1936 Žy. 0.266 5 Žy. 0.358 9 Žy. 0.46517 Žy. 0.59717 Žy. 0.748 27 Žy. 0.914 8 Žy. 1.111 25 Žy. 1.33547 Žy. 1.579 39 Žy. 1.849 9 Žy. 2.14765 Žy. 2.460 51Žy. 2.84580 Žy.
Ž10y2 .
Ž10y3 .
0.132 1 0.262 1 0.4472 0.6952 1.051 5 1.5596 2.089 1 2.801 6 3.67710 4.873 22 6.08517 7.51456 9.569 28 11.008 44 14.74965 17.43444 20.34771
0.100 2 Žy. 0.191 6 Žy. 0.328 3 Žy. 0.5158 Žy. 0.7659 Žy. 1.080 20 Žy. 1.451 37 Žy. 1.91950 Žy. 2.47578 Žy. 3.040 39 Žy. 3.662 73 Žq. 4.325163 Žq. 4.95555 Žq. 5.571148 Žq. 6.145125 Žq. 6.876 183 Žq. 7.293 213 Žq.
Ž10y2 .
Ž10y3 .
0.3571 0.850 2 1.6644 2.830 6
0.223 3 Žy. 0.498 6 Žy. 0.95316 Žy. 1.56737 Žy.
B.J. Yoon et al.r Chemical Physics Letters 301 (1999) 481–486
484 Table 1 Žcontinued. d
Solvent
Solute
p 0ex VrRT
ex
yS rR
ex
m rkT
p ex VrRT
yS exrR
Ž10y2 .
Ž10y3 .
4.538 20 6.91743 9.42317 13.452 50 17.670 50 25.336 102 31.523125 44.655169 63.985183
2.366 80 Žy. 3.342 50 Žy. 4.336 62 Žy. 5.445167 Žq. 6.268 93 Žq. 7.320 190 Žq. 7.971180 Žq. 8.884 102 Žq. 9.864 200 Žq.
mexrkT
mexrkT b
mexrkT c
6.904 10.259 13.759 18.897 23.938 32.686 43.494 53.539 73.849
6.676 10.157 15.4 – – – – – –
6.630 10.153 14.781 20.669 27.971 36.838 47.425 59.885 74.370
3
rs s 0.8 Ž10y3 .
Ž10y4 .
0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4
3.3396 Žy. 3.383 24 Žy. 3.399 10 Žy. 3.46311Žy. 3.520 15 Žy. 3.60716 Žy. 3.676 26 Žy. 3.803 24 Žy. 3.956 45 Žy.
6.656 8 6.736 8 6.852 9 7.0955 7.108 15 7.448 20 7.692 7 8.14510 8.702 20
9.995 10.119 10.251 10.558 10.706 11.055 11.368 11.948 12.658
a
The standard deviations of bunches of 2 = 10 5 samplings for solvent and 1 = 10 5 samplings for solute are given by the small numbers, and the orders of the standard deviations are noted on the top of each column. b Data from Ref. w6x. c Calculated from BMCSL equation. d The sign of the second coefficient Ž B . of the function yln z Ž r . s Ar q Br 3.
given in Eq. Ž6. has been used. For the pressure of solvent, p 0ex , the following equation has been used: p 0ex V RT
s
ž
Ny1 2 N
/
3
prs 3 g 00 Ž 1 . q
1 p ex V N RT
sity rs 3 s 0.8 we sampled twice as many times a for low densities.
Ž 9.
where N is the number of molecules used in the simulation and g 00 Ž1. is the RDF of solvent–solvent at contact and the last term is the contribution of solvent–solute interaction in Eq. Ž8.. The major difference between this method and the insertion method is that the insertion method does not disturb the solvent structure whereas this method involves equilibration of one solute and the solvent molecule. Thus the chemical potential from the insertion method is an infinitely diluted one while the chemical potential from this method is not infinitely diluted but has a finite density of solute, 1rŽnumber of solvent used in the simulation.. The structure of solvent near the solute is disturbed in this method and thus the pressure of system increases with an increase in the size of solute Žsee Table 1.. We have used 108 molecules Ž107 solvent molecules and one solute molecule. in MC calculations and have taken an average of 2 = 10 6 samplings Žtrials of move. for the solvent molecule and 10y6 samplings for the solute molecule at densities of rs 3 s 0.1 and 0.4 after discarding 5 = 10 5 configurations starting from the initial face-centered-cubic structure. For the den-
3. Result and discussion In this work, the excess entropy and the excess pressure over an ideal gas have been calculated for both solute and solvent. The RDFs of solvent molecules around the solute molecule at rs 3 s 0.4 for different sizes of solute are shown in Fig. 1 and compared with the analytic expression of Boublik w13x 2 . When the size of solute increases the first peak in the RDF increases, i.e. the pressure in Eq. Ž8. increases. In Fig. 2, we also show the RFSDFs of the solute molecule in the cavity surrounded by solvent molecules at rs 3 s 0.4. From Fig. 2 we can see that the RFSDF of the solute drops off from that of an ideal gas Ži.e., unity for all r . as the sizes of solute become large, and thus the free volume decreases. In Fig. 2, the RFSDFs using Eq. Ž2. and the analytic two-cavity function w13x derived for infinite solute dilution are compared. The RFSDFs from the simplified two-cavity function by Ben-Amotz w14x
2
The coefficient Ž p r16. in Eq. Ž51. of w13x must be Ž p r12..
B.J. Yoon et al.r Chemical Physics Letters 301 (1999) 481–486
Fig. 1. The radial distribution functions of solvent molecules around a solute molecule at rs 3 s 0.4 for different sizes of solute. The solid lines are simulated ones and the dotted lines are calculated using the analytical expression of Ref. w13x. s is the diameter of solvent molecule.
have been also considered and overlap with Boublik’s form except for d s 0, but are not shown in Fig. 2. However, the equations have a singular value for d s 0. The pressure and entropy calculated using
Fig. 2. The radial free-space distribution functions of a solute surrounded by solvent molecules at rs 3 s 0.4 for different sizes of solute. For an ideal gas, the function is unity for all r values. The solid lines are simulated results and those from the analytic equation of Ref. w13x are represented by dotted lines. For calculating RFSDF at ds 0.0 from the analytic equation, ds 0.001 has been substituted in order to avoid a singular value. s is the diameter of solvent molecule.
485
these two functions ŽRDF and RFSDF, respectively. at rs 3 s 0.1, 0.4, and 0.8 are given in Table 1 and compared with the results by the insertion method and BMCSL equation. The standard deviations for the entropy are smaller than those for the pressure since the entropy is calculated by the integration of RFSDF as in Eq. Ž4. while the pressure is calculated from the single-point value of RDF as in Eqs. Ž8. and Ž9.. The excess chemical potentials of solutes in Table 1 are plotted in Fig. 3; however, the standard deviations are not plotted since they are smaller than the size of symbols used in the plot. The present method generally gives the slightly higher chemical potentials than those of the insertion method and BMCSL analytic-equation at rs 3 s 0.1. This general phenomenon at low densities has appeared in previous calculations w9,10x. This might be due to the fact that the definition of cut-off distance assuming a simple cubic lattice in Eqs. Ž5. and Ž6. is not correct Ža bit too large. at low densities. We inserted the solute molecule at the beginning and therefore the pressure of a system for large solute molecules is a bit larger than that of the system used in the insertion method w6x. At the densities rs 3 s 0.4 and 0.8, the chemical potentials for large solutes of the present method are less than those of the insertion method and BMCSL estimations. It might
Fig. 3. The excess chemical potentials are compared with those of the insertion method and BMCSL predictions. The results of this work are denoted by diamonds while circles and lines denote the results of the insertion method and BMCSL predictions, respectively.
486
B.J. Yoon et al.r Chemical Physics Letters 301 (1999) 481–486
be also due to the fact that the cut-off distance in Eq. Ž6. using the arithmatic mean is not properly defined for large solute molecule. Another possibility is that the large solute at high density involves something like a fluid-to-solid phase transition w11x of the disturbed solvent molecules near the solute, and thus the free volume of the solute becomes large. Such a phase transition is plausible since the solute and solvent molecules are initially mixed in this work and, for large solute, the function yln z Ž r . is concave, i.e., when yln z Ž r . is fitted to yln z Ž r . s Ar q Br 3 , the coefficient B is positive; while for small solute the function is convex Žsee Table 1.. In a previous work w11x, it was found that, for the fluid phase, yln z Ž r . was convex and, for the solid phase, it was concave. The deviations of analytical forms for large solutes in RDF and RFSDF Žin Figs. 1 and 2. from simulated values could also be for the same reason. If these calculations are reasonable, the entropies of insertion method are smaller Ži.e., larger negative value. and not valid for large solute at high densities. The insertion probabilities are evaluated to be smaller than the true values because the simulation generates equilibrated ensembles of solvent molecules but does not generate the unequilibrated, large vacancies that are suitable for a large solute to occupy, and they are not found effectively by checking the insertion probability of solute in equilibrated ensembles. One of the main advantages of this method is that we can calculate the pressure and the entropy separately in one simulation. Another importance is that we can estimate the chemical potentials of large solutes that cannot be calculated by the insertion method. The chemical potentials at high densities can also be calculated by the present method, and the RFSDFs are easily calculated during MC simulations by keeping the trajectories such as in Eq. Ž1. without requiring any other extra runs. When a small number of molecules are used in the simulation, the effect of the solute size on pres-
sure could be larger than in the case in which an infinite number of molecules is used in the simulation. Consequently, we are studying the effect of the number of molecules in the simulations. The reliability of the fluid-to-solid phase transition of solvent molecules near solute is going to be checked by calculating the RFSDF of solvent near a solute separately from that far from the solute. The results of these further studies will be presented in the near future.
Acknowledgements This study is supported by the academic research fund of Ministry of Education, Republic of Korea ŽBSRI-97-3414.. One of the authors ŽBJY. appreciates the partial support of Professor Ben-Amotz when he stayed at Purdue University.
References w1x R. Balescu, Equilibrium and Non-Equilibrium Statistical Mechanics, John Wiley, New York, 1975, p. 237. w2x B. Widom, J. Chem. Phys. 39 Ž1963. 2808. w3x B. Widom, J. Phys. Chem. 86 Ž1982. 869. w4x M.P. Allen, D.J. Tildesley, Computer Simulation of Liquid, Oxford, New York, 1991. w5x L.E.S. de Souza, A. Stamatopoulou, D. Ben-Amotz, J. Chem. Phys. 100 Ž1994. 1456. w6x A. Stamatopoulou, L.E.S. de Souza, D. Ben-Amotz, J. Talbot, J. Chem. Phys. 102 Ž1995. 2109. w7x T. Boublik, J. Chem. Phys. 53 Ž1970. 471. w8x G.A. Mansoori, N.F. Carnahan, K.E. Starling, T.W. Leland Jr., J. Chem. Phys. 54 Ž1971. 1523. w9x B.J. Yoon, H.A. Scheraga, J. Mol. Struct. ŽTHEOCHEM. 199 Ž1989. 33. w10x B.J. Yoon, S.D. Hong, M.S. Jhon, H.A. Scheraga, Chem. Phys. Lett. 181 Ž1991. 73. w11x B.J. Yoon, Chem. Phys. Lett. 234 Ž1995. 35. w12x A. Rotenberg, J. Chem. Phys. 43 Ž1965. 4377. w13x T. Boublik, Mol. Phys. 59 Ž1986. 775. w14x D. Ben-Amotz, J. Phys. Chem. 97 Ž1993. 2314.
Chemical Physics Letters 301 Ž1999. 481–486
Calculation of the chemical potentials of hard-sphere solutes solvated in hard-sphere solvents using the radial free-space distribution function Byoung Jip Yoon ) , Kyoung K. Baeck, Sang Il Jeon Department of Chemistry, Kangnung National UniÕersity, Kangnung, Kangwondo 210-702, South Korea Received 5 October 1998; in final form 24 November 1998
Abstract The chemical potentials of hard-sphere solutes solvated in hard-sphere solvents were calculated by a Monte Carlo method using the radial free-space distribution function. This method is based on calculating the entropy by comparing the free volume of a molecule with that of an ideal gas, and is applicable even when the diameter of solute is very large and the density of system is high. The results are in good agreement with previous simulations using the insertion method Žfor small solutes. and analytical expressions derived from Boublic–Mansoori–Carnahan–Starling–Leland equation of state. q 1999 Elsevier Science B.V. All rights reserved.
1. Introduction The chemical potential is one of the most important thermodynamic properties since it determines the directions of spontaneous processes. However, its theoretical estimation cannot simply be made as an ensemble average w1x. The insertion method w2–4x, i.e., calculating the probability of inserting a particle into an equilibrated ensemble has been considered as a practical and exact method that has been used in Monte Carlo ŽMC. simulations. The chemical potentials of hard-sphere solutes in hard-sphere solvents have been calculated w5,6x by the insertion method and compared with the analytical expression derived from the equation of state of Boublic–Mansoori–
) Corresponding author. Fax: q82 391 647 1183; e-mail: [email protected].
Carnahan–Starling–Leland ŽBMCSL. w7,8x. The insertion method for calculating the chemical potential has a limit when the density of a system is high and the sizes of solutes are very large. This limit is reached when the chemical potentials are greater than about 14 kT, as this corresponds to insertion probabilities that are less than 10y6 , i.e., the lowest probability detected when the usual number of simulation trial steps are of the order of millions. In this work, we use the different method employing the radial free-space distribution function ŽRFSDF. which gives information regarding the probability that a molecule can move freely Žin the cavity formed by surrounding nearest-neighbor molecules.. The method employing RFSDF overcomes the difficulties that arise when the density is high and the size of solute is very large. This method has provided good results for estimating the free energies of non-solvated single component systems
0009-2614r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 Ž 9 9 . 0 0 0 7 8 - 0
B.J. Yoon et al.r Chemical Physics Letters 301 (1999) 481–486
482
with hard-sphere and Lennard-Jones potentials w9,10x. In this work, the method is applied to solvated systems with a hard-sphere potential, and the results have been compared with the data of insertion method w6x and BMCSL predictions w7,8x.
as half of the average nearest-neighbor distance. The definition of the cut-off distance means that the free volume for an ideal gas is simply Vidf s 4r3p Ž r f . 3 because RFSDF for an ideal gas is unity at all values of r. The cut-off distance for a hard-sphere fluid solvent, r 0f , has been defined w9,11x as follows:
2. Method
r 0f s
The RFSDF has been defined in order to calculate the entropy andror chemical potential of fluids and solids by MC simulation method w9x. The entropy is easily calculated from the RFSDF without any of the difficulties which might arise in the usual method 1. The RFSDF, z Ž r ., is given in the MC procedure by the following ratio:
z Ž r. s
Acceptances of displacement of r Trials of displacement of r
.
Ž 1.
z Ž r . starts from unity at r s 0 and decreases at large r. The RFSDF is the distribution function of the particle-cavity Žincluding the cavity formed by removing the particle. structure. In the case of hardsphere fluids and solids, the particle and cavity are interchangeable and thus the function is related to the usual two-cavity distribution function w11x, y Ž r ., i.e., yŽ r .
z Ž r. s
y Ž 0.
.
Ž 2.
R
s ln
Vf
ž / Vidf
Ž 3.
in which R is the gas constant, V f is the free volume of a molecule in the system, and Vidf is that of an ideal gas. The free volume V f is calculated from the RFSDF. f
V s4p
rf 2
H0
r z Ž r. dr ,
Ž 4.
where r f is a cut-off distance to specify the molecular volume, which has been formally defined w9,11x
1
See the references cited in Ref. w9x.
2 r 1r3s
,
Ž 5.
where r is the number density and s is the diameter of the hard sphere. This definition is equal to half of the nearest-neighbor distance assuming the structure to have a simple cubic lattice. In this work, it is not easy to define the cut-off distance for the solute molecule surrounded by solvent molecules, but we used the following mean value to evaluate the pressure in Eq. Ž8. below, rfs
Ž1qd. 2
r 0f ,
Ž 6.
where d is the ratio of the diameter of solute to that of solvent. To determine the excess chemical potential over an ideal gas, mex , then the entropy term is added to the pressure term,
mex
p ex V s
kT
The excess entropy over an ideal gas is calculated from the free volume as follows w9x, S ex
1
S ex
Ž 7.
y RT
R
in which k is the Boltzmann constant and T is temperature. The excess pressure, p ex , of solute over an ideal gas is calculated from the radial distribution function ŽRDF. as follows w12x: p ex V
2 1qd s
RT
3
ž
2
3
/
prs 3 g 10
ž
1qd 2
/
Ž 8.
In Eq. Ž8., g 10 is the radial distribution function of solvent molecules around a solute molecule, and g 10 ŽŽ1 q d .r2. is the value when the surfaces of solute and solvent molecule are in contact. For the entropy of solvent molecule, the RFSDF of solvent has been used and the cut-off distance given in Eq. Ž5. has been used. For the entropy of the solute molecule, the RFSDF of the solute has been calculated separately and the cut-off distance
B.J. Yoon et al.r Chemical Physics Letters 301 (1999) 481–486
483
Table 1 The excess pressures, entropies and chemical potentials of solvent and solute of hard-sphere system at different solute–solvent ratios. Small numbers are the standard deviationa d
Solvent p 0ex VrRT
Solute ex
yS rR
ex
m rkT
p ex VrRT
yS exrR
Ž10y3 .
Ž10y3 .
mexrkT
mexrkT b
mexrkT c
0.059 0.104 0.169 0.255 0.368 0.502 0.678 0.886 1.146 1.334 1.753 2.163 2.606 3.056 3.589 4.230 4.845 5.563
0.054 0.095 0.153 0.232 0.334 0.464 0.624 0.817 1.047 1.317 1.630 1.989 2.397 2.858 3.378 3.963 4.610 5.315
0.054 0.095 0.153 0.232 0.334 0.464 0.623 0.815 1.044 1.311 1.621 1.977 2.381 2.837 3.348 3.916 4.546 5.240
0.232 0.453 0.776 1.210 1.816 2.639 3.540 4.720 6.151 7.913 9.748 11.839 14.524 16.579 20.894 24.309 27.610
0.235 0.443 0.764 1.224 1.850 2.669 3.711 5.004 6.586 8.457 10.621 13.122 15.4 – – – –
0.235 0.442 0.762 1.219 1.839 2.648 3.669 4.929 6.452 8.264 10.389 12.853 15.680 18.897 22.527 26.597 31.130
0.580 1.348 2.617 4.397
0.543 1.175 2.309 4.096
0.543 1.171 2.293 4.062
3
rs s 0.1 Ž10y4 .
Ž10y4 .
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4
0.356 3 Žy. d 0.266 2 Žy. 0.266 1Žy. 0.2652 Žy. 0.268 3 Žy. 0.268 2 Žy. 0.2674 Žy. 0.270 2 Žy. 0.271 2 Žy. 0.271 3 Žy. 0.274 3 Žy. 0.274 2 Žy. 0.276 3 Žy. 0.278 3 Žy. 0.280 1Žy. 0.282 3 Žy. 0.2844 Žy. 0.2852 Žy.
0.239 8 0.239 1 0.23510 0.241 4 0.2396 0.239 14 0.2338 0.242 9 0.2438 0.244 7 0.246 3 0.25114 0.25112 0.2559 0.256 9 0.261 9 0.26710 0.2695
rs 3 s 0.4 Ž10y3 .
Ž10y4 .
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2
1.134 3 Žy. 1.072 7 Žy. 1.071 4 Žy. 1.074 7 Žy. 1.0794 Žy. 1.083 4 Žy. 1.092 4 Žy. 1.100 7 Žy. 1.110 8 Žy. 1.120 5 Žy. 1.134 14 Žy. 1.1535 Žy. 1.171 6 Žy. 1.194 11Žy. 1.219 7 Žy. 1.2455 Žy. 1.290 6 Žy.
1.516 2 1.5236 1.499 3 1.534 2 1.520 4 1.512 3 1.541 4 1.561 4 1.578 5 1.560 2 1.613 2 1.6555 1.658 3 1.7254 1.811 5 1.842 10 1.920 7
rs 3 s 0.8 Ž10y3 .
Ž10y4 .
0.0 0.2 0.4 0.6
3.316 9 Žy. 3.299 21Žy. 3.313 26 Žy. 3.32314 Žy.
6.6079 6.604 9 6.650 9 6.658 10
0.595 0.505 0.501 0.506 0.507 0.507 0.500 0.512 0.514 0.515 0.520 0.525 0.527 0.533 0.535 0.543 0.550 0.554
2.650 2.596 2.569 2.608 2.599 2.595 2.633 2.671 2.687 2.680 2.747 2.808 2.855 2.920 3.030 3.087 3.210
9.923 9.903 9.963 9.981
0.028 4 0.049 2 0.0794 0.120 4 0.1755 0.236 4 0.321 9 0.42117 0.54948 0.686 6 0.834 17 1.052 28 1.271 56 1.47711 1.740 91 2.083133 2.40583 2.758 97
0.032 3 Žy. d 0.054 2 Žy. 0.090 3 Žy. 0.1354 Žy. 0.1936 Žy. 0.266 5 Žy. 0.358 9 Žy. 0.46517 Žy. 0.59717 Žy. 0.748 27 Žy. 0.914 8 Žy. 1.111 25 Žy. 1.33547 Žy. 1.579 39 Žy. 1.849 9 Žy. 2.14765 Žy. 2.460 51Žy. 2.84580 Žy.
Ž10y2 .
Ž10y3 .
0.132 1 0.262 1 0.4472 0.6952 1.051 5 1.5596 2.089 1 2.801 6 3.67710 4.873 22 6.08517 7.51456 9.569 28 11.008 44 14.74965 17.43444 20.34771
0.100 2 Žy. 0.191 6 Žy. 0.328 3 Žy. 0.5158 Žy. 0.7659 Žy. 1.080 20 Žy. 1.451 37 Žy. 1.91950 Žy. 2.47578 Žy. 3.040 39 Žy. 3.662 73 Žq. 4.325163 Žq. 4.95555 Žq. 5.571148 Žq. 6.145125 Žq. 6.876 183 Žq. 7.293 213 Žq.
Ž10y2 .
Ž10y3 .
0.3571 0.850 2 1.6644 2.830 6
0.223 3 Žy. 0.498 6 Žy. 0.95316 Žy. 1.56737 Žy.
B.J. Yoon et al.r Chemical Physics Letters 301 (1999) 481–486
484 Table 1 Žcontinued. d
Solvent
Solute
p 0ex VrRT
ex
yS rR
ex
m rkT
p ex VrRT
yS exrR
Ž10y2 .
Ž10y3 .
4.538 20 6.91743 9.42317 13.452 50 17.670 50 25.336 102 31.523125 44.655169 63.985183
2.366 80 Žy. 3.342 50 Žy. 4.336 62 Žy. 5.445167 Žq. 6.268 93 Žq. 7.320 190 Žq. 7.971180 Žq. 8.884 102 Žq. 9.864 200 Žq.
mexrkT
mexrkT b
mexrkT c
6.904 10.259 13.759 18.897 23.938 32.686 43.494 53.539 73.849
6.676 10.157 15.4 – – – – – –
6.630 10.153 14.781 20.669 27.971 36.838 47.425 59.885 74.370
3
rs s 0.8 Ž10y3 .
Ž10y4 .
0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4
3.3396 Žy. 3.383 24 Žy. 3.399 10 Žy. 3.46311Žy. 3.520 15 Žy. 3.60716 Žy. 3.676 26 Žy. 3.803 24 Žy. 3.956 45 Žy.
6.656 8 6.736 8 6.852 9 7.0955 7.108 15 7.448 20 7.692 7 8.14510 8.702 20
9.995 10.119 10.251 10.558 10.706 11.055 11.368 11.948 12.658
a
The standard deviations of bunches of 2 = 10 5 samplings for solvent and 1 = 10 5 samplings for solute are given by the small numbers, and the orders of the standard deviations are noted on the top of each column. b Data from Ref. w6x. c Calculated from BMCSL equation. d The sign of the second coefficient Ž B . of the function yln z Ž r . s Ar q Br 3.
given in Eq. Ž6. has been used. For the pressure of solvent, p 0ex , the following equation has been used: p 0ex V RT
s
ž
Ny1 2 N
/
3
prs 3 g 00 Ž 1 . q
1 p ex V N RT
sity rs 3 s 0.8 we sampled twice as many times a for low densities.
Ž 9.
where N is the number of molecules used in the simulation and g 00 Ž1. is the RDF of solvent–solvent at contact and the last term is the contribution of solvent–solute interaction in Eq. Ž8.. The major difference between this method and the insertion method is that the insertion method does not disturb the solvent structure whereas this method involves equilibration of one solute and the solvent molecule. Thus the chemical potential from the insertion method is an infinitely diluted one while the chemical potential from this method is not infinitely diluted but has a finite density of solute, 1rŽnumber of solvent used in the simulation.. The structure of solvent near the solute is disturbed in this method and thus the pressure of system increases with an increase in the size of solute Žsee Table 1.. We have used 108 molecules Ž107 solvent molecules and one solute molecule. in MC calculations and have taken an average of 2 = 10 6 samplings Žtrials of move. for the solvent molecule and 10y6 samplings for the solute molecule at densities of rs 3 s 0.1 and 0.4 after discarding 5 = 10 5 configurations starting from the initial face-centered-cubic structure. For the den-
3. Result and discussion In this work, the excess entropy and the excess pressure over an ideal gas have been calculated for both solute and solvent. The RDFs of solvent molecules around the solute molecule at rs 3 s 0.4 for different sizes of solute are shown in Fig. 1 and compared with the analytic expression of Boublik w13x 2 . When the size of solute increases the first peak in the RDF increases, i.e. the pressure in Eq. Ž8. increases. In Fig. 2, we also show the RFSDFs of the solute molecule in the cavity surrounded by solvent molecules at rs 3 s 0.4. From Fig. 2 we can see that the RFSDF of the solute drops off from that of an ideal gas Ži.e., unity for all r . as the sizes of solute become large, and thus the free volume decreases. In Fig. 2, the RFSDFs using Eq. Ž2. and the analytic two-cavity function w13x derived for infinite solute dilution are compared. The RFSDFs from the simplified two-cavity function by Ben-Amotz w14x
2
The coefficient Ž p r16. in Eq. Ž51. of w13x must be Ž p r12..
B.J. Yoon et al.r Chemical Physics Letters 301 (1999) 481–486
Fig. 1. The radial distribution functions of solvent molecules around a solute molecule at rs 3 s 0.4 for different sizes of solute. The solid lines are simulated ones and the dotted lines are calculated using the analytical expression of Ref. w13x. s is the diameter of solvent molecule.
have been also considered and overlap with Boublik’s form except for d s 0, but are not shown in Fig. 2. However, the equations have a singular value for d s 0. The pressure and entropy calculated using
Fig. 2. The radial free-space distribution functions of a solute surrounded by solvent molecules at rs 3 s 0.4 for different sizes of solute. For an ideal gas, the function is unity for all r values. The solid lines are simulated results and those from the analytic equation of Ref. w13x are represented by dotted lines. For calculating RFSDF at ds 0.0 from the analytic equation, ds 0.001 has been substituted in order to avoid a singular value. s is the diameter of solvent molecule.
485
these two functions ŽRDF and RFSDF, respectively. at rs 3 s 0.1, 0.4, and 0.8 are given in Table 1 and compared with the results by the insertion method and BMCSL equation. The standard deviations for the entropy are smaller than those for the pressure since the entropy is calculated by the integration of RFSDF as in Eq. Ž4. while the pressure is calculated from the single-point value of RDF as in Eqs. Ž8. and Ž9.. The excess chemical potentials of solutes in Table 1 are plotted in Fig. 3; however, the standard deviations are not plotted since they are smaller than the size of symbols used in the plot. The present method generally gives the slightly higher chemical potentials than those of the insertion method and BMCSL analytic-equation at rs 3 s 0.1. This general phenomenon at low densities has appeared in previous calculations w9,10x. This might be due to the fact that the definition of cut-off distance assuming a simple cubic lattice in Eqs. Ž5. and Ž6. is not correct Ža bit too large. at low densities. We inserted the solute molecule at the beginning and therefore the pressure of a system for large solute molecules is a bit larger than that of the system used in the insertion method w6x. At the densities rs 3 s 0.4 and 0.8, the chemical potentials for large solutes of the present method are less than those of the insertion method and BMCSL estimations. It might
Fig. 3. The excess chemical potentials are compared with those of the insertion method and BMCSL predictions. The results of this work are denoted by diamonds while circles and lines denote the results of the insertion method and BMCSL predictions, respectively.
486
B.J. Yoon et al.r Chemical Physics Letters 301 (1999) 481–486
be also due to the fact that the cut-off distance in Eq. Ž6. using the arithmatic mean is not properly defined for large solute molecule. Another possibility is that the large solute at high density involves something like a fluid-to-solid phase transition w11x of the disturbed solvent molecules near the solute, and thus the free volume of the solute becomes large. Such a phase transition is plausible since the solute and solvent molecules are initially mixed in this work and, for large solute, the function yln z Ž r . is concave, i.e., when yln z Ž r . is fitted to yln z Ž r . s Ar q Br 3 , the coefficient B is positive; while for small solute the function is convex Žsee Table 1.. In a previous work w11x, it was found that, for the fluid phase, yln z Ž r . was convex and, for the solid phase, it was concave. The deviations of analytical forms for large solutes in RDF and RFSDF Žin Figs. 1 and 2. from simulated values could also be for the same reason. If these calculations are reasonable, the entropies of insertion method are smaller Ži.e., larger negative value. and not valid for large solute at high densities. The insertion probabilities are evaluated to be smaller than the true values because the simulation generates equilibrated ensembles of solvent molecules but does not generate the unequilibrated, large vacancies that are suitable for a large solute to occupy, and they are not found effectively by checking the insertion probability of solute in equilibrated ensembles. One of the main advantages of this method is that we can calculate the pressure and the entropy separately in one simulation. Another importance is that we can estimate the chemical potentials of large solutes that cannot be calculated by the insertion method. The chemical potentials at high densities can also be calculated by the present method, and the RFSDFs are easily calculated during MC simulations by keeping the trajectories such as in Eq. Ž1. without requiring any other extra runs. When a small number of molecules are used in the simulation, the effect of the solute size on pres-
sure could be larger than in the case in which an infinite number of molecules is used in the simulation. Consequently, we are studying the effect of the number of molecules in the simulations. The reliability of the fluid-to-solid phase transition of solvent molecules near solute is going to be checked by calculating the RFSDF of solvent near a solute separately from that far from the solute. The results of these further studies will be presented in the near future.
Acknowledgements This study is supported by the academic research fund of Ministry of Education, Republic of Korea ŽBSRI-97-3414.. One of the authors ŽBJY. appreciates the partial support of Professor Ben-Amotz when he stayed at Purdue University.
References w1x R. Balescu, Equilibrium and Non-Equilibrium Statistical Mechanics, John Wiley, New York, 1975, p. 237. w2x B. Widom, J. Chem. Phys. 39 Ž1963. 2808. w3x B. Widom, J. Phys. Chem. 86 Ž1982. 869. w4x M.P. Allen, D.J. Tildesley, Computer Simulation of Liquid, Oxford, New York, 1991. w5x L.E.S. de Souza, A. Stamatopoulou, D. Ben-Amotz, J. Chem. Phys. 100 Ž1994. 1456. w6x A. Stamatopoulou, L.E.S. de Souza, D. Ben-Amotz, J. Talbot, J. Chem. Phys. 102 Ž1995. 2109. w7x T. Boublik, J. Chem. Phys. 53 Ž1970. 471. w8x G.A. Mansoori, N.F. Carnahan, K.E. Starling, T.W. Leland Jr., J. Chem. Phys. 54 Ž1971. 1523. w9x B.J. Yoon, H.A. Scheraga, J. Mol. Struct. ŽTHEOCHEM. 199 Ž1989. 33. w10x B.J. Yoon, S.D. Hong, M.S. Jhon, H.A. Scheraga, Chem. Phys. Lett. 181 Ž1991. 73. w11x B.J. Yoon, Chem. Phys. Lett. 234 Ž1995. 35. w12x A. Rotenberg, J. Chem. Phys. 43 Ž1965. 4377. w13x T. Boublik, Mol. Phys. 59 Ž1986. 775. w14x D. Ben-Amotz, J. Phys. Chem. 97 Ž1993. 2314.