- Email: [email protected]
Activity coefficient of the dipolar Yukawa and Stockmayer fluids from the inverse grand-canonical Monte Carlo simulation
Fluid Phase Equilibria 317 (2012) 29–35
Contents lists available at SciVerse ScienceDirect
Fluid Phase Equilibria journal homepage: www.elsevier.com/locate/fluid
Activity coefficient of the dipolar Yukawa and Stockmayer fluids from the inverse grand-canonical Monte Carlo simulation S. Lamperski ∗ , R. Górniak Department of Physical Chemistry, Faculty of Chemistry, A. Mickiewicz University, Grunwaldzka 6, 60-780 Pozna´ n, Poland
a r t i c l e
i n f o
Article history: Received 5 October 2011 Received in revised form 10 December 2011 Accepted 16 December 2011 Available online 24 December 2011 Keywords: Dipolar Yukawa fluid Stockmayer fluid Inverse grand-canonical Monte Carlo simulation (IGCMC) Activity coefficient Mean spherical approximation (MSA) Perturbation theory (PT)
a b s t r a c t The inverse grand-canonical Monte Carlo (IGCMC) simulation is applied to calculate the activity coefficient of the dipolar Yukawa and Stockmayer fluids. The dependence of the activity coefficient and its dipolar contribution on the number density of a fluid at the reduced temperature T* = 3.0 and reduced electric dipole moment (* )2 = 1.0, 3.0, and 5.0 is presented. The results are compared with predictions of: (i) the perturbation theory (PT), (ii) the free energy and virial routes of the mean spherical approximation (MSA), and (iii) Kronome–Liszi–Szalai (KLS) approach (J. Chem. Soc., Faraday Trans. 93 (1997), 3053–3059). Comparison shows that at (* )2 = 5.0 the PT results are in a very good agreement with the IGCMC data. The KLS approach yields the results convergent with the simulation one for low number densities of the fluid while the MSA free energy and virial routes underestimate and overestimate the dipolar interactions, respectively. © 2012 Elsevier B.V. All rights reserved.
1. Introduction The activity is one of the fundamental thermodynamic properties. Its component, the activity coefficient, gives insight into intermolecular interactions. Also, the activity coefficient is needed in the grand-canonical Monte Carlo (GCMC) simulations of inhomogeneous systems to hold constant the bulk concentration [1]. In the previous paper [2], we have applied the inverse grand-canonical Monte Carlo (IGCMC) simulation, developed recently by Lamperski [3], to calculate the activity coefficient of a dipolar hard-sphere (DHS) fluid. Good agreement between the IGCMC and Widom method [4] results was observed. Comparison with the mean spherical approximation (MSA) showed that this theory underestimated the dipolar interactions while the agreement with the perturbation theory (PT) was good. The dipolar Yukawa (DY) and Stockmayer (STM) fluids are more reasonable models of polar fluids than DHS as in these models the non-electrostatic attractions are included. The models were widely used in the calculation of the dielectric constant [5–7], the internal energy and pressure [5,8,9], the heat capacity [10,11], phase equilibrium in an external field [12] and system with a high
∗ Corresponding author. Tel.: +48 61 82 91 454; fax: +48 61 82 91 505. E-mail address: [email protected] (S. Lamperski). 0378-3812/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2011.12.017
dipole moment [13,14]. Computer simulations of the free energy and activity coefficient have been rarely performed. The umbrella sampling [15] and the self-consistent histogram [16] methods were applied to calculate the free energy. The chemical potential of STM fluid was evaluated from the MC simulation in the Gibbs ensemble [17–19]. Here we apply the IGCMC method for the calculation of activity coefficients of DY and STM fluids. The advantage of IGCMC over the classical GCMC simulation [20] is that the former allows evaluation of the activity coefficient for a specified concentration (number density) of a system. A similar possibility is available by applying the algorithm propose by Malasics et al. [21]. The new technique has been successfully applied to the primitive model electrolyte [3], to the solvent primitive model electrolyte [22], and to DHS fluid [2]. The IGCMC mean activity coefficients of the model electrolyte were applied in the GCMC simulations of different inhomogeneous systems [23–25]. The new technique allows us to calculate not only the mean activity coefficient of an electrolyte, but also individual coefficients. For polar fluids, most of the thermodynamic functions had been studied earlier except for the activity coefficients. So, our interest is to fill this gap. Recently, we have considered DHS fluid [2]. In this paper, we analyze the activity coefficients of DY and STM potentials at different electric dipole moments and densities of a fluid, calculated by the IGCMC technique. Also, we present a comparison with the MSA and PT results.
30
S. Lamperski, R. Górniak / Fluid Phase Equilibria 317 (2012) 29–35
the point electric dipole moment is located at the centre of a hard sphere. The potential of the dipole–dipole interactions is
Nomenclature Add Add , Addd 2 3 Aex A* B b2 C d dT d2 G Idd , Iddd kB k1 N N0 n p r T u V y zdd zex ˛ ε0 ε εY C HS ex *
dipole contribution to the Helmholtz energy second and third perturbation term excess free energy reduced free energy term in the Ewald sum MSA parameter term in Ewald sum hard-sphere diameter temperature dependent hard-sphere diameter MSA parameter reciprocal lattice vector integrals in PT Boltzmann constant coefficient in IGCMC technique number of molecules in the simulation box number of molecules in the simulation box corresponding to the specified number density integer vector in Ewald sum pressure distance between two molecules absolute temperature potential energy volume dipole strength compression factor excess compression factor Gaussian distribution parameter vacuum permittivity the depth of well in the Lennard–Jones potential the depth of well in the Yukawa potential activity coefficient Coulombic contribution to the activity coefficient hard-sphere contribution to the activity coefficient packing fraction the range of attractions in the Yukawa potential electric dipole moment, chemical potential excess chemical potential number density reduced number density coupling parameter the separation at which the Lennard–Jones potential uLJ = 0
2. The model
uD (rij ) =
[i · j rij2 − 3(i · rij )(j · rij )] 4 ε0 rij5
(2)
where ε0 is the vacuum permittivity. Thus the DY potential, uDY , is the sum of uY and uD potentials. In the STM model, hard spheres are replaced by soft spheres. The point electric dipole moment is located at the centre of the soft sphere. The soft and dispersion interactions are given by the LJ potential
12 rij
uLJ (rij ) = 4ε
6
−
rij
(3)
where ε and are the LJ potential parameters (the depth of the potential well and the separation at which uLJ = 0, respectively). The STM potential, uSTM , consists of uLJ and uD . The LJ soft-core potential makes the STM model more realistic than DY. We use the following reduced properties for the DY potential: (∗ )2 =
2 , 4 ε0 εY d3
T∗ =
kB T , εY
∗ = d3
(4)
∗ = 3 .
(5)
and for the STM model: (∗ )2 =
2 , 4 ε0 ε 3
T∗ =
kB T , ε
where kB is the Boltzmann constant, T is temperature and is the number density of molecules. 3. Simulation methods and theories 3.1. The inverse grand-canonical Monte Carlo method The details of the GCMC simulation of polar fluids have been described for example in our previous paper [2]. The disadvantage of the GCMC technique is that it does not allow evaluation of the activity coefficient for a specified concentration (number density). This disadvantage is overcome in the IGCMC method [3]. In IGCMC, we divide the whole simulation into m short GCMC simulations in each of which the number n of configurations is generated. The simulation starts with a reasonable initial value of ln 0 . After each simulation the average number of molecules is calculated. If this number is greater than the number of molecules N0 in the simulation box corresponding to the specified number density , then a new value of ln i is obtained by subtracting from the previous ln i − 1 a small value k (e.g. 0.02). Otherwise ln i − 1 is increased by k. After the equilibration is reached, the average ln i over the m short simulations is calculated 1 ln i m m
We consider the dipolar Yukawa (DY) and Stockmayer (STM) fluids. The Yukawa potential, uY , is the sum of the hard-sphere and attractive (dispersion) interactions. For two molecules at positions ri and rj , uY is given by
uY (rij ) =
∞, εY d exp[−(rij − d)], − rij
rij < d rij ≥ d
(1)
where rij is the distance between the molecules i and j, d is the diameter of a hard sphere, and εY and are the depth and range of the attraction interactions, respectively. The properties of Yukawa fluid are similar to those of the Lennard–Jones (LJ) fluid when = 1.8/d [26] and this value for we will use in this work. In the DY model,
ln =
(6)
i=1
3.2. Thermodynamics The activity coefficient is related to the excess chemical potential ex by the formula ex = kB T ln
(7)
The excess chemical potential can be obtained from the excess free energy, Aex , using the thermodynamic relation
ex =
∂Aex ∂N
(8) V,T
S. Lamperski, R. Górniak / Fluid Phase Equilibria 317 (2012) 29–35
Let us introduce the reduced free energy, A* , per molecule A∗ =
and
A NkB T
(9)
Using Eqs. (7)–(9) we obtain
1 ln = kB T
∂Aex ∂N
=
A∗ex
+N
V,T
∂A∗ex ∂N
∂A∗ex ∂N
= V,T
∂A∗ex ∂
= ∗ T
∂A∗ex ∂∗
(10) V,T
= T
∂A∗ex ∂
(11) T
where is the packing fraction =
∗ 6
(12)
The activity coefficient can be also calculated from the excess contribution to the compression factor, zex (the excess against the perfect gas). The thermodynamic relation, which tells that pressure is a negative value of the derivative of free energy against volume, can be written in terms of the derivative of the number density
−p =
∂A ∂V
∂A ∂
= N,T
T
∂ ∂V
2 =− N
N
∂A ∂
(13) T
∂ ∂V
=− N
2 N
(14)
From the equality of the first and the last term of Eq. (13), and using Eq. (9) we obtain the expression for A∗ex
A∗ex
zex d =
= 0
∗
zex ∗ d = ∗
0
0
zex d
∂A∗ex ∂
ln =
(16)
dT =
zex d + zex
(17)
APT,dd DHS NkB T
=
Add /NkB T 2 1 − Addd /Add 3 2
(18)
Here Add 2 NkB T Addd 3 NkB T
1 − exp
−uLJ (rij )
kB T
drij
(23)
Cotterman et al. [31] approximated this equation by the expression which for the spherical LJ potential takes the form dT =
1 + 0.29770T ∗ 1 + 0.33163T ∗ + 0.0010477T ∗2
rmin
0 rmin
rij2 yHS (rij ) exp
rij2 yHS (rij ) exp
−
(24)
−uHS (rij )
kB T
drij
−(uLJ (rij ) + ε)
drij = 0
kB T
(25)
0
uHS (rij ) =
∞, 0,
rij < d rij ≥ d
(26)
and yHS is the Percus–Yevick y function [33–35] yHS (rij ) = 1 + 62
r ij
d
+ 12 1
3 rij d
(27)
with 1 =
In the PT theory of DY fluid, the Yukawa potential is considered as the reference potential and the dipole–dipole interactions as a perturbation. The free energy of the dipole–dipole interactions, Add , DY can be approximated by the corresponding free energy Add of the DHS DHS fluid. Henderson et al. [27] obtained the analytical solution of PT with the Padé approximation given by Stell et al. [28]
NkB T
0
3.3. The perturbation theory
≈
T
0
APT,dd DY
Using the above two results, the expression (10) for the activity coefficient written in terms of zex takes the form
where rmin = 21/6 , uHS is the hard-sphere potential
= zex
(22)
When APT,dd is known, the activity coefficient can be calculated DY from Eq. (10) with the numerically evaluated derivative. In the STM fluid the LJ potential is treated as the reference system, and the dipole–dipole interaction as a perturbation [28,29]. According to Barker and Henderson (BH) [30] the free energy of a soft-sphere fluid can be approximated by that of a hard-sphere fluid with the temperature dependent hard-sphere diameter, dT . For the LJ potential the Authors suggested calculating this diameter from the equation
(15)
and for its derivative
5(1 + 1.12754d3 + 0.561922 d6 ) 24(1 − 0.05495d3 + 0.133322 d6 )
(21)
and is valid for 0 < T* < 15. Another approach for dT was proposed by Weeks, Chandler, and Anderson (WCA) [32]. It requires solution of the equation
as at N hold constant
1 + 0.18158d3 − 0.114672 d6 3(1 − 0.49303d3 + 0.062932 d6 )
I dd =
I ddd =
The last term of Eq. (10) can be rewritten in terms of , * , or N
31
2 =
(1 + 2)2
(28)
(1 − )4 −(1 + /2)
2
(1 − )4
(29)
The diameter dT calculated from Eq. (25) is both temperature and density dependent. Now the dipolar free energy of STM, APT,dd , can be obtained from STM Eq. (18) where the hard-sphere d in Eqs. (19)–(22) is replaced with dT calculated from Eq. (23) or (25). The temperature dependence of diameter means that * and also depend on temperature, but does not. 3.4. The MSA theory
=−
=
2 4 3d3 (4 ε0 kB T )
2
4 2 2 6 27d3 (4 ε0 kB T )
3
I
dd
(19)
The MSA free energy of DY fluid, AMSA , can be separated into the DY
hard-sphere Yukawa AMSA and dipolar AMSA,dd contributions [36] Y DY I
ddd
(20)
MSA AMSA + AMSA,dd DY = AY DY
(30)
32
S. Lamperski, R. Górniak / Fluid Phase Equilibria 317 (2012) 29–35
The dipolar contribution AMSA,dd is approximated by the dipoDY lar hard-sphere free energy,
MSA, dd ADHS .
7
The analytical solution for
6
was given by Wertheim [37] (the free energy route) AMSA,dd DHS AMSA,dd DY NkB T
AMSA,dd DHS
≈
NkB T
3 = − I(y)
5 (31)
y=
ln
where y is the dipole strength 2
(32)
9ε0 kB T
8 2 3
(1 + )
2
(1 − 2)
+
4
(2 − )
2
8(1 + )
(33)
4
(1 + 4)
2
(1 − 2)
4
−
(1 − 2) (1 + )
2
4
= 3y
(34)
3 3 I(y) − y
(35)
, This expression is obtained by differentiating the free energy, AMSA DHS so the ln dd results obtained from the free energy (Eq. (31)) and from the above formula are exactly the same. Another expression for the compression factor is from the virial (pressure) route [38]. It gives dd zDY, v
≈
dd zDHS, v
3 = − y
(36)
As we see, different routes give different results for the compression factor. This is the example of the so-called thermodynamic inconsistency which stems from the fact that the theory is approximate. An interesting MSA approach based on the perturbation theory was given by Kronome, Liszi and Szalai (KLS) [39]. By adding Eqs. (35) and (36) with weights of 2/3 and 1/3, respectively, like Carnahan and Starling [40] did, the Authors proposed the following expression for the dipolar contribution to the compression factor: dd dd zDY,KLS ≈ zDHS,KLS =
2 3 I(y) − y
(37)
In this work, we use the above three approaches to calculate the activity coefficient from the compression factor using Eq. (17). As there is no analytical solution of MSA for the STM fluid [5] dd from Eqs. (31) and (37), where one can evaluate AMSA,dd and zSTM STM is replaced with the temperature dependent T as described in Section 3.3. 3.5. Long-range interactions The long-range dipolar interactions in IGCMC were handled with the Ewald sum (ES) [41] methods. In the ES method, the total potential energy of dipole–dipole interactions consists of a real space sum (the first term of Eq. (38)), reciprocal space sum (the second term) reduced by the self-interaction term (the third term) 1 1 uij = 3 2 rij + n i,j=1 n N
/ 0 i,j=1G =
0.0
0.2
0.4
0.6
0.8
Fig. 1. Dependence of the IGCMC activity coefficient ln on * for the dipolar Yukawa fluid at T* = 3.0 and for the following (* )2 : 1.0 (circles), 3.0 (squares), and 5.0 (triangles). The line through the IGCMC points is given as a visual aid.
2˛r B(r) = erfc(˛r) + √ exp(−˛2 r 2 )
(39)
2˛r C(r) = erfc(˛r) + √ exp(−˛2 r 2 )(1 + 23 ˛2 r 2 )
(40)
where n = (nx , ny , nz ), parameter ˛ regulates the relative convergence of the sum. Reciprocal lattice vectors are given by G = ((2 /L)nx , (2 /L)ny , (2 /L)nz ). Eqs. (38)–(40) are taken from [42]. 4. Results and discussion The IGCMC simulations for DY and STM fluids were carried out for a reduced temperature T* = 3.0, reduced electric dipole moment (* )2 = 1.0, 3.0, and 5.0, and reduced density * ranging from 0.1 to 0.8. The temperature T* is higher than the critical temperature Tc∗ for STM fluid which for the values of (* )2 considered is 1.41 [17], 1.82 [18], and 2.31 [19], respectively. For T = 298.15 K and d = = 400 pm, ε = εY = 826.33 J/mol, while the dipole moment in D corresponding to the considered values of (* )2 is 0.9371, 1.6231, and 2.0955. Earlier, we have found [2] that the activity coefficient results for dipolar hard-sphere fluid are weakly dependent on the number of molecules, N0 . Thus in the present simulations we used N0 = 256. The ES method was applied to describe the long-range dipole–dipole interactions. The value of parameter ˛ was 5.1. The sum in the reciprocal space extended over all values of G such that |n|2 ≤ 27. The number of simulation steps applied to equilibrate the system was 6,000,000 while in the main part of the simulation this number amounted to 10,000,000. To reach faster a value of ln close to the expected one, we applied k1 = 0.05 in the first part of the simulation, while in the second part k1 = 0.02. The standard deviation of ln depended on * and varied from 0.02 to 0.03. The IGCMC results for the DY fluid are shown in Fig. 1 and listed in Table 1. With increasing value of (* )2 the activity coefficient decreases. The curve for (* )2 = 1.0 has the exponential-like shape similar to that for a the hard-sphere fluid [22]. For higher (* )2 a shallow minimum with negative values is observed at low * .
i · (rij + n) · j · (rij + n) i · j B ˛ rij + n − 3 C ˛ rij + n rij + n2
2 −G2 /(4˛2 ) (i · G)(j · G) 2˛3 + e exp(iG · r ) − √ ij V G2 3
N
0
According to Eq. (17) the activity coefficient ln dd can be also calculated from the dipolar contribution to the compression factor, zdd . In the free energy route [38] dd dd zDY,A ≈ zDHS,A =
2
-1
The parameter is obtained by numerical solution of the following equation:
3
1
and the function I(y) is given by I(y) =
4
i=1....N
2i
(38)
S. Lamperski, R. Górniak / Fluid Phase Equilibria 317 (2012) 29–35 Table 1 The IGCMC activity coefficient ln for the dipolar Yukawa fluid at T* = 3.0. *
4
ln
3
(* )2 = 0.0
(* )2 = 1.0
(* )2 = 3.0
(* )2 = 5.0
0.06768 0.2124 0.4643 0.8700 1.5016 2.4607 3.9215 6.0987
0.05006 0.1748 0.4058 0.7843 1.3919 2.30869 3.7566 5.9270
−0.09116 −0.09450 0.00972 0.2551 0.7157 1.4861 2.7280 4.7140
−0.3718 −0.5751 −0.6373 −0.5452 −0.2508 0.3386 1.4532 3.3095
2
ln
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
33
1 0 -1 -2
0.0
0.0
0.2
0.4
0.6
0.8
-0.5 Fig. 4. Dependence of the IGCMC activity coefficient ln on * for the Stockmayer fluid at T* = 3.0 and for the following (* )2 : 1.0 (circles), 3.0 (squares), and 5.0 (triangles). The line through the IGCMC points is given as a visual aid.
ln
dd
-1.0 -1.5 -2.0 -2.5 -3.0 0.0
0.2
0.4
0.6
0.8
Fig. 2. Dependence of the dipolar contribution to the activity coefficient ln dd on * for the dipolar Yukawa fluid. Physical parameters and symbols as in Fig. 2.
The activity coefficient of DY can be separated, like the free energy, into the Yukawa ln Y and dipolar ln dd contributions. The ln dd results are shown in Fig. 2. These results take negative values and a nearly linear dependence of ln dd on * is observed. In Fig. 3, we compare the IGCMC dipolar contribution (triangles- the solid line is given as a visual aid) with the corresponding PT (a dashed line) and MSA (short dashed, dotted-dashed, and dotted-dotteddashed lines) data for (* )2 = 5.0. The ln dd results calculated by the MSA virial route (a dotted-dashed line) are evidently overestimated. Better agreement with the computer simulation results
0
ln
dd
-1
-2
is observed for ln dd obtained from the MSA free energy route (a short dashed line) although this time the theoretical results are underestimated. The curve following from the KLS approximation (a dotted–dotted-dashed line) runs close to the IGCMC line at * < 0.5 while at high * some divergence is observed which intensifies with increasing * . The PT results are in excellent agreement with the IGCMC data in the whole range of * . Earlier, for a dipolar hard-sphere fluid, we also observed [2] a great similarity of PT and IGCMC results. In Fig. 4 and in Table 2, the IGCMC results for the STM fluid are presented. The curves shown in Fig. 4 are similar to those in Fig. 1 although the value of the activity coefficient in Fig. 4 is lower and the minimum of STM curves is deeper. It is due to a stronger dispersion attraction in LJ than in Y fluid (compare the second columns of Tables 1 and 2). The ln dd curves shown in Fig. 5 exhibit similar behaviour like DY curves but they are slightly more bent. In the next two graphs, we compare the ln dd data for the STM fluid obtained by the IGCMC simulation with the PT and MSA results for (* )2 = 5.0. Fig. 6 shows the theoretical results for dT calculated from the BH formula (Eq. (23)) which for T* = 3.0 yields dT / = 0.94505. Again the MSA virial and free energy route curves run far from the simulation data. The KLS curve is in much better agreement with the IGCMC results but the curvatures of these curves are different. The PT curve is very close to the IGCMC one in the whole range of * . The value of dT / calculated from the WCA theory (Eq. (25)) depends on the density and for T* = 3.0 it varies from 0.97312 for * = 0.10 to 0.96900 for * = 0.80. The MSA and PT activity coefficients calculated with theses values of dT are shown in Fig. 7. As seen, the WCA theory does not improve the convergence of the MSA virial and free energy route results with the simulation data. However, the KLS approach at * < 0.4 and the PT theory at * > 0.4 give nearly perfect agreement with the IGCMC results. Table 2 The IGCMC activity coefficient ln for the Stockmayer fluid at T* = 3.0.
-3
*
-4 0.0
0.2
0.4
0.6
0.8
Fig. 3. Dependence of the dipolar contribution to the activity coefficient ln dd on * for the dipolar Yukawa fluid at T* = 3.0, and (* )2 : = 5.0. The symbols represent the IGCMC data while the lines the different theories results (PT – dashed line, MSA (free energy route) – short dashed line, MSA (virial route) – dotted–dashed line, and KLS – dotted–dotted–dashed line).
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
ln (* )2 = 0.0
(* )2 = 1.0
(* )2 = 3.0
(* )2 = 5.0
−0.02312 0.01136 0.12181 0.3381 0.7059 1.2886 2.1796 3.4833
−0.04667 −0.03467 0.05276 0.2431 0.5855 1.1394 1.9774 3.2424
−0.2077 −0.3368 −0.3818 −0.3239 −0.1310 0.2616 0.9345 2.0245
−0.5497 −0.8999 −1.1150 −1.2217 −1.1897 −0.9684 −0.4798 0.4017
34
S. Lamperski, R. Górniak / Fluid Phase Equilibria 317 (2012) 29–35
0.0
0
-0.5
dd
-1
-1.5
ln
ln
dd
-1.0
-2.0
-2
-2.5 -3.0
-3
-3.5 0.0
0.2
0.4
0.6
0.8
Fig. 5. Dependence of the dipolar contribution to the activity coefficient ln dd on * for the Stockmayer fluid. Physical parameters and symbols as in Fig. 4.
0.0
0.2
0.4
0.6
0.8
Fig. 8. Dependence of the IGCMC dipolar contribution to the activity coefficient ln dd on * at T* = 3.0, and (* )2 = 5.0 for the dipolar hard-sphere (circles, solid line), dipolar Yukawa (squares, dashed line) and Stockmayer (triangles, short dashed line) fluids.
0
ln
dd
-1
-2
-3
-4 0.0
0.2
0.4
0.6
0.8
Fig. 6. Dependence of the dipolar contribution to the activity coefficient ln dd on * for the Stockmayer fluid at T* = 3.0, and (* )2 : = 5.0. The symbols represent the IGCMC data while the lines the different theories results with dT given by Eq. (23) (PT – dashed line, MSA (free energy route) – short dashed line, MSA (virial route) – dotted–dashed line, and KLS – dotted–dotted–dashed line).
Although the DHS fluid is not a subject of this paper, it is interesting to compare the properties of the following dipolar fluids: DHS, DY, and STM. For theses models, the IGCMC results of ln dd at *2 = 5.0 are collected in Fig. 8. The DHS and DY results are nearly identical. It is not surprising as in both cases the repulsion interactions are described by the same hard-sphere potential. The STM curve runs below the DHS and DY lines. As follows from the BH formula or WCA theory, the average closest approach of two STM (or LJ) molecules is smaller than . Thus for = d, as we assumed in this paper, the point dipole moments of STM molecules can approach closer to each other than the dipoles of DHS or DY fluids and therefore the dipole–dipole interactions are stronger. This lowers ln dd of STM in comparison to DHS or DY. A separate problem is whether the dispersion terms in DY and STM potentials, which give rise to the intermolecular attractions, influence the ln dd results. An indepth analysis of the DHS and DY curves implies the positive answer to this question. At higher densities a subtle but permanent separation of the curves is seen with the DY curve running below the DHS one. As dispersion is a short range interaction, it should be manifested at a high density, as shown in Fig. 8. However, the separation is small and can be comprised in the error bar, so the problem needs further, much more precise investigation.
0
5. Conclusions
ln
dd
-1
-2
-3
-4 0.0
0.2
0.4
0.6
0.8
Fig. 7. Dependence of the dipolar contribution to the activity coefficient ln dd on * for the Stockmayer fluid at T* = 3.0, and (* )2 : = 5.0. The symbols represent the IGCMC data while the lines the different theories results with dT given by Eq. (25) (PT – dashed line, MSA (free energy route) – short dashed line, MSA (virial route) – dotted–dashed line, and KLS – dotted–dotted–dashed line).
In this paper, we have shown that the IGCMC method can be successfully applied to calculate the activity coefficient of the DY and STM fluids. Subtraction of the activity coefficients of the Y and LJ fluids from these results gave the dipolar contributions to the overall activity coefficient. Comparison of the dipolar contribution yielded by the IGCMC simulation with that obtained by the PT, MSA and KLS theories showed that for the DY fluid the PT results are in a very good agreement with the simulation data in the whole range of * . At low * the KLS approach also gives the correct results. For the STM fluid with dT calculated according to the WCA theory, the KLS approach at * < 0.4 and the PT theory at * > 0.4 give nearly a perfect agreement with the IGCMC results. Finally, we should note that our conclusion that the MSA virial route underestimates the magnitude of the activity coefficient while PT is much more accurate is consistent with the results for other thermodynamic properties. A great similarity of the DHS and DY results shown in Fig. 8 supports the assumption that the free energy of the dipole–dipole interactions of the DY fluid can be approximated by the
S. Lamperski, R. Górniak / Fluid Phase Equilibria 317 (2012) 29–35
corresponding free energy of DHS fluid (Eqs. (18) and (31)). Explanation of a subtle difference observed at * > 0.4 requires low temperature investigation which will better separate the dipolar and dispersion attraction from the temperature independent hardsphere interactions. It could throw light on the role of dispersion attraction in the dipole–dipole interaction of DY fluid. Acknowledgments
[14] [15] [16] [17] [18] [19] [20] [21] [22] [23]
Financial support from Adam Mickiewicz University, Faculty of Chemistry, is appreciated.
[24] [25] [26]
References
[27]
[1] D. Frenkel, B. Smit, Understanding Molecular Simulations: From Algorithms to Applications, Academic Press, San Diego, 1996, p. 121. [2] S. Lamperski, R. Górniak, Fluid Phase Equilib. 295 (2010) 255–263. [3] S. Lamperski, Mol. Simul. 33 (2007) 1193–1198. [4] B. Widom, J. Chem. Phys. 39 (1963) 2808–2812. [5] I. Szalai, D. Henderson, D. Boda, K.-Y. Chan, J. Chem. Phys. 111 (1999) 337–344. [6] I. Szalai, K.-Y. Chan, Y.W. Tang, Mol. Phys. 101 (2003) 1819–1828. [7] L.E. Johnson, R. Barnes, T.W. Draxler, B.E. Eichinger, B.H. Robinson, J. Phys. Chem. B 114 (2010) 8431–8440. [8] C. Kriebel, Winkelmann, J. Mol. Phys. 88 (1996) 559–578. [9] Y.V. Kalyuzhnyi, I.A. Protsykevytch, P.T. Cummings, Europhys. Lett. 80 (2007) 56002. [10] Z. Máté, I. Szalai, Fluid Phase Equilib. 54 (2010) 54–60. [11] Z. Máté, I. Szalai, D. Boda, D. Henderson, Mol. Phys. 109 (2011) 203–208. [12] A. Gábor, I. Szalai, Mol. Phys. 106 (2008) 1094–1106. [13] Yu.V. Kalyuzhnyi, Y.V. Protsykevitch, Condens. Mater. Phys. 10 (2007) 553–562.
[28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42]
35
Yu.V. Kalyuzhnyi, G. Stell, Mol. Phys. 78 (1993) 1247–1258. G.M. Torrie, J.P. Valleau, J. Comp. Phys. 23 (1977) 187–199. C.H. Bennet, J. Comp. Phys. 22 (1976) 245–268. B. Smit, C.P. Williams, E.M. Hendriks, S.W. de Leeuw, Mol. Phys. 68 (1989) 765–769. M.E. van Leeuwen, B. Smit, E.M. Hendriks, Mol. Phys. 78 (1993) 271–283. M.E. van Leeuwen, Mol. Phys. 82 (1994) 383–392. J.P. Valleau, L.K. Cohen, J. Chem. Phys. 72 (1980) 5935–5942. A. Malasics, D. Gillespie, D. Boda, J. Chem. Phys. 128 (2008) 124102–124108. S. Lamparski, M. Płuciennik, Mol. Simul. 36 (2010) 111–117. S. Lamperski, C.W. Outhwaite, L.B. Bhuiyan, J. Chem. Phys. B 113 (2009) 8925–8929. S. Lamperski, C.W. Outhwaite, J. Colloid Interface Sci. 328 (2008) 458–462. S. Lamparski, J. Kłos, J. Chem. Phys. 129 (2008) 164503–164510. D. Henderson, E. Waisman, J.L. Lebowitz, L. Blum, Mol. Phys. 35 (1978) 241–255. D. Henderson, L. Blum, A. Tani, Advances in Chemistry Series 300, American Chemical Society, Washington, DC, 1986, pp. 281–296. G. Stell, J.C. Rasaiah, H. Narang, Mol. Phys. 27 (1974) 559–564. G. Stell, J.C. Rasaiah, H. Narang, Mol. Phys. 27 (1974) 1393–1414. J.A. Barker, D. Henderson, J. Chem. Phys. 47 (1967) 2856–2862. R.L. Cotterman, B.J. Schwarz, J.M. Prausnitz, AIChE J. 32 (1986) 1787–1798. J.D. Weeks, D. Chandler, H.C. Anderson, J. Chem. Phys. 54 (1971) 5237–5247. M.S. Wertheim, Phys. Rev. Lett. 10 (1963) 321–323. M.S. Wertheim, J. Math. Phys. 5 (1964) 643–652. E. Thiele, J. Chem. Phys. 39 (1963) 474–480. D. Henderson, D. Boda, I. Szalai, K.Y. Chan, J. Chem. Phys. 110 (1999) 7348–7358. S.M. Wertheim, J. Chem. Phys. 55 (1971) 4291–4298. C.G. Gray, K.E. Gibbson, Theory of molecular fluids Fundamentals, vol. 1, Oxford University Press, New York, 1984, p372. G. Kronome, J. Liszi, I. Szalai, J. Chem. Soc. Faraday Trans. 93 (1997) 3053–3059. F. Carnahan, E. Starling, J. Chem. Phys. 53 (1969) 635–637. P. Ewald, Ann. Phys. 64 (1921) 253–287. W. Smith, CCP5 Newsletter 46 (1998) 18–23.
Contents lists available at SciVerse ScienceDirect
Fluid Phase Equilibria journal homepage: www.elsevier.com/locate/fluid
Activity coefficient of the dipolar Yukawa and Stockmayer fluids from the inverse grand-canonical Monte Carlo simulation S. Lamperski ∗ , R. Górniak Department of Physical Chemistry, Faculty of Chemistry, A. Mickiewicz University, Grunwaldzka 6, 60-780 Pozna´ n, Poland
a r t i c l e
i n f o
Article history: Received 5 October 2011 Received in revised form 10 December 2011 Accepted 16 December 2011 Available online 24 December 2011 Keywords: Dipolar Yukawa fluid Stockmayer fluid Inverse grand-canonical Monte Carlo simulation (IGCMC) Activity coefficient Mean spherical approximation (MSA) Perturbation theory (PT)
a b s t r a c t The inverse grand-canonical Monte Carlo (IGCMC) simulation is applied to calculate the activity coefficient of the dipolar Yukawa and Stockmayer fluids. The dependence of the activity coefficient and its dipolar contribution on the number density of a fluid at the reduced temperature T* = 3.0 and reduced electric dipole moment (* )2 = 1.0, 3.0, and 5.0 is presented. The results are compared with predictions of: (i) the perturbation theory (PT), (ii) the free energy and virial routes of the mean spherical approximation (MSA), and (iii) Kronome–Liszi–Szalai (KLS) approach (J. Chem. Soc., Faraday Trans. 93 (1997), 3053–3059). Comparison shows that at (* )2 = 5.0 the PT results are in a very good agreement with the IGCMC data. The KLS approach yields the results convergent with the simulation one for low number densities of the fluid while the MSA free energy and virial routes underestimate and overestimate the dipolar interactions, respectively. © 2012 Elsevier B.V. All rights reserved.
1. Introduction The activity is one of the fundamental thermodynamic properties. Its component, the activity coefficient, gives insight into intermolecular interactions. Also, the activity coefficient is needed in the grand-canonical Monte Carlo (GCMC) simulations of inhomogeneous systems to hold constant the bulk concentration [1]. In the previous paper [2], we have applied the inverse grand-canonical Monte Carlo (IGCMC) simulation, developed recently by Lamperski [3], to calculate the activity coefficient of a dipolar hard-sphere (DHS) fluid. Good agreement between the IGCMC and Widom method [4] results was observed. Comparison with the mean spherical approximation (MSA) showed that this theory underestimated the dipolar interactions while the agreement with the perturbation theory (PT) was good. The dipolar Yukawa (DY) and Stockmayer (STM) fluids are more reasonable models of polar fluids than DHS as in these models the non-electrostatic attractions are included. The models were widely used in the calculation of the dielectric constant [5–7], the internal energy and pressure [5,8,9], the heat capacity [10,11], phase equilibrium in an external field [12] and system with a high
∗ Corresponding author. Tel.: +48 61 82 91 454; fax: +48 61 82 91 505. E-mail address: [email protected] (S. Lamperski). 0378-3812/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2011.12.017
dipole moment [13,14]. Computer simulations of the free energy and activity coefficient have been rarely performed. The umbrella sampling [15] and the self-consistent histogram [16] methods were applied to calculate the free energy. The chemical potential of STM fluid was evaluated from the MC simulation in the Gibbs ensemble [17–19]. Here we apply the IGCMC method for the calculation of activity coefficients of DY and STM fluids. The advantage of IGCMC over the classical GCMC simulation [20] is that the former allows evaluation of the activity coefficient for a specified concentration (number density) of a system. A similar possibility is available by applying the algorithm propose by Malasics et al. [21]. The new technique has been successfully applied to the primitive model electrolyte [3], to the solvent primitive model electrolyte [22], and to DHS fluid [2]. The IGCMC mean activity coefficients of the model electrolyte were applied in the GCMC simulations of different inhomogeneous systems [23–25]. The new technique allows us to calculate not only the mean activity coefficient of an electrolyte, but also individual coefficients. For polar fluids, most of the thermodynamic functions had been studied earlier except for the activity coefficients. So, our interest is to fill this gap. Recently, we have considered DHS fluid [2]. In this paper, we analyze the activity coefficients of DY and STM potentials at different electric dipole moments and densities of a fluid, calculated by the IGCMC technique. Also, we present a comparison with the MSA and PT results.
30
S. Lamperski, R. Górniak / Fluid Phase Equilibria 317 (2012) 29–35
the point electric dipole moment is located at the centre of a hard sphere. The potential of the dipole–dipole interactions is
Nomenclature Add Add , Addd 2 3 Aex A* B b2 C d dT d2 G Idd , Iddd kB k1 N N0 n p r T u V y zdd zex ˛ ε0 ε εY C HS ex *
dipole contribution to the Helmholtz energy second and third perturbation term excess free energy reduced free energy term in the Ewald sum MSA parameter term in Ewald sum hard-sphere diameter temperature dependent hard-sphere diameter MSA parameter reciprocal lattice vector integrals in PT Boltzmann constant coefficient in IGCMC technique number of molecules in the simulation box number of molecules in the simulation box corresponding to the specified number density integer vector in Ewald sum pressure distance between two molecules absolute temperature potential energy volume dipole strength compression factor excess compression factor Gaussian distribution parameter vacuum permittivity the depth of well in the Lennard–Jones potential the depth of well in the Yukawa potential activity coefficient Coulombic contribution to the activity coefficient hard-sphere contribution to the activity coefficient packing fraction the range of attractions in the Yukawa potential electric dipole moment, chemical potential excess chemical potential number density reduced number density coupling parameter the separation at which the Lennard–Jones potential uLJ = 0
2. The model
uD (rij ) =
[i · j rij2 − 3(i · rij )(j · rij )] 4 ε0 rij5
(2)
where ε0 is the vacuum permittivity. Thus the DY potential, uDY , is the sum of uY and uD potentials. In the STM model, hard spheres are replaced by soft spheres. The point electric dipole moment is located at the centre of the soft sphere. The soft and dispersion interactions are given by the LJ potential
12 rij
uLJ (rij ) = 4ε
6
−
rij
(3)
where ε and are the LJ potential parameters (the depth of the potential well and the separation at which uLJ = 0, respectively). The STM potential, uSTM , consists of uLJ and uD . The LJ soft-core potential makes the STM model more realistic than DY. We use the following reduced properties for the DY potential: (∗ )2 =
2 , 4 ε0 εY d3
T∗ =
kB T , εY
∗ = d3
(4)
∗ = 3 .
(5)
and for the STM model: (∗ )2 =
2 , 4 ε0 ε 3
T∗ =
kB T , ε
where kB is the Boltzmann constant, T is temperature and is the number density of molecules. 3. Simulation methods and theories 3.1. The inverse grand-canonical Monte Carlo method The details of the GCMC simulation of polar fluids have been described for example in our previous paper [2]. The disadvantage of the GCMC technique is that it does not allow evaluation of the activity coefficient for a specified concentration (number density). This disadvantage is overcome in the IGCMC method [3]. In IGCMC, we divide the whole simulation into m short GCMC simulations in each of which the number n of configurations is generated. The simulation starts with a reasonable initial value of ln 0 . After each simulation the average number of molecules is calculated. If this number is greater than the number of molecules N0 in the simulation box corresponding to the specified number density , then a new value of ln i is obtained by subtracting from the previous ln i − 1 a small value k (e.g. 0.02). Otherwise ln i − 1 is increased by k. After the equilibration is reached, the average ln i over the m short simulations is calculated 1 ln i m m
We consider the dipolar Yukawa (DY) and Stockmayer (STM) fluids. The Yukawa potential, uY , is the sum of the hard-sphere and attractive (dispersion) interactions. For two molecules at positions ri and rj , uY is given by
uY (rij ) =
∞, εY d exp[−(rij − d)], − rij
rij < d rij ≥ d
(1)
where rij is the distance between the molecules i and j, d is the diameter of a hard sphere, and εY and are the depth and range of the attraction interactions, respectively. The properties of Yukawa fluid are similar to those of the Lennard–Jones (LJ) fluid when = 1.8/d [26] and this value for we will use in this work. In the DY model,
ln =
(6)
i=1
3.2. Thermodynamics The activity coefficient is related to the excess chemical potential ex by the formula ex = kB T ln
(7)
The excess chemical potential can be obtained from the excess free energy, Aex , using the thermodynamic relation
ex =
∂Aex ∂N
(8) V,T
S. Lamperski, R. Górniak / Fluid Phase Equilibria 317 (2012) 29–35
Let us introduce the reduced free energy, A* , per molecule A∗ =
and
A NkB T
(9)
Using Eqs. (7)–(9) we obtain
1 ln = kB T
∂Aex ∂N
=
A∗ex
+N
V,T
∂A∗ex ∂N
∂A∗ex ∂N
= V,T
∂A∗ex ∂
= ∗ T
∂A∗ex ∂∗
(10) V,T
= T
∂A∗ex ∂
(11) T
where is the packing fraction =
∗ 6
(12)
The activity coefficient can be also calculated from the excess contribution to the compression factor, zex (the excess against the perfect gas). The thermodynamic relation, which tells that pressure is a negative value of the derivative of free energy against volume, can be written in terms of the derivative of the number density
−p =
∂A ∂V
∂A ∂
= N,T
T
∂ ∂V
2 =− N
N
∂A ∂
(13) T
∂ ∂V
=− N
2 N
(14)
From the equality of the first and the last term of Eq. (13), and using Eq. (9) we obtain the expression for A∗ex
A∗ex
zex d =
= 0
∗
zex ∗ d = ∗
0
0
zex d
∂A∗ex ∂
ln =
(16)
dT =
zex d + zex
(17)
APT,dd DHS NkB T
=
Add /NkB T 2 1 − Addd /Add 3 2
(18)
Here Add 2 NkB T Addd 3 NkB T
1 − exp
−uLJ (rij )
kB T
drij
(23)
Cotterman et al. [31] approximated this equation by the expression which for the spherical LJ potential takes the form dT =
1 + 0.29770T ∗ 1 + 0.33163T ∗ + 0.0010477T ∗2
rmin
0 rmin
rij2 yHS (rij ) exp
rij2 yHS (rij ) exp
−
(24)
−uHS (rij )
kB T
drij
−(uLJ (rij ) + ε)
drij = 0
kB T
(25)
0
uHS (rij ) =
∞, 0,
rij < d rij ≥ d
(26)
and yHS is the Percus–Yevick y function [33–35] yHS (rij ) = 1 + 62
r ij
d
+ 12 1
3 rij d
(27)
with 1 =
In the PT theory of DY fluid, the Yukawa potential is considered as the reference potential and the dipole–dipole interactions as a perturbation. The free energy of the dipole–dipole interactions, Add , DY can be approximated by the corresponding free energy Add of the DHS DHS fluid. Henderson et al. [27] obtained the analytical solution of PT with the Padé approximation given by Stell et al. [28]
NkB T
0
3.3. The perturbation theory
≈
T
0
APT,dd DY
Using the above two results, the expression (10) for the activity coefficient written in terms of zex takes the form
where rmin = 21/6 , uHS is the hard-sphere potential
= zex
(22)
When APT,dd is known, the activity coefficient can be calculated DY from Eq. (10) with the numerically evaluated derivative. In the STM fluid the LJ potential is treated as the reference system, and the dipole–dipole interaction as a perturbation [28,29]. According to Barker and Henderson (BH) [30] the free energy of a soft-sphere fluid can be approximated by that of a hard-sphere fluid with the temperature dependent hard-sphere diameter, dT . For the LJ potential the Authors suggested calculating this diameter from the equation
(15)
and for its derivative
5(1 + 1.12754d3 + 0.561922 d6 ) 24(1 − 0.05495d3 + 0.133322 d6 )
(21)
and is valid for 0 < T* < 15. Another approach for dT was proposed by Weeks, Chandler, and Anderson (WCA) [32]. It requires solution of the equation
as at N hold constant
1 + 0.18158d3 − 0.114672 d6 3(1 − 0.49303d3 + 0.062932 d6 )
I dd =
I ddd =
The last term of Eq. (10) can be rewritten in terms of , * , or N
31
2 =
(1 + 2)2
(28)
(1 − )4 −(1 + /2)
2
(1 − )4
(29)
The diameter dT calculated from Eq. (25) is both temperature and density dependent. Now the dipolar free energy of STM, APT,dd , can be obtained from STM Eq. (18) where the hard-sphere d in Eqs. (19)–(22) is replaced with dT calculated from Eq. (23) or (25). The temperature dependence of diameter means that * and also depend on temperature, but does not. 3.4. The MSA theory
=−
=
2 4 3d3 (4 ε0 kB T )
2
4 2 2 6 27d3 (4 ε0 kB T )
3
I
dd
(19)
The MSA free energy of DY fluid, AMSA , can be separated into the DY
hard-sphere Yukawa AMSA and dipolar AMSA,dd contributions [36] Y DY I
ddd
(20)
MSA AMSA + AMSA,dd DY = AY DY
(30)
32
S. Lamperski, R. Górniak / Fluid Phase Equilibria 317 (2012) 29–35
The dipolar contribution AMSA,dd is approximated by the dipoDY lar hard-sphere free energy,
MSA, dd ADHS .
7
The analytical solution for
6
was given by Wertheim [37] (the free energy route) AMSA,dd DHS AMSA,dd DY NkB T
AMSA,dd DHS
≈
NkB T
3 = − I(y)
5 (31)
y=
ln
where y is the dipole strength 2
(32)
9ε0 kB T
8 2 3
(1 + )
2
(1 − 2)
+
4
(2 − )
2
8(1 + )
(33)
4
(1 + 4)
2
(1 − 2)
4
−
(1 − 2) (1 + )
2
4
= 3y
(34)
3 3 I(y) − y
(35)
, This expression is obtained by differentiating the free energy, AMSA DHS so the ln dd results obtained from the free energy (Eq. (31)) and from the above formula are exactly the same. Another expression for the compression factor is from the virial (pressure) route [38]. It gives dd zDY, v
≈
dd zDHS, v
3 = − y
(36)
As we see, different routes give different results for the compression factor. This is the example of the so-called thermodynamic inconsistency which stems from the fact that the theory is approximate. An interesting MSA approach based on the perturbation theory was given by Kronome, Liszi and Szalai (KLS) [39]. By adding Eqs. (35) and (36) with weights of 2/3 and 1/3, respectively, like Carnahan and Starling [40] did, the Authors proposed the following expression for the dipolar contribution to the compression factor: dd dd zDY,KLS ≈ zDHS,KLS =
2 3 I(y) − y
(37)
In this work, we use the above three approaches to calculate the activity coefficient from the compression factor using Eq. (17). As there is no analytical solution of MSA for the STM fluid [5] dd from Eqs. (31) and (37), where one can evaluate AMSA,dd and zSTM STM is replaced with the temperature dependent T as described in Section 3.3. 3.5. Long-range interactions The long-range dipolar interactions in IGCMC were handled with the Ewald sum (ES) [41] methods. In the ES method, the total potential energy of dipole–dipole interactions consists of a real space sum (the first term of Eq. (38)), reciprocal space sum (the second term) reduced by the self-interaction term (the third term) 1 1 uij = 3 2 rij + n i,j=1 n N
/ 0 i,j=1G =
0.0
0.2
0.4
0.6
0.8
Fig. 1. Dependence of the IGCMC activity coefficient ln on * for the dipolar Yukawa fluid at T* = 3.0 and for the following (* )2 : 1.0 (circles), 3.0 (squares), and 5.0 (triangles). The line through the IGCMC points is given as a visual aid.
2˛r B(r) = erfc(˛r) + √ exp(−˛2 r 2 )
(39)
2˛r C(r) = erfc(˛r) + √ exp(−˛2 r 2 )(1 + 23 ˛2 r 2 )
(40)
where n = (nx , ny , nz ), parameter ˛ regulates the relative convergence of the sum. Reciprocal lattice vectors are given by G = ((2 /L)nx , (2 /L)ny , (2 /L)nz ). Eqs. (38)–(40) are taken from [42]. 4. Results and discussion The IGCMC simulations for DY and STM fluids were carried out for a reduced temperature T* = 3.0, reduced electric dipole moment (* )2 = 1.0, 3.0, and 5.0, and reduced density * ranging from 0.1 to 0.8. The temperature T* is higher than the critical temperature Tc∗ for STM fluid which for the values of (* )2 considered is 1.41 [17], 1.82 [18], and 2.31 [19], respectively. For T = 298.15 K and d = = 400 pm, ε = εY = 826.33 J/mol, while the dipole moment in D corresponding to the considered values of (* )2 is 0.9371, 1.6231, and 2.0955. Earlier, we have found [2] that the activity coefficient results for dipolar hard-sphere fluid are weakly dependent on the number of molecules, N0 . Thus in the present simulations we used N0 = 256. The ES method was applied to describe the long-range dipole–dipole interactions. The value of parameter ˛ was 5.1. The sum in the reciprocal space extended over all values of G such that |n|2 ≤ 27. The number of simulation steps applied to equilibrate the system was 6,000,000 while in the main part of the simulation this number amounted to 10,000,000. To reach faster a value of ln close to the expected one, we applied k1 = 0.05 in the first part of the simulation, while in the second part k1 = 0.02. The standard deviation of ln depended on * and varied from 0.02 to 0.03. The IGCMC results for the DY fluid are shown in Fig. 1 and listed in Table 1. With increasing value of (* )2 the activity coefficient decreases. The curve for (* )2 = 1.0 has the exponential-like shape similar to that for a the hard-sphere fluid [22]. For higher (* )2 a shallow minimum with negative values is observed at low * .
i · (rij + n) · j · (rij + n) i · j B ˛ rij + n − 3 C ˛ rij + n rij + n2
2 −G2 /(4˛2 ) (i · G)(j · G) 2˛3 + e exp(iG · r ) − √ ij V G2 3
N
0
According to Eq. (17) the activity coefficient ln dd can be also calculated from the dipolar contribution to the compression factor, zdd . In the free energy route [38] dd dd zDY,A ≈ zDHS,A =
2
-1
The parameter is obtained by numerical solution of the following equation:
3
1
and the function I(y) is given by I(y) =
4
i=1....N
2i
(38)
S. Lamperski, R. Górniak / Fluid Phase Equilibria 317 (2012) 29–35 Table 1 The IGCMC activity coefficient ln for the dipolar Yukawa fluid at T* = 3.0. *
4
ln
3
(* )2 = 0.0
(* )2 = 1.0
(* )2 = 3.0
(* )2 = 5.0
0.06768 0.2124 0.4643 0.8700 1.5016 2.4607 3.9215 6.0987
0.05006 0.1748 0.4058 0.7843 1.3919 2.30869 3.7566 5.9270
−0.09116 −0.09450 0.00972 0.2551 0.7157 1.4861 2.7280 4.7140
−0.3718 −0.5751 −0.6373 −0.5452 −0.2508 0.3386 1.4532 3.3095
2
ln
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
33
1 0 -1 -2
0.0
0.0
0.2
0.4
0.6
0.8
-0.5 Fig. 4. Dependence of the IGCMC activity coefficient ln on * for the Stockmayer fluid at T* = 3.0 and for the following (* )2 : 1.0 (circles), 3.0 (squares), and 5.0 (triangles). The line through the IGCMC points is given as a visual aid.
ln
dd
-1.0 -1.5 -2.0 -2.5 -3.0 0.0
0.2
0.4
0.6
0.8
Fig. 2. Dependence of the dipolar contribution to the activity coefficient ln dd on * for the dipolar Yukawa fluid. Physical parameters and symbols as in Fig. 2.
The activity coefficient of DY can be separated, like the free energy, into the Yukawa ln Y and dipolar ln dd contributions. The ln dd results are shown in Fig. 2. These results take negative values and a nearly linear dependence of ln dd on * is observed. In Fig. 3, we compare the IGCMC dipolar contribution (triangles- the solid line is given as a visual aid) with the corresponding PT (a dashed line) and MSA (short dashed, dotted-dashed, and dotted-dotteddashed lines) data for (* )2 = 5.0. The ln dd results calculated by the MSA virial route (a dotted-dashed line) are evidently overestimated. Better agreement with the computer simulation results
0
ln
dd
-1
-2
is observed for ln dd obtained from the MSA free energy route (a short dashed line) although this time the theoretical results are underestimated. The curve following from the KLS approximation (a dotted–dotted-dashed line) runs close to the IGCMC line at * < 0.5 while at high * some divergence is observed which intensifies with increasing * . The PT results are in excellent agreement with the IGCMC data in the whole range of * . Earlier, for a dipolar hard-sphere fluid, we also observed [2] a great similarity of PT and IGCMC results. In Fig. 4 and in Table 2, the IGCMC results for the STM fluid are presented. The curves shown in Fig. 4 are similar to those in Fig. 1 although the value of the activity coefficient in Fig. 4 is lower and the minimum of STM curves is deeper. It is due to a stronger dispersion attraction in LJ than in Y fluid (compare the second columns of Tables 1 and 2). The ln dd curves shown in Fig. 5 exhibit similar behaviour like DY curves but they are slightly more bent. In the next two graphs, we compare the ln dd data for the STM fluid obtained by the IGCMC simulation with the PT and MSA results for (* )2 = 5.0. Fig. 6 shows the theoretical results for dT calculated from the BH formula (Eq. (23)) which for T* = 3.0 yields dT / = 0.94505. Again the MSA virial and free energy route curves run far from the simulation data. The KLS curve is in much better agreement with the IGCMC results but the curvatures of these curves are different. The PT curve is very close to the IGCMC one in the whole range of * . The value of dT / calculated from the WCA theory (Eq. (25)) depends on the density and for T* = 3.0 it varies from 0.97312 for * = 0.10 to 0.96900 for * = 0.80. The MSA and PT activity coefficients calculated with theses values of dT are shown in Fig. 7. As seen, the WCA theory does not improve the convergence of the MSA virial and free energy route results with the simulation data. However, the KLS approach at * < 0.4 and the PT theory at * > 0.4 give nearly perfect agreement with the IGCMC results. Table 2 The IGCMC activity coefficient ln for the Stockmayer fluid at T* = 3.0.
-3
*
-4 0.0
0.2
0.4
0.6
0.8
Fig. 3. Dependence of the dipolar contribution to the activity coefficient ln dd on * for the dipolar Yukawa fluid at T* = 3.0, and (* )2 : = 5.0. The symbols represent the IGCMC data while the lines the different theories results (PT – dashed line, MSA (free energy route) – short dashed line, MSA (virial route) – dotted–dashed line, and KLS – dotted–dotted–dashed line).
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
ln (* )2 = 0.0
(* )2 = 1.0
(* )2 = 3.0
(* )2 = 5.0
−0.02312 0.01136 0.12181 0.3381 0.7059 1.2886 2.1796 3.4833
−0.04667 −0.03467 0.05276 0.2431 0.5855 1.1394 1.9774 3.2424
−0.2077 −0.3368 −0.3818 −0.3239 −0.1310 0.2616 0.9345 2.0245
−0.5497 −0.8999 −1.1150 −1.2217 −1.1897 −0.9684 −0.4798 0.4017
34
S. Lamperski, R. Górniak / Fluid Phase Equilibria 317 (2012) 29–35
0.0
0
-0.5
dd
-1
-1.5
ln
ln
dd
-1.0
-2.0
-2
-2.5 -3.0
-3
-3.5 0.0
0.2
0.4
0.6
0.8
Fig. 5. Dependence of the dipolar contribution to the activity coefficient ln dd on * for the Stockmayer fluid. Physical parameters and symbols as in Fig. 4.
0.0
0.2
0.4
0.6
0.8
Fig. 8. Dependence of the IGCMC dipolar contribution to the activity coefficient ln dd on * at T* = 3.0, and (* )2 = 5.0 for the dipolar hard-sphere (circles, solid line), dipolar Yukawa (squares, dashed line) and Stockmayer (triangles, short dashed line) fluids.
0
ln
dd
-1
-2
-3
-4 0.0
0.2
0.4
0.6
0.8
Fig. 6. Dependence of the dipolar contribution to the activity coefficient ln dd on * for the Stockmayer fluid at T* = 3.0, and (* )2 : = 5.0. The symbols represent the IGCMC data while the lines the different theories results with dT given by Eq. (23) (PT – dashed line, MSA (free energy route) – short dashed line, MSA (virial route) – dotted–dashed line, and KLS – dotted–dotted–dashed line).
Although the DHS fluid is not a subject of this paper, it is interesting to compare the properties of the following dipolar fluids: DHS, DY, and STM. For theses models, the IGCMC results of ln dd at *2 = 5.0 are collected in Fig. 8. The DHS and DY results are nearly identical. It is not surprising as in both cases the repulsion interactions are described by the same hard-sphere potential. The STM curve runs below the DHS and DY lines. As follows from the BH formula or WCA theory, the average closest approach of two STM (or LJ) molecules is smaller than . Thus for = d, as we assumed in this paper, the point dipole moments of STM molecules can approach closer to each other than the dipoles of DHS or DY fluids and therefore the dipole–dipole interactions are stronger. This lowers ln dd of STM in comparison to DHS or DY. A separate problem is whether the dispersion terms in DY and STM potentials, which give rise to the intermolecular attractions, influence the ln dd results. An indepth analysis of the DHS and DY curves implies the positive answer to this question. At higher densities a subtle but permanent separation of the curves is seen with the DY curve running below the DHS one. As dispersion is a short range interaction, it should be manifested at a high density, as shown in Fig. 8. However, the separation is small and can be comprised in the error bar, so the problem needs further, much more precise investigation.
0
5. Conclusions
ln
dd
-1
-2
-3
-4 0.0
0.2
0.4
0.6
0.8
Fig. 7. Dependence of the dipolar contribution to the activity coefficient ln dd on * for the Stockmayer fluid at T* = 3.0, and (* )2 : = 5.0. The symbols represent the IGCMC data while the lines the different theories results with dT given by Eq. (25) (PT – dashed line, MSA (free energy route) – short dashed line, MSA (virial route) – dotted–dashed line, and KLS – dotted–dotted–dashed line).
In this paper, we have shown that the IGCMC method can be successfully applied to calculate the activity coefficient of the DY and STM fluids. Subtraction of the activity coefficients of the Y and LJ fluids from these results gave the dipolar contributions to the overall activity coefficient. Comparison of the dipolar contribution yielded by the IGCMC simulation with that obtained by the PT, MSA and KLS theories showed that for the DY fluid the PT results are in a very good agreement with the simulation data in the whole range of * . At low * the KLS approach also gives the correct results. For the STM fluid with dT calculated according to the WCA theory, the KLS approach at * < 0.4 and the PT theory at * > 0.4 give nearly a perfect agreement with the IGCMC results. Finally, we should note that our conclusion that the MSA virial route underestimates the magnitude of the activity coefficient while PT is much more accurate is consistent with the results for other thermodynamic properties. A great similarity of the DHS and DY results shown in Fig. 8 supports the assumption that the free energy of the dipole–dipole interactions of the DY fluid can be approximated by the
S. Lamperski, R. Górniak / Fluid Phase Equilibria 317 (2012) 29–35
corresponding free energy of DHS fluid (Eqs. (18) and (31)). Explanation of a subtle difference observed at * > 0.4 requires low temperature investigation which will better separate the dipolar and dispersion attraction from the temperature independent hardsphere interactions. It could throw light on the role of dispersion attraction in the dipole–dipole interaction of DY fluid. Acknowledgments
[14] [15] [16] [17] [18] [19] [20] [21] [22] [23]
Financial support from Adam Mickiewicz University, Faculty of Chemistry, is appreciated.
[24] [25] [26]
References
[27]
[1] D. Frenkel, B. Smit, Understanding Molecular Simulations: From Algorithms to Applications, Academic Press, San Diego, 1996, p. 121. [2] S. Lamperski, R. Górniak, Fluid Phase Equilib. 295 (2010) 255–263. [3] S. Lamperski, Mol. Simul. 33 (2007) 1193–1198. [4] B. Widom, J. Chem. Phys. 39 (1963) 2808–2812. [5] I. Szalai, D. Henderson, D. Boda, K.-Y. Chan, J. Chem. Phys. 111 (1999) 337–344. [6] I. Szalai, K.-Y. Chan, Y.W. Tang, Mol. Phys. 101 (2003) 1819–1828. [7] L.E. Johnson, R. Barnes, T.W. Draxler, B.E. Eichinger, B.H. Robinson, J. Phys. Chem. B 114 (2010) 8431–8440. [8] C. Kriebel, Winkelmann, J. Mol. Phys. 88 (1996) 559–578. [9] Y.V. Kalyuzhnyi, I.A. Protsykevytch, P.T. Cummings, Europhys. Lett. 80 (2007) 56002. [10] Z. Máté, I. Szalai, Fluid Phase Equilib. 54 (2010) 54–60. [11] Z. Máté, I. Szalai, D. Boda, D. Henderson, Mol. Phys. 109 (2011) 203–208. [12] A. Gábor, I. Szalai, Mol. Phys. 106 (2008) 1094–1106. [13] Yu.V. Kalyuzhnyi, Y.V. Protsykevitch, Condens. Mater. Phys. 10 (2007) 553–562.
[28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42]
35
Yu.V. Kalyuzhnyi, G. Stell, Mol. Phys. 78 (1993) 1247–1258. G.M. Torrie, J.P. Valleau, J. Comp. Phys. 23 (1977) 187–199. C.H. Bennet, J. Comp. Phys. 22 (1976) 245–268. B. Smit, C.P. Williams, E.M. Hendriks, S.W. de Leeuw, Mol. Phys. 68 (1989) 765–769. M.E. van Leeuwen, B. Smit, E.M. Hendriks, Mol. Phys. 78 (1993) 271–283. M.E. van Leeuwen, Mol. Phys. 82 (1994) 383–392. J.P. Valleau, L.K. Cohen, J. Chem. Phys. 72 (1980) 5935–5942. A. Malasics, D. Gillespie, D. Boda, J. Chem. Phys. 128 (2008) 124102–124108. S. Lamparski, M. Płuciennik, Mol. Simul. 36 (2010) 111–117. S. Lamperski, C.W. Outhwaite, L.B. Bhuiyan, J. Chem. Phys. B 113 (2009) 8925–8929. S. Lamperski, C.W. Outhwaite, J. Colloid Interface Sci. 328 (2008) 458–462. S. Lamparski, J. Kłos, J. Chem. Phys. 129 (2008) 164503–164510. D. Henderson, E. Waisman, J.L. Lebowitz, L. Blum, Mol. Phys. 35 (1978) 241–255. D. Henderson, L. Blum, A. Tani, Advances in Chemistry Series 300, American Chemical Society, Washington, DC, 1986, pp. 281–296. G. Stell, J.C. Rasaiah, H. Narang, Mol. Phys. 27 (1974) 559–564. G. Stell, J.C. Rasaiah, H. Narang, Mol. Phys. 27 (1974) 1393–1414. J.A. Barker, D. Henderson, J. Chem. Phys. 47 (1967) 2856–2862. R.L. Cotterman, B.J. Schwarz, J.M. Prausnitz, AIChE J. 32 (1986) 1787–1798. J.D. Weeks, D. Chandler, H.C. Anderson, J. Chem. Phys. 54 (1971) 5237–5247. M.S. Wertheim, Phys. Rev. Lett. 10 (1963) 321–323. M.S. Wertheim, J. Math. Phys. 5 (1964) 643–652. E. Thiele, J. Chem. Phys. 39 (1963) 474–480. D. Henderson, D. Boda, I. Szalai, K.Y. Chan, J. Chem. Phys. 110 (1999) 7348–7358. S.M. Wertheim, J. Chem. Phys. 55 (1971) 4291–4298. C.G. Gray, K.E. Gibbson, Theory of molecular fluids Fundamentals, vol. 1, Oxford University Press, New York, 1984, p372. G. Kronome, J. Liszi, I. Szalai, J. Chem. Soc. Faraday Trans. 93 (1997) 3053–3059. F. Carnahan, E. Starling, J. Chem. Phys. 53 (1969) 635–637. P. Ewald, Ann. Phys. 64 (1921) 253–287. W. Smith, CCP5 Newsletter 46 (1998) 18–23.