Accurate equations of state for CF4, CF4–Ar, and CF4–CH4 fluids using two-body and three-body intermolecular potentials from molecular dynamics simulation

Journal of Fluorine Chemistry 168 (2014) 81–92

Contents lists available at ScienceDirect

Journal of Fluorine Chemistry journal homepage: www.elsevier.com/locate/fluor

Accurate equations of state for CF4, CF4–Ar, and CF4–CH4 fluids using two-body and three-body intermolecular potentials from molecular dynamics simulation Mohsen Abbaspour *, Maryam Sheykh Department of Chemistry, Hakim Sabzevari University, Sabzevar, Khorasan Razavi, Iran

A R T I C L E I N F O

A B S T R A C T

Article history: Received 22 May 2014 Received in revised form 17 August 2014 Accepted 7 September 2014 Available online 16 September 2014

Molecular dynamics simulations have been performed to obtain pressures and equations of state of CF4, CF4–Ar, and CF4–CH4 fluids using different inversion and ab initio pair-potentials. To take many-body forces into account, the three-body potentials of Hauschild and Prausnitz, Mol. Simul. 11 (1993) 177–185, Wang and Sadus, J. Chem. Phys. 125 (2006) 144509–144513, and Guzman et al., Mol. Phys. 109 (2011) 955–967 have been used with the pair-potentials. The significance of this work is that the many-body potential of Hauschild and Prausnitz is extended as a function of density, temperature, and molar fraction and is used with the HFD-like pair-potentials of CF4, CF4–Ar, and CF4–CH4 systems to improve the prediction of the pressure values without requiring an expensive three-body calculation. We have also simulated the self-diffusion coefficient of CF4 in good agreement with experimental data. ß 2014 Elsevier B.V. All rights reserved.

Keywords: CF4 Many-body interaction Molecular dynamics simulation Equation of state Self-diffusion coefficient

1. Introduction Theoretical calculations and computer simulations, usually applied for the prediction of bulk properties, still prefer small, highly symmetric molecules, like CF4, for their investigations. Carbon tetrafluoride, one of the halocarbon refrigerants is a rather simple, tetrahedral molecule with octopol–octopol interactions and represents a model system for several theoretical studies. It is also used in electronics microfabrication alone or in combination with other gases as a plasma etchant for silicon, silicon dioxide, and silicon nitride [1,2]. The practical applicability of molecular simulations in process engineering requires appropriate molecular models for real fluids. The molecular models of the individual fluids and fluid mixtures must accurately describe thermodynamic data of them [3]. For carbon tetrafluoride, several intermolecular potentials are available in the literature [4–12]. Some of them fit ab initio data (including basis set effects) to an analytical functional form [7,9–12] and other works use one or two center Lennard-Jones (LJ) interaction model (including quadrupole or octopole moments) on each of the interaction sites parametrized on the basis

* Corresponding author. Tel.: +98 5144013322. E-mail addresses: [email protected], [email protected] (M. Abbaspour). http://dx.doi.org/10.1016/j.jfluchem.2014.09.005 0022-1139/ß 2014 Elsevier B.V. All rights reserved.

of experimental data for second virial coefficient, self-diffusion coefficient, or vapor–liquid equilibrium data. Goharshadi et al. [13] determined the potential energy function of CF4–CF4 system via the inversion of reduced viscosity collision integrals at zero pressure and fitted to obtain a Hartree–Fock dispersion (HFD)-like potential which accurately produced second virial coefficient and transport properties of carbon tetrafluoride over wide temperature ranges. Therefore the two-body HFD-like potential of Goharshadi et al. [13] has been used in this work. Recently, Vayner et al. [14] performed ab initio calculations to characterize the CF4–Ar intermolecular potential. They investigated both the size of the basis set and the basis set superposition error (BSSE) and proposed analytical expressions as a sum of two-body potentials. We have also used the ab initio potential of Vayner et al. [14] in our simulation. Clifford et al. [15] determined an effective spherically averaged intermolecular pair potential for CF4–CH4 system using inversion technique from transport properties. We have also used the inversion potential of Clifford et al. [15] in our simulation. Besides the two-body interactions, it is also well known [16,17] that three-body interactions can make a small but significant contribution to the energy of fluids such as rare gases, diatomic and polyatomic molecules. To obtain a quantitative agreement with the experiment, pair potentials must be used in conjunction with three-body interactions. Hauschild and Prausnitz (HP) [18] have proposed a many-body correction term to the true two-body

M. Abbaspour, M. Sheykh / Journal of Fluorine Chemistry 168 (2014) 81–92

82

viscosity, self-diffusion coefficient, and second virial coefficient of CF4 in very good agreement with experimental data over a wide temperature range. Therefore, the pair-potentials of Goharshadi et al. [13] has been selected for simulations of the systems in this work. It should be also noted that these inversion potentials have been obtained at short and intermediate separation ranges (and not in the whole range). Knowing the inner branch of the potential well from the viscosity in conjunction with the second virial coefficient helps to determine the outer branch of the potential (long range) [26]. The equations used for this purpose are as follows:

potential to evaluate phase behavior of the fluids such as methane and argon using simulation without requiring an expensive threebody calculation. Wang and Sadus (WS) [19] also showed that there is a simple and accurate relationship between the two-body and three-body interactions that allows us to obtain the threebody effect without any additional computational cost. Guzman et al. (GZ) [20] also developed a model for an effective triple-dipole interaction for some monatomic and polyatomic molecules based on accurate calculations of three-body effects on the third virial coefficient. The thermodynamics and transport properties of the heavy globular molecules such as CF4 (and its mixtures) are strongly requested by contemporary technologies (chemical vapor deposition, for example) for optimization of the processes and design of the devices. The experimental determination of these toxic and aggressive gases in a wide range of temperatures is complicated and expensive [21]. Therefore, the purpose of the present paper is to perform Molecular dynamics (MD) simulations to obtain pressures and equations of states of CF4, CF4–Ar, and CF4–CH4 fluids using the different pair-potentials. To take many-body forces into account, the HP [18], WS [19], and GZ [20] three-body potentials have been used with the pair-potentials. We have also simulated the self-diffusion coefficient of CF4 at different densities.

U 2B

e

¼ T  1

(1)

rR3  rL3 ¼ ðB  BHS ÞNðT  Þ

(2)

where rR and rL are the coordinates of the outer and inner wall of the potential well, respectively. B is the experimental second virial coefficient [27] and BHS is the second virial coefficient of a hard sphere, namely: 2pNAs3/3. The values of N(T*) has been given in Ref. [26]. At long range, only the well width of the potential obtained from the second virial coefficient data is available. This has been used in conjunction with the inner coordinates of the well obtained in the viscosity inversion to give the potential energy in the whole separation range. Therefore, the isotropic (spherically averaged) one-site pairpotentials of CF4–CF4, Ar–Ar, and CH4–CH4 systems, which have been obtained using the inversion method and Eqs. (1) and (2) at the whole separation ranges, have been presented in Figs. 1 and 2 and fitted to the following HFD-like model which have been used in our simulations:

2. Theory 2.1. Intermolecular pair-potentials The monatomic and some of polyatomic Systems have been intensively and successfully studied over a broad range of temperatures and densities using pair interactions obtained from inversion method [13,22–24]. Goharshadi et al. [13,22,23] determined the potential energy functions of CF4–CF4, Ar–Ar, and CH4–CH4 systems via the inversion of reduced viscosity collision integrals at zero pressure. The aim of an inversion method is to obtain the potential function by considering the experimental data as a functional form instead of fitting the data to a constrained potential form having a few parameters. The direct inversion procedure for the viscosity is based on the idea that at a given reduced temperature (T*) the value of reduced collision integral (V(2,2)*) is determined by the potential over only a small range of separation distance around a value of r. It should be mentioned that for non-spherical molecules, the presence of the internal structure leads to the existence of inelastic collisions and therefore, the inversion scheme produces an isotropic pair potential energy and so the collision integrals must be averaged over all possible relative orientations occurring in collisions [13,25]. The isotropic (spherically averaged) and one site potentials of Goharshadi et al. [13,22,23] for CF4–CF4, Ar–Ar, and CH4–CH4 systems from the inversion procedure have been reproduced the

   C  C  C10   ¼ A exp a y þ b y2  66  88  10 U2B y y y

(3)

 where U2B is the reduced two-body (pair) potential,   U2B ¼ U 2B =eU2B , and y = r/s (s is the distance at which the intermolecular potential has zero value and e is the potential well depth). The values of the parameters of the pair-potentials have been given in Table 1. Borodin et al. [10] developed an ab initio quantum chemistry based atom–atom force field for perfluoroalkanes performed at the MP2/aug-cc-pvDz level for different orientations. In the classical force field for CF4, the total potential energy of the ensemble of atoms, represented by the coordinate vector r, is denoted as UBORD(r). The latter is represented as a sum of nonbonded interactions UNB(rij) as well as energy contributions due to bonds UBOND(rij) having bond length rij, and bends UBEND(qijk) having bending angle uijk which is given by:  X NB   X BOND   X BEND  (4) U ri j þ U ri j þ U qi jk U BORD ðr Þ ¼ i< j

ij

i jk

Table 1 The coefficients of the one-site (spherically averaged) two-body HFD-like potentials for different interactions. Interaction CF4–CF4a CF4–CF4b CF4–CF4c Ar–Ar CH4–CH4 CF4–Ar CF4–CH4 a b c

e/k (K) 258.545 292.5 152.1 143.2240 161.3997 187.54 199.9976

Tsuzuki et al. pot. [11]. Mahlanen et al. pot. [12]. Goharshadi et al. pot. [13].

s (A˚) 4.248 3.8852 4.7 3.7565 3.721 3.81 3.975

A* 5

1.4039  10 3.86401  1013 1.4039  108 1.6167  107 1.4039  105 1.4039  105 5.3974  104

a*

b*

C6 *

C8*

C10*

4.0000 26.0780 40.0286 57.4722 4.0000 4.0000 0.9880

50.0000 6.6486 73.8642 305.0846 50.0000 50.0000 10.9924

8.7540 5.2820 1.8518 1.8359 2.2239 4.7700 4.4845

34.0740 25.0930 5.8184 1.9973 6.5193 23.2300 5.2494

25.1590 19.9600 7.1632 3.1696 8.7242 18.5300 3.0178

M. Abbaspour, M. Sheykh / Journal of Fluorine Chemistry 168 (2014) 81–92

83

Fig. 1. The different reduced two-body and total (two-body plus three-body) potentials used in the fluid CF4 simulations. All of the potentials have been reduced by e/ k = 152.1 K and s = 4.7 A˚.

84

M. Abbaspour, M. Sheykh / Journal of Fluorine Chemistry 168 (2014) 81–92

Fig. 2. The different reduced two-body and total potentials used in the fluid CF4–Ar and. CF4–CH4 simulations. All of the potentials have been reduced by e/k = 152.1 K and s = 4.7 A˚. (The meaning of the different lines and symbols are the same as Fig. 1 with the exception of the solid lines described in the plots.)

M. Abbaspour, M. Sheykh / Journal of Fluorine Chemistry 168 (2014) 81–92

85

Table 2 The coefficients of the atom–atom pair-potential of Borodin et al. [10]. Interaction

A (kj/mol)

B (A˚1)

C (kj A˚6/mol)

KBOND (kj A˚2/mol)

r0 (A˚)

F–F C–F C–C F–C–F

265,994.300 28,938.976 62,599.680

4.261 3.084 3.090

519.574 1320.044 2678.544

3017.960

1.3391

The nonbonded energy UNB(rij) consists of the sum of two-body (pair) repulsion and dispersion energy terms between atoms i and j having atom types a and b represented by the Buckingham potential:     C ab U NB r i j ¼ Aab exp Bab r i j  6 ri j

(5)

whereas contributions due to bonds UBOND(rij) and bends UBEND(uijk) for atoms i, j, and k having atom types a, b, and g, respectively, are given by:  2   1 0 r i j  rab U BOND r i j ¼ kBOND 2 ab

(6)

 2   1 ui jk  u0abg U BEND ui j ¼ kBEND 2 abg

(7)

0 is an equilibrium bond length for atom types a and b; where rab u0abg is an equilibrium bend angle for atom types a, b, and g; and kabBOND and kabgBEND are the bond and bend force constants, respectively. The values of the parameters of the pair-potentials have been given in Table 2. We have also presented the nonbonded ab initio atom–atom interaction potential of Borodin et al. [10] in Fig. 1 and used in our simulations for pure CF4. Tsuzuki et al. [11] calculated intermolecular interaction energies of different orientation CF4 dimers using high level ab initio calculations at the MP2/aug-cc-pvTz level. They concluded that dispersion interaction is mainly responsible for the attraction between the CF4 dimers. We have fitted the HFD-like potential (Eq. (3)) to the calculated interaction energies at eight different orientations and presented in Fig. 1 and used in our simulations. The values of the resulted one-site HFD-like model have been given in Table 1. According to the work of Bock et al. [28], in the fitting procedure, the points are weighted with the following function:  2   E0  E E0 (8) þ ½1  QðE0  EÞ wðEÞ ¼ QðE0  EÞexp 2E0 E

where Q(x) is Heaviside’s step function. The parameter E0 has been chosen to +200 K [28]. The weight function (Eq. (8)) gives the points in the attractive region a higher weight. Mahlanen et al. [12] determined potential energy surface of CF4 dimers with MP2/aug(df)-6-311G* ab initio calculations for 11 different orientations. We have fitted the HFD-like potential (Eq. (3)) to the weighted interaction energies at different orientations using the exponential weight function of Bock et al. [28] and presented in Fig. 1. The values of the resulted one-site HFD-like model have been given in Table 1 and used in our simulations for pure CF4. Vayner et al. [14] obtained the CF4–Ar interactions using ab initio calculations at CCSD(T)/aug-cc-pVTZ level including BSSE correction. We have also fitted the HFD-like model (Eq. (3)) to the weighted interaction energies at different orientations using the exponential function of Bock et al. [28] (Eq. (8)) and presented in Fig. 2. The values of the resulted (spherically averaged) one-site HFD-like model of Vayner et al. [14] have been also given in Table 1 and used in our simulations for the CF4–Ar two-body interactions.

KBEND (kj rad2/mol)

u0 (deg)

1003.2

108.54

Clifford et al. [15] determined an effective spherically averaged intermolecular pair potential for CF4–CH4 system using the inversion technique. In order to use the spherically averaged (one site) inversion potential of Clifford et al. [15] in the simulation, we have fitted the potential with the HFD-like model and shown in Fig. 2. The values of the parameters of the fitted potential have been given in Table 1. In order to investigate the difference between the pairpotentials of CF4, we have compared the one-site pair-potentials of Tsuzuki et al. [11], Mahlanen et al. [12], and Goharshadi et al. [13] in Fig. 3. According to this figure, the pair-potential of Mahlanen et al. [12] has smaller mean diameter and deeper potential well (which make the potential more steeper and more attractive) than other pair-potentials. In the other words, we expect that the pair-potential of Mahlanen et al. [12] produces the pressure values smaller than other pair-potentials. 2.2. Three-body potentials Hauschild and Prausnitz (HP) [18] have proposed a three-body correction term in conjunction with the two-body Kihara potential for methane which was a field term proportional to the 9/10 power of the overall density (r0.9) as well as to the attractive-energy

CF4-CF4 10 The inversion pair-pot. of Goharshadi et al. [13] The ab initio pair-pot. of Mahlanen et al. [12] The ab initio pair-pot. of Tsuzuki et al. [11]

1

U2B* 0.1 0.0

-0.5

-1.0

-1.5

-2.0

-2.5 0.8

1.0

1.2

1.4

1.6

1.8

y (r/σ) Fig. 3. Comparison between the different one-site HFD-like pair-potentials.

M. Abbaspour, M. Sheykh / Journal of Fluorine Chemistry 168 (2014) 81–92

86

contribution from the two-body Kihara potential (U2B-ATT): U 3B  ¼ a



r rC

0:9

 U2B -ATT

(9)

where U3B* is the reduced three-body potential, rc is the density at the vapor–liquid critical point, and a is an adjustable constant. The power 9/10 of the density dependence is chosen to obtain better agreement with experimental data at high densities. The total (two-body plus three-body) reduced potential is given by:   UT ¼ U2B  U3B

(10)

The significance of this equation is that it allows us to use twobody (pair) potentials to predict accurately the properties of real fluids without incurring the computational cost of three-body calculations. Eq. (10) obeys the necessary boundary condition that, as r ! 0 the two-body potential is recovered. Hauschild and Prausnitz [18] have adjusted the parameter a by comparison between prediction of liquid densities and experimental data. Their results indicated that a is temperature-dependent. We have determined the a values for fluid CF4 from the inversion HFD-like pair-potential (using the similar method of Hauschild and Prausnitz [18]) at different temperature ranges and fitted to the following equations: 9 8 0:45 T > 1:3T C > > > > > > 3 > >X > > 1 i > = < ðai þ bi T Þ r T C < T  1:3T C > (11) aCF4 ¼ i¼0 > > 3 > > > > > >X 1 i > > ðci þ di T Þ r T  T C > > ; : i¼0

where the density is in (g/cm3) and x is the molar faction of argon or methane. It is shown that the a values (for the mixtures) depend on density and molar fraction at a given temperature. According to previous works [24,29], the following extended total reduced potential using the two-body HFD-like potential (from the inversion method) and the HP many-body potential (Eq. (4)) has been used in our simulation:

    C  C  C10  UT ¼ A exp a y þ b y2  66  88  10 y y y "   #  r 0:9 C6 C8 C10  6  8  10  a (14) rC y y y where y = r/s (s is the distance at which the intermolecular potential has zero value) and rc is the critical density from Ref. [30] and it has been calculated by the following formula for the mixture systems [29]: X rC;mix ¼ mi rC;i (15) i

where mi and rc,i are the mass fraction and critical density of the exponent i in the mixture, respectively. In this work, the parameter a replaced with Eq. (11) and Eqs. ((12) and (13)) for the pure and mixtures, respectively. Wang and Sadus (WS) [19] showed that there is a simple and accurate relationship between the two-body and three-body potential energies of a pure fluid as: U 3B ¼ 

0:85nrU 2B

(16)

es 6

3

where the temperature and density are in (K) and (g/cm ), respectively. The values of the parameters of this equation have been given in Table 3. Eq. (11) shows that at low temperatures the a values not only depend on temperature but also depend on density for a pure substance. The a values for the CF4–Ar and CF4–CH4 interactions have been determined at different densities and molar fractions at constant temperature (T = 373.15 K) and fitted to the following equations:

aCF4 Ar ¼ ð1:1979 þ 1:4280xÞ

   1:6718 þ 5:6485x  5:0987x2 r 

aCF4 CH4 ¼ 6:3450  34:1295x þ 36:1688x2

(12)



   7:9344  40:7580x þ 43:8338x2   2  1  eð69:79þ348:04x352:71x Þr

Table 3 The values of the coefficients of Eq. (11). Parameter

Value

a0 b0 c0 d0 a1 b1 c1 d1 a2 b2 c2 d2 a3 b3 c3 d3

55.9497 0.189714 0.953935 0.00553221 153.124 0.511836 0.595409 0.00344565 137.992 0.455472 0.914857 0.00529576 41.1547 0.133035 3.87997 0.0224892

(13)

where n is the nonadditive coefficient and r is the number density. The total reduced interaction potential, two-body plus three-body (UT ), for a pure fluid is as:

 0:85nr  1 (17) UT ¼ U3B 6

es

We have also applied the expression of Smit et al. [31] for calculating the configuration pressure when the total potential is used: *   + X X 1 dU 2B r i j P con f ¼  ri j 3V dr i j i< j *   +

X X 2=9nr dU 2B r i j 0:85nr2 U 2B  ri j  (18) þ 6 6 dri j es es i< j where the angle brackets represent ensemble averages. To apply Eq. (16) to the mixtures, Wang and Sadus [32] proposed a simple one-fluid approximation: U 3B ¼ 

0:85rU 2B 

e11 s

6 11

x21 n111 þ x1 x2 n112 þ x1 x2 n221 þ x22 n222



(19)

where xi is the molar fraction of component i and nijk is the threebody potential for the three different components i, j, and k. The cross intermolecular parameters can be estimated from the values of pure components as:

n112 ¼ n122 ¼

qffiffiffiffiffiffiffiffiffiffi 3

n21 n2

qffiffiffiffiffiffiffiffiffiffi 3

n1 n22

(20)

(21)

In these equations, 1 stands for CF4 and 2 stands for Ar or CH4. The values of the nonadditivity coefficients for the pure components can be found in Ref. [20]. The total reduced interaction

M. Abbaspour, M. Sheykh / Journal of Fluorine Chemistry 168 (2014) 81–92

potential for a mixture is as follows:  1 UT ¼ U2B

0:85r  2 x n þ x1 x2 n112 þ x1 x2 n221 þ x22 n222 e11 s 611 1 111

!  (22)

We have also applied the expression of Smit et al. [31] for calculating the configuration pressure in the mixtures: Pcon f ¼ 

*   + X X 1 dU 2B r i j ri j 3V dr i j i< j

* X X 2=9r 

þ

i< j

e11 s 611

* 0:85r2 U 2B 



e11 s 611

  +  dU 2B r i j x21 n111 þ x1 x2 n112 þ x1 x2 n221 þ x22 n222  ri j dr i j

x21 n111 þ x1 x2 n112 þ x1 x2 n221 þ x22 n222



+

(23) The significance of the WS three-body equation is that it allows us to use the two-body potentials to predict accurately the properties of real fluids without incurring the computational cost of three-body calculations. Eqs. (17) and (22) have been used with the two-body HFD-like potentials of CF4, CF4–Ar, and CF4–CH4 in this work. Recently, Guzman et al. (GZ) [20] developed a model for an effective Axilrod–Teller–Muto triple-dipole interaction for pure systems based on accurate calculations of three-body effects on the third virial coefficient. The effective three-body interaction is written as a two-body density-dependent potential and is obtained by averaging the function over the position of the third particle.     þ U3B (24) UT ¼ U2B where UT is the reduced total (two-body plus three-body) potential    and U3B is the mean three-body potential and can be approximated by the following function:    e2ð y1Þ U3B ¼ n aðT; sÞr1 z

(25)

where n*(=n/ed9) is the reduced non-additive coefficient and d is the position of the well depth of the potential, r1 = rd3, z(=r/d), and a (T,s) is given by:

aðT; sÞ ¼ a0 ðsÞ þ a1 ðsÞb þ a2 ðsÞb2 þ a3 ðsÞb3 þ a4 ðsÞb4 *

(26)

*

where b = 1/T (T = kT/e) and an(s) vary with s very smoothly and have been fitted by quadratic polynomials in s:

an ðsÞ ¼ an0 þ an1 s þ an2 s2

(27)

87

three-body interactions dominate (the non-additive coefficient of CF4 is greater than Ar and CH4 [20]). It is also shown that the GZ three-body potential [20] has very small effect on the two-body potentials (with the exception of the ab initio potential of Mahlanen et al. [12]). According to Figs. 1 and 2, the effect of adding the HP threebody potential (Eq. (9)) to the HFD-like two-body potentials of the pure and mixture systems is to decrease the minimum of the twobody interactions which increase the effective energy depth (and makes the potential more attractive). In the other words, we expect that the HP three-body term decrease the two-body pressure values. It is also shown that these effects are greater in higher densities (r* = 0.9) in the pure system. But such concentration effect (which is observed for the WS potential) is not for the HP three-body potential for the mixtures. According to Fig. 2, the HP three-body effects are greater at concentrations high in CF4 in the CF4–Ar mixture and are greater at concentrations low in CF4 in the CF4–CH4 mixture. These different behaviors of the modified isotropic density-dependent expression is due to the fact that the HP many-body potential not only accounts for the triple–dipole interaction but also for other effects such as the three-body repulsion or many-body terms. This is mainly due to the adjustable parameter a which is obtained by comparison with the experimental data. In the other words, the HP many-body effects (especially the parameter a ˛ ) depend only on the difference between the two-body results and the experiment. It should be noted that, according to previous works [24,29], the extended HP many-body interaction has been used only with the two-body HFD-like potential from the inversion method. It is also possible to use the HP model with other pair-potentials but the a values should be determined by the same way. In the other words, the different pair-potentials result the different a values. 2.3. Simulation details The MD simulations using MOLDY software [33] have been performed for a system of total 500 molecules in a cubic box and the conventional periodic boundary condition has been applied. The NVT ensemble has been used using a Nose–Hoover thermostat for molecules interacting via the different two-body and total potentials. The size of time steps, Dt, and the number of time steps, nt, and the cutoff radius, rc, have been chosen as 0.001 ps, 500,000 (the equilibration time was 250,000 steps and production runs were 250,000 steps), and 3s, respectively. The long-range correction terms have been evaluated using MOLDY to recover the contribution of the long-range cut-off of the intermolecular potentials on the pressure. 3. Results and discussion

where s is simply a dimensionless measure of the volume of the potential well. The values of the parameters of these equations can be found in Ref. [20]. Eq. (24) has been also used with the two-body HFD-like potential of CF4 in this work. In order to investigate the effect of the different three-body equations to the total potentials, we have presented the different two-body and total potentials of CF4, CF4–Ar, and CF4–CH4 at both low and high densities in Figs. 1 and 2. It is shown that the effect of adding the WS three-body potential (Eqs. (16) and (19)) to the different two-body potentials of the pure and mixture systems, in order to obtain the effective interactions, is to raise the minimum of the two-body interactions which decrease the effective energy depth (and makes the potential more repulsive). In the other words, we expect that the WS three-body term increase the two-body pressure values. It is also shown that these effects are greater in higher densities (r* = 0.9) in the pure system and at concentrations high in CF4 in the mixture systems where the

3.1. Pure CF4 3.1.1. Pressure and EOS The MD simulations have been performed to obtain reduced pressures of fluid carbon tetrafluoride using the different two-body (pair) and total potentials from the HP [18], WS [19], and GZ [20] three-body potentials. The normal conventions have been adopted for the reduced density (r* = rs3), reduced temperature (T* = kT/e), and reduced pressure (P* = Ps3/e) using the HFD-like parameters of CF4 (e/k = 152.1 K, s = 4.7 A˚). Our results of reduced pressure for fluid carbon tetrafluoride in the NVT ensemble have been compared at T* = 1.14 (in the liquid state) and T* = 3 (in the supercritical region) and different densities with the experimental data [30] in Fig. 4. As this figure shows, at low densities, the atom–atom pair-potential of Borodin et al. [10] shows better agreement with the experimental values but at

M. Abbaspour, M. Sheykh / Journal of Fluorine Chemistry 168 (2014) 81–92

88 20

T*= 1.14 16

12

P* 8

4

0 1.02

1.04

1.06

1.08

1.10

1.12

1.14

1.16

1.18

ρ∗ T*= 3

4

3

effective energy depth (and makes the potential more repulsive). But it is shown that the HP many-body potential has decreased the two-body pressure results using the pair-potential of Goharshadi et al. [13] and make them closer to the experiment. This is due to the fact that the effect of adding the HP three-body potential (Eq. (9)) to the HFD-like two-body potential is to decrease the minimum of the two-body interactions (Fig. 1) which increase the effective energy depth (and makes the potential more attractive). It is also expected and shown that the different potentials represent smaller deviations at lower densities and higher temperatures. According to Fig. 4, the total potential from the HP many-body potential shows the best agreement with the experiment. Therefore the MD simulation has been used to determine carbon tetrafluoride EOS using the total potential (Eq. (14)). The reduced pressure values of carbon tetrafluoride have been simulated for the reduced temperatures from 1 to 10 (152.1–1521 K) and reduced density from 0.1 to 1 (140–1400 kg/m3). We have then fitted the simulated results to the EOS of Parsafar et al. [34]. This EOS gives a good description for all types of fluid, including nonpolar, polar, hydrogen-bonded, and metallic, for temperatures ranging from the triple point to the highest temperature for which experimental data are reported [34]. According to the EOS of Parsafar et al. (Z  1)/r*2 in terms of r may be given as: ð Z  1Þ

P*

r2 2

¼eþ

f

r

þ g r2

(28)

where Z = P/rRT and e, f, and g are temperature dependent parameters as:

1

e ¼ e1 þ e2 T  þ

0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

ρ∗ Fig. 4. Comparison between our simulated reduced pressure and the experimental values of fluid CF4 at different reduced temperatures and densities using the different potentials. The solid curve represents the experiment [30], the open circles (*) are for the results using the HFD-like pair-potential of Goharshadi et al. [13], the open squares (&) are for the results using the total potential from the GZ threebody equation, the up open triangles (D) are for the results using total potential from the WS three-body equation, the black circles (*) are for the results using the total potential from the HP three-body equation. The green, blue, and red symbols are for the potentials of Borodin et al. [10], Tsuzuki et al. [11], and Mahlanen et al. [12] too, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

higher densities, the one-site (spherically averaged) pair-potential of Mahlanen et al. [12] produces smaller pressure values and shows better agreement with the experiment than other pairpotentials. This is due to the fact that the pair-potential of Mahlanen et al. [12] has deeper potential well (which make the potential more attractive) than other pair-potentials (Fig. 3). It is also shown that the one-site pair-potential of Goharshadi et al. [13] presents greater pressure values than other pair-potentials which is due to the fact that the pair-potential of Goharshadi et al. [13] has shallower potential well (which make the potential more repulsive) than other pair-potentials (Fig. 3). According to Fig. 4, the two body potentials used for the simulations do not accurately reproduce the two-body interaction between molecules in general, as the potentials used for simulations are too simple and some correction terms need to be added. According to Fig. 4, the both WS and GZ three-body potentials have small effects to the two-body results and increase them to deviate from the experiment. This is due to the fact that the net effect of the WS and GZ three-body potentials is to raise the minimum of the two-body interactions which decrease the

e3 T

f ¼ f 1 þ f 2T  þ g ¼ g1 þ g2T  þ

(29)

f3 T

(30)

g3 T

(31)

where the values of the parameters ei, fi, and gi have been given in Table 4. The temperature dependences of the parameters of the EOS (e, f, and g) have been presented in Fig. 5. According to this figure, these parameters are almost linear versus 1/T*. We have presented the graphs of (Z  1)/r*2 versus r* at different temperatures in Fig. 6. According to this figure, the simulated results are in good agreement with the EOS. It is also shown that these graphs are almost linear but the deviations from linearity become significant for the smaller densities. We have also compared the pressure values calculated from the EOS (Eq. (28)) with those of NIST EOS [30], Mohsen-Nia et al. EOS [35], and Martin and Bhada EOS [36] at two reduced temperatures in Fig. 7. It is shown that our EOS (Eq. (28)) shows better agreement with the NIST EOS than other EOSs. It also should be noted that the EOS (Eq. (28)) covers more extensive temperature and pressure ranges than the NIST EOS. Table 4 The parameters of the EOS (Eqs. (29)–(31)). Parameter

Value

e1 e2 e3 f1 f2 f3 g1 g2 g3

1.1065 0.0621 7.1512 1.9809 0.0559 5.1369 4.4348 0.2643 4.8468

M. Abbaspour, M. Sheykh / Journal of Fluorine Chemistry 168 (2014) 81–92 6

8

g

T* = 2

5

NIST EOS [30] Our EOS (Eq. 28) EOS of Mohsennia et al. [35] EOS of Martin and Bhada [36]

6 4

4

89

3

P* 2

2 1

f

0

0

e

-2

-1 -2

-4

0.0

0.0

0.1

0.2

0.3

0.4

0.5

0.2

0.4

0.6

0.6

ρ

1/Τ∗ Fig. 5. Temperature dependence of the parameters, e, f, and g which are obtained from fitting to the simulation results.

0.8

1.0

1.2

∗

5 *

T=4 3.1.2. Self-diffusion coefficient By using the Einstein formula [37], the self-diffusion coefficient, D, is obtained from the time-dependent mean-square displacement of particles, expected linear at large times: E 1 XD jr i ðtÞ  r i ð0Þj2 t ! 1 6tN i¼1;N

(32)

D ¼ lim

4

3

P* 2

where ri(0) and ri(t) are the positions of the particles at time t = 0 and t, respectively. reduced self-diffusion coefficient, D ¼  The calculated  D ðm=eÞ1=2 =s of carbon tetrafluoride from the time-dependent mean-square displacement using different potentials has been compared at T* = 1.6 at different densities with the experimental values of Khoury and Kobayashi [8] obtained by the pulsed nuclear magnetic resonance method in Fig. 8. According to this figure, the self-diffusion results using the atom–atom pair-potential of Borodin et al. [10] shows better agreement with the experimental values than other pair-potentials. It is also shown that the both WS and GZ three-body potentials have very small effect to the 15

1

0 0.0

0.1

0.2

0.3

0.4

ρ

0.5

0.6

Fig. 7. Comparison between the different EOSs at different reduced temperatures and densities.

self-diffusion and improve the simulated values to give better agreement with the experiment. According to Fig. 8, the HP many-body potential has decreased the two-body pressure results using the pair-potential of

10 4

(Z-1)/ρ*2

5

T*=1.6

0 3

-5

-10

D*

2

-15 -20 1

-25 0.0

0.2

0.4

0.6

0.8

1.0

1.2

ρ∗ Fig. 6. The graphs of (Z  1)/r2 versus density which are obtained from fitting to the simulation results at different temperatures (From bottom to top: T* = 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6, 7, 8, 9, and 10).

0.7

∗

0 0.0

0.2

0.4

ρ∗

0.6

0.8

Fig. 8. Same as Fig. 4 but for the self-diffusion coefficient of CF4.

M. Abbaspour, M. Sheykh / Journal of Fluorine Chemistry 168 (2014) 81–92

90

Goharshadi et al. [13] and make them closest to the experiment. The careful investigation indicated that the a values (used in Eq. (14)) for the self diffusion coefficients of fluid CF4 at constant temperature (T* = 1.6) obey the following equation:

potential (Eq. (22)) to the two-body potentials of the mixture systems is to raise the minimum of the two-body interactions (Fig. 2) which decrease the effective energy depth (and makes the potential more repulsive).

aCF4 ¼ 150:0443  166:5924r þ 392:7079r2  147:9203r3

(33)

5

where r is in (g/cm3). According to Eqs. (11) and (33), the a values are dependent to the physical properties (which are different for the thermodynamic properties and the transport properties).

3

The MD simulation has been performed to obtain reduced pressures of CF4–Ar and CF4–CH4 mixtures using two-body HFDlike potential (Eq. (3)) and total potentials using the HP and WS three-body expressions. The normal conventions have been also adopted for the reduced temperature, density, and pressure using the HFD-like parameters of CF4 ((e/k = 152.1 K, s = 4.7 A˚). In order to compare our simulated results with the experiment, we have used the virial expansion for low and moderate pressure truncated after the third term [38]: P

r T 

Experiment (Eq. 34) Using total pot. from the HP three-body pot. Calculated using the pair-pot. Using total pot. from the WS three-body pot.

4

3.2. CF4–Ar and CF4–CH4 mixtures



T* = 2.45, xAr=0.3

¼ 1 þ B r þ C  r2

P* 2

1

0 0.0

0.1

0.2

0.3

0.4

ρ

(34)

with

0.5

0.6

0.5

0.6

0.7

0.8

∗

5

B ¼ x21 B11 þ 2 x1 x2 B12 þ x22 B22

(35)

    þ 3x21 x2 C112 þ 3x22 x1 C122 þ x32 C222 C  ¼ x31 C111

(36)

xAr=0.5 4

where xi is the molar fraction of component in the mixture and  B11 ð¼ B=s 3 Þ and C111 ð¼ C=s 6 Þ are the reduced second and third virial coefficients, respectively, of pure component i. In these equations, 1 stands for argon or methane and 2 stands for carbon tetrafluoride. The parameters C112 and C122 computed using: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3  2 C  (37) ¼ C111 C112 222  C122 ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3  C 2 C111 222

3

P* 2

1

(38)

0 0.0

 The experimental values of parameters B11 , B12 , B22 , C111 , and *  C222 for CF4–Ar and CF4–CH4 mixtures at T = 2.45 (T = 373.15 K) presented in Table 5 [27]. We have compared our simulated reduced pressure results of CF4–Ar and CF4–CH4 mixtures at constant temperature in the NVT ensemble and three different molar fractions and different densities using the different potentials with the experiment (Eq. (34)) in Figs. 9 and 10, respectively. It is shown that the WS three body potential (Eq. (22)) has increased the two-body results to deviate from the experiment. It is also shown that these effects are greater in higher densities and at concentrations high in CF4 in the mixture systems where the three-body interactions dominate. This is due to the fact that the effect of adding the WS three-body

0.1

0.2

0.3

0.4

0.7

0.8

ρ∗ 6

xAr=0.7 5

4

P*

3

2

Table 5 Reduced virial coefficients of CF4–Ar and CF4–CH4 mixtures at T* = 2.45. Coefficient

CF4–Ara

CF4–CH4a

B11 B12 B22  C111  C222

0.0640 0.2239 0.6894 0.2865 1.0895

0.3359 0.3268 0.6958 0.5710 1.0895

a

Ref. [27].

1

0 0.0

0.1

0.2

0.3

0.4

ρ

0.5

0.6

0.7

0.8

∗

Fig. 9. Comparison between the simulated results of reduced pressure and experimental data of CF4–Ar mixture at different molar fractions and reduced densities using the different potentials.

M. Abbaspour, M. Sheykh / Journal of Fluorine Chemistry 168 (2014) 81–92 6

* T = 2.45, xCH =0.3 5

Experiment (Eq. 34) Using total pot. from the HP three-body pot. Calculated using the pair-pot. Using total pot. from the WS three-body pot.

4

P*

3

2

1

0 0.0

0.1

0.2

0.3

ρ∗

0.4

0.5

0.6

0.7

xCH4=0.5 8

6

P* 4

2

0 0.1

0.2

0.3

ρ

∗

0.4

0.5

0.6

0.7

12

xCH4=0.7

10

8

P*

6

4

2

0 0.0

0.1

0.2

0.3

potentials more attractive). It is also shown that the improvement in agreement is mostly at high densities. This makes sense because three-body effects would be expected to make a greater contribution at higher densities. For non-polar molecules it is expected (not necessarily always) that the triple-dipole three-body potentials (such as the WS potential) increase the two-body results. In the other words, it seems that the three-body potentials can improve the two-body potentials that produce the results a bit below the experiment. But it is shown that the HP potential decreases the two-body results and make them closer to the experiment (Fig. 4, Figs. 9 and 10). Therefore, the HP isotropic density-dependent expression not only accounts for the triple-dipole interaction but also for other effects such as the three-body repulsion or many-body terms. This is mainly due to the adjustable parameter a which is obtained by comparison with the experimental data. 4. Conclusions

10

0.0

91

0.4

0.5

0.6

ρ∗ Fig. 10. Comparison between the simulated results of reduced pressure and experimental data of CF4–CH4 mixture at different molar fractions and reduced.

According to Figs. 9 and 10, the HP many-body potential has decreased the two-body results and make them closer to the experiment. This is due to the fact that the HP three-body potential decreases the minimum of the two-body interactions (Fig. 2) which increases the effective energy depth (and makes the

The MD simulations have been performed to obtain pressures of CF4 and CF4–Ar and CF4–CH4 mixtures using the different (inversion and ab initio) two-body potentials. To take many-body forces into account, the HP [18], WS [19], and GZ [20] expressions have been used with the two-body potentials. The pressure results of pure CF4 indicated that the both WS and GZ three-body potentials have small effects to all of the two-body results and increase them to deviate from the experiment which is due to the raising the minimum of the two-body interactions by the both three-body potentials which make the potential more repulsive. The HP many-body potential has decreased the pressure results (using the HFD-like pair-potential) and make them closer to the experiment which is due to the decreasing the minimum of the two-body interactions by the three-body potential (which make the potential more attractive). The MD simulation using the total potential from the HP many-body potential has been also used to determine an equation of state which may be used as a reference for fluid carbon tetrafluoride. The obtained pressure results of the CF4–Ar and CF4–CH4 mixtures indicated that the WS three body potential has increased the two-body results to deviate from the experiment. It is also shown that these effects are greater in higher densities and at concentrations high in CF4 in the mixture systems where the three-body interactions dominate. This is due to the fact that the effect of adding the WS three-body potential (Eq. (22)) to the two-body potentials of the mixture systems is to raise the minimum of the two-body interactions (which makes the potential more repulsive). We have also extended the HP many-body potential as a function of density, temperature, and molar fraction and used with the two-body HFD-like potentials of CF4–Ar and CF4–CH4 mixtures to improve the prediction of the pressure values (which decreases the two-body results by decreasing the minimum of the two-body potentials) without requiring an expensive three-body calculation. We have also concluded that the HP isotropic density-dependent expression not only accounts for the triple-dipole interaction but also for the other effects such as the three-body repulsion or many-body terms. We have also compared the simulated self-diffusion coefficient of fluid carbon tetrafluoride using the different two-body and total potentials with the experiment. The results indicated that the atom–atom pair-potential of Borodin et al. [10] shows better agreement with the experimental values than other pair-potentials. It is also shown that the both WS and GZ three-body potentials have very small effect to the self-diffusion but the HP many-body potential has decreased the two-body pressure results using the HFD-like pair-potential and make them closest to the

92

M. Abbaspour, M. Sheykh / Journal of Fluorine Chemistry 168 (2014) 81–92

experiment. We also concluded that the a values are different for the thermodynamic properties and the transport properties. References [1] I. Waldner, A. Bassen, H. Bertagnoli, K. Todheide, G. Strauss, A.K. Soper, J. Chem. Phys. 107 (1997) 10667–10674. [2] K. Williams, K. Gupta, M. Wasilik, J. Microelectromech. Syst. 12 (2003) 761–777. [3] J. Vrabec, J. Stoll, H. Hasse, J. Phys. Chem. B 105 (2001) 12126–12133. [4] K.E. MacCormack, W.G. Schneider, J. Chem. Phys. 19 (1951) 849–855. [5] T.H. Spurling, A.G. De Roco, T.S. Stornick, J. Chem. Phys. 48 (1968) 1006–1008. [6] C.A. Gough, S.E. Debolt, P.A. Kollman, J. Comput. Chem. 13 (1992) 963–970. [7] B.J. Palmer, J.L. Anchell, J. Phys. Chem. 99 (1995) 12239–12248. [8] F. Khoury, R. Kobayashi, J. Chem. Phys. 55 (1971) 2439–2445. [9] R.D. Parra, X.C. Zeng, J. Mol. Struct. (Theochem.) 503 (2000) 213–220. [10] O. Borodin, G.D. Smith, D. Bedrov, J. Phys. Chem. B 106 (2002) 9912–9922. [11] S. Tsuzuki, T. Uchimaru, M. Mikami, S. Urata, J. Chem. Phys. 116 (2002) 3309–3315. [12] R. Mahlanen, J. Jalkanen, T.A. Pakkanen, Chem. Phys. 313 (2005) 271–277. [13] E.K. Goharshadi, M. Abbaspour, A. Morsali, Ind. Eng. Chem. Res. 42 (2003) 2256–2261. [14] G. Vayner, Y. Alexeev, J. Wang, T.L. Windus, W.L. Hase, J. Phys. Chem. A 110 (2006) 3174–3178. [15] A.A. Clifford, E. Dickinson, G.P. Matthews, E.B. Smith, J. Chem. Soc., Faraday Trans. 172 (1976) 2917–2922. [16] G. Marcelli, B.D. Todd, R.J. Sadus, J. Chem. Phys. 115 (2001) 9410–9413. [17] M.T. Oakley, H. Do, R.J. Wheatley, Fluid Phase Equilib. 290 (2010) 48–54.

[18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38]

T. Hauschild, J.M. Prausnitz, Mol. Simul. 11 (1993) 177–185. L. Wang, R.J. Sadus, J. Chem. Phys. 125 (2006) 144509–144513. O. Guzman, F. del Rio, J.E. Ramos, Mol. Phys. 109 (2011) 955–967. L. Zarkova, P. Pirgov, Vacuum 48 (1997) 21–27. E.K. Goharshadi, M. Abbaspour, Fluid Phase Equilib. 212 (2003) 53–65. E.K. Goharshadi, M. Jami-Alahmadi, B. Najafi, Can. J. Chem. 81 (2003) 866–871. S. Salemi, M. Abbaspour, M. Ghabdian, J. Supercrit. Fluids 89 (2014) 119–127. G.C. Maitland, M. Rigby, E.B. Smith, W.A. Wakeham, Intermolecular Forces: Their Origin and Determination, Clarendon Press, Oxford, 1986. P. Clancy, D.W. Cough, G.P. Matthews, E.B. Smith, G.C. Maitland, Mol. Phys. 30 (1975) 1397–1407. J.H. Dymond, E.B. Smith, The Second Virial Coefficients of Pure Gases and Mixtures, Oxford University, Oxford, 1980. S. Bock, E. Bich, E. Vogel, Chem. Phys. 257 (2000) 147–156. M. Abbaspour, E. Nameni, J. Supercrit. Fluids 74 (2013) 61–69. National Institute of Standard and Technologie (NIST). Avalable at http://webbook.nist.gov/chemistry/fluid/. B. Smit, T. Hauschild, M. Prausnitz, Mol. Phys. 77 (1992) 1021–1031. L. Wang, R.J. Sadus, J. Chem. Phys. 125 (2006) 74503–74508. The MOLDY program was coded by K. Refson (2001), and can be downloaded from the internet at hhttp://www.ccp5.ac.uk/moldy/moldy.htmli. G.A. Parsafar, H.V. Spohr, G.N. Patey, J. Phys. Chem. B 113 (2009) 11977. M. Mohsen-Nia, H. Modarress, G.A. Mansoori, Fluid Phase Equilib. 206 (2003) 27–39. J.J. Martin, R.K. Bhada, AlChE J. 17 (1971) 683–688. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. E.A. Mason, T.H. Spurling, The Virial Equation of State, Pergamon, London, 1969.