J. Phys. Chem. B 2008, 112, 13885–13897
13885
Surface Tensions of Linear and Branched Alkanes from Monte Carlo Simulations Using the Anisotropic United Atom Model F. Biscay,†,‡ A. Ghoufi,‡ F. Goujon,† V. Lachet,‡ and P. Malfreyt*,† Laboratoire de Thermodynamique et Interactions Mole´culaires, FRE 3099 CNRS, UniVersite´ Blaise Pascal, 63177 Aubiere Cedex, and IFP, 1-4 AV. de Bois Pre´au, 92852 Rueil Malmaison Cedex, France ReceiVed: July 11, 2008; ReVised Manuscript ReceiVed: September 2, 2008
The anisotropic united atoms (AUA4) model has been used for linear and branched alkanes to predict the surface tension as a function of temperature by Monte Carlo simulations. Simulations are carried out for n-alkanes (n-C5,n-C6,n-C7, and n-C10) and for two branched C7 isomers (2,3-dimethylpentane and 2,4dimethylpentane). Different operational expressions of the surface tension using both the thermodynamic and the mechanical definitions have been applied. The simulated surface tensions with the AUA4 model are found to be consistent within both definitions and in good agreement with experiments. 1. Introduction A good knowledge of interfacial properties is of great use to optimize many industrial processes. For instance, in the chemical industry, such properties are useful to design efficient product formulations. In the oil and gas industry, interfacial forces play an important role, especially in oil recovery schemes or in drilling fluid designs. In oil recovery processes, interfacial behaviors strongly influence saturation, distribution, and multiphase flow of fluids in porous media. In such operations, up to three fluid phases can be involved as follows: oil, water, and gas (water being present either at its irreducible saturation or in bulk after secondary recovery by waterflooding and gas being possibly formed during the first depletion or being injected during enhanced recovery processes). To model such oil recovery schemes and to predict their efficiency, accurate knowledge of the three-phase fluid distribution and the threephase flow characteristics is required, meaning, among others, accurate knowledge of oil-gas, oil-water, and water-gas interfacial tensions (γog, γow, and γwg). These three interfacial tensions are, for instance, required to estimate the spreading coefficient S ) γwg - γog - γow that governs the ability of the oil to form spreading films on water in the presence of gas.1 Interfacial tensions are strongly dependent on thermodynamic conditions such as temperature, pressure, and fluid compositions and are often delicate to assess experimentally due to the high pressure conditions that exist in oil and gas reservoirs or in wells. The oil and gas industry and, more generally, the chemical industry therefore need reliable models to predict these interfacial properties. Some modeling methods are classically used to complement available experimental measurements in the case of simple fluids. Most of them are correlative methods using the parachor theory.2,3 However, this approach has a limited predictive capacity because it depends on the empirical relationships that are used. In recent years, the square-gradient (SG) theory4-10 has been used with the statistical associating fluid theory (SAFT) to calculate the interfacial properties of some pure fluids and mixtures.11-16 Whereas the surface tension of pure fluids is reasonably predicted with the SAFT-SG approach, * To whom correspondence should
[email protected]. † Universite ´ Blaise Pascal. ‡ IFP.
be
addressed.
E-mail:
the predictive capabilities are significantly reduced in the case of mixtures due to the need of empirical adjustable parameters. Recently, Gloor et al. established an accurate Helmholtz free energy density functional theory (DFT) based upon the SAFT theory for attractive potentials of variable range (SAFT-VR). The resulting SAFT-VR DFT free energy density functional has been successfully applied for the calculation of the surface tension of linear alkanes.17 By contrast, the temperature dependence of the surface tension is much less reproduced in the case of associating compounds; one solution consists of using the surface tension in developing the intermolecular potential model. For a more comprehensive review of the SAFT approach and its many extensions in the case of liquid-vapor interfaces, the reader is directed to recent papers17,18 and references therein. These modern molecular theories can provide useful interfacial properties but cannot describe the interface at the molecular level. The ability to associate an accurate calculation of the surface tension with a molecular description of the interface for systems involving different types of interactions is one of the most important outcomes of advances in molecular simulation techniques over the past couple of decades. As a result, we start to use molecular simulation, and more specifically Monte Carlo (MC) simulation, to assess the interfacial study of different systems encountered in the petroleum industry and to get more fundamental insight of these interfacial phenomena. In this work, preliminary studies involving pure hydrocarbon compounds (linear and branched alkanes) are proposed using MC simulations for liquid-vapor surface tension calculations. To allow the application of these simulation methods to such industrial problems, it is essential that a consistent intermolecular potential is available for a large range of compounds and thermophysical properties. During the past few years, a lot of efforts have been invested in the development of the anisotropic united atoms (AUA) intermolecular potential proposed by Toxvaerd.19,20 This potential, in which hydrogen atoms are not explicitly considered, is more realistic and flexible than classical united atoms potentials21-24 and less demanding in computing time than the all atoms models classically used for the simulation of liquids.25,26 Thanks to an original optimization method using vapor-liquid phase equilibrium data (coexisting densities, vapor pressures, and vaporization enthalpies), the AUA potential has been
10.1021/jp806127j CCC: $40.75 2008 American Chemical Society Published on Web 10/11/2008
13886 J. Phys. Chem. B, Vol. 112, No. 44, 2008 extended with good success to assess thermodynamic bulk properties of a variety of n-alkanes,27 branched alkanes,28 and cyclic alkanes.29 The set of reference molecules that have been used to perform this optimization is composed of three linear alkanes (ethane, n-pentane, and n-dodecane) for the optimization of CH2 and CH3 group parameters, one branched alkane (isobutane) for the optimization of the CH group in branched alkanes, and two cyclic alkanes (cyclopentane and cyclohexane) for the optimization of the saturated cyclic CH2 group parameters.30 No interfacial tension value has been included in this optimization procedure. After little readjustment, these AUA models have also lead to very satisfying results for the calculation of transport properties of n-alkanes, such as shear viscosity and self-diffusion coefficient.31 To the best of our knowledge, no attempt to apply these AUA models for interfacial properties calculations has yet been proposed. As a result, we propose here to calculate the surface tension of linear and branched alkanes with the AUA model. The calculation of the surface tension by molecular simulation can be performed using different methodologies. The most natural method consists of modeling two coexisting phases in physical contact in a liquid-slab arrangement. This methodology takes advantage of providing a molecular description of the interfacial region. However, it depends on the set up conditions (truncation procedures,32-34 cutoff radius,32,35 longrange corrections,33,35 (LRCs) and size of the simulation cell35,36). Other methodologies avoid the need to establish an interface. These approaches37-42 are based upon the thermodynamic definition of the surface tension and express the surface tension from the contribution to the free energy due to the formation of the interface. Some of these methods have been successfully applied on Lennard-Jones38,40,43 (LJ), square well,41 n-alkanes,42 and LJ mixtures.44 Throughout this paper, we have adopted the route that consists of explicitly modeling the liquid-vapor interface by using a liquid-slab arrangement. This leads to heterogeneous systems where a nonuniformity of the density number takes place along the direction normal to the interface. These two phase simulations require one to pay special attention to a certain number of factors that can impact the results of surface tension: the size effects due to the number of particles,36,45 the range of interaction,32,35,46 the truncation effects involved in the calculation of the potential and corresponding force,32-34 the method (mechanical or thermodynamic routes) used for the surface tension calculation,37,47,48 and the LRCs to the surface tension.33-35,48,49 Concerning the calculation of the surface tension, a certain number of operational expressions have been established. The most general working expression uses macroscopic normal and tangential pressures that can be related to the derivative of the pair potential. The final form has been obtained by Kirkwood and Buff50-53 (KB) and gives a macroscopic scalar surface tension. In the case of a planar liquid-vapor interface lying in the x, y plane, the density gradient takes place in the z direction normal to the surface. The surface tension can be expressed as ∞ ∫ -∞ [pN(z) - pT(z)] dz, where pN(z) and pT(z) are local values of the normal and tangential components of the pressure tensor, respectively. There are many ways of expressing the local components of the pressure, which depend on the contour joining of two interacting molecules. Irving and Kirkwood6,54-56 (IK) use a straight line to join the two particules. This choice is the most natural and the one generally made. Nevertheless, the scalar value of surface tension is invariant to the choice of the pressure tensor. Until recently, only the method of IK was
Biscay et al. TABLE 1: LJ Parameters of the Interaction Sites for the Linear and Branched Alkanes CH CH2 CH3
σ (Å)
ε/k (K)
δ (Å)
3.363 3.461 3.607
50.98 86.29 120.15
0.646 0.384 0.216
designed to provide a profile of the surface tension. For this reason, we established the local version (KBZ)48 associated with the KB expression. The KB, KBZ, and IK methods are based upon the mechanical definition of the surface tension. The testarea47 (TA) method based upon the perturbation formalism was proposed by Gloor and co-workers for the calculation of the surface tension. This method takes advantage of expressing the surface tension as a difference of energy between a reference state and a perturbed state characterized by an infinitesimal increase or decrease of the surface. However, the aspect of the local surface tension was missing within the TA method. This led us to establish the local expression of the surface tension calculated from the TA approach.34 The paper is organized as follows. In the next two sections, we describe the potential model, the computational procedures, and the different definitions of the surface tensions with their corresponding operational expressions. In section 4, we discuss the impact of the size effects on the surface tension calculation, and we present the results of surface tensions of some linear and branched alkanes. We conclude in section 5. 2. Simulation Methodology 2.1. Potential Model. Following the choice made by Ungerer and co-workers,27-29 linear and branched alkanes are described using the AUA4 model where CH3, CH2, and CH groups are represented with a single center of force located near the geometric center of the atoms of each group. The AUA model consists of a displacement of the LJ centers of force toward the hydrogen atoms. The center of force i of a CH2 group is located on the external bissector of the Ci-1-Ci-Ci+1 angle and the CH3 force center i is located on the Ci-Ci-1 axis. In the case of a CH group28 linked with three carbon atoms C1, C2, and f C3, the force center is placed on the vector b T ) C1CH + f
f
C2CH + C3CH. The magnitude of the shift between the carbon center and the interaction site is represented by a third adjustable parameter δ (see Table 1). The total configurational energy U is defined by
U ) UINTRA + UINTER + ULRC
(1)
where UINTRA, UINTER, and ULRC are the intramolecular, intermolecular, and LRCs energy contributions. The intramolecular interactions include contributions from the bending energy, torsion energy, and nonbonded interactions. No stretching energy was included as the bond lengths are considered to be constant and equal to 1.535 Å (see Table 2). The expressions for the bending potential energy
1 Ubend ⁄ k ) kbend(cos θ - cos θ0)2 2
(2)
and for the torsion potential energy 8
Utors ⁄ k )
∑ ai(cos φ)i
(3)
i)0
are identical to those of Toxvaerd,20 where k is the Boltzmann constant, θ is the bond angle, and φ is the dihedral angle. The
MC Simulations Using the AUA Model
J. Phys. Chem. B, Vol. 112, No. 44, 2008 13887
TABLE 2: Parameters for Molecular Weight, Carbon-Carbon Distance, Bending Potential, and Torsion Potential MW (g/mol)
CH CH2 CH3
C-C distance (Å) bending
θ0 (deg) kbend (K) θ0 (deg) kbend (K) a0 (K) a1 (K) a2 (K) a3 (K) a4 (K) a5 (K) a6 (K) a7 (K) a8 (K) a0 (K) a1 (K) a2 (K) a3 (K) a4 (K) a5 (K) a6 (K) a7 (K) a8 (K) a0 (K) a1 (K) a2 (K) a3 (K) a4 (K) a5 (K) a6 (K) a7 (K) a8 (K)
C-CH-C C-CH2-C C-CH2-CH2-C
torsion
13.03 14.03 15.03 1.535 112 72700 114 74900 1001.35 2129.52 -303.06 -3612.27 2226.71 1965.93 -4489.34 -1736.22 2817.37 373.05 919.04 268.15 -1737.21 0 0 0 0 0 373.05 919.04 268.15 -1737.21 0 0 0 0 0
C-CH-CH-C
C-CH-CH2-C
parameters of the bending and torsion potentials are given in Table 2. The intramolecular nonbonded interactions occur between groups separated by more than three bonds and were computed using the same AUA4 LJ 6-12 potential as intermolecular interactions (see just below). The intermolecular interactions due to the repulsion-dispersion interactions were computed using the truncated LJ potential Ni
N-1 N
UINTER )
Nj
∑ ∑ ∑ ∑ uLJ(riajb) i)1 j>i a)1 b)1 Ni
N-1 N
)
Nj
[( ) ( ) ]
∑ ∑ ∑ ∑ 4εab i)1 j>i a)1 b)1
σab riajb
12
-
σab riajb
Figure 1. Direct (γTA,D), reverse (|γTA,I|), and average (γTA) surface tension values calculated from the TA method for the liquid-vapor interface of n-pentane at T ) 300 K as a function of ξ. The dotted line corresponds to the average value of γ calculated from ξ ) 5 × 10-4.
TABLE 3: Number of Alkane Molecules, Box Dimensions, and Cutoff Value system
N
n-pentane 500 n-hexane 300 n-heptane 300 n-decane 250 2,3-dimethylpentane (2,3-DMP) 300 2,4-dimethylpentane (2,4-DMP) 300
Lx (Å) Ly (Å) Lz (Å) rc (Å) 35.0 34.7 34.6 34.6 35.5 33.1
35.0 34.7 34.6 34.6 35.5 33.1
235.0 150.0 225.0 250.0 200.0 200.0
12.0 12.0 12.0 12.0 12.0 12.0
TABLE 4: Relative Probabilities of the Different Types of MC Moves linear alkanes (%)
branched alkanes (%)
45 35 20 0
35 35 20 10
translation rotation regrowth by configurational bias pivot
ULRC was then calculated by summing up all of the local contributions of each slab. The ULRC term was then added in the total energy of the system to be used in the Metropolis scheme. The LRCs to the total energy within each kth slab are defined by two parts49
6
Nz
(4)
where riajb is the distance between atom a in molecule i and atom b in molecule j, εab is the energy parameter of the interaction, and σab is the LJ core diameter. Ni is the number of atoms in the molecule i. The LJ parameters for the interactions between unlike sites were calculated using the Lorentz-Berthelot combining rules
ULRC )
∑
Nz
uLRC(zk) )
i)1
(1) (2) (zk) + uLRC (zk)] ∑ [uLRC
where Nz is the number of slabs and (1) uLRC (zk) )
(6)
i)1
[ ( ) ( )]
Ni Nj σ12 σ6ab ab 8π εab 1/3 9 - 3 F(zk)2Vs 3 rc rc a)1 b)1
∑∑
(7)
and
εab ) (εaaεbb)
1⁄2
1 σab ) (σaa + σbb) 2
(5)
The list of the LJ parameters for the AUA4 model is given in Table 1. As the geometry of the system shows a heterogeneity along the axis normal to the interface (z-axis), we calculated the LRC to the repulsion-dispersion energy as a function of zk by splitting the cell into slabs of width δz. The total LRC energy
∫
∞ (2) uLRC (zk) ) πF(zk)Vs r c
dr
Ns
∫-r d∆z∑ [F(zk+i) r
i)1
F(zk+i-1)]rULJ,m(r) (8) F(zk) and Vs are, respectively, the density and the volume of the slab k. ∆z is defined as the difference z - zk. Ns is the number of slabs between z and zk. rc is the cutoff radius (see Table 3),
13888 J. Phys. Chem. B, Vol. 112, No. 44, 2008
Biscay et al.
TABLE 5: Surface Tension Values (mN m-1) of n-Pentane Calculated at T ) 300 K and T ) 360 K from MC Simulations as a Function of the Interfacial Areaa γKB
γIK
Lx ) Ly
Lz
γLRC
γTOT
γLRC
35.0 37.8 40.9 43.8 46.8
235.0 235.0 235.0 235.0 235.0
5.01 4.81 4.81 3.81 4.71
16.519 17.019 17.018 15.627 16.519
4.31 4.51 4.61 4.51 4.61
35.0 39.2 42.4 45.4 47.7
235.0 235.0 235.0 235.0 235.0
2.91 3.11 3.11 3.11 3.11
9.410 9.916 8.414 9.616 9.519
2.71 3.01 2.91 3.11 3.01
γTA γTOT
γKBZ
γLRC
γTOT
γLRC
γTOT
γEXP
T ) 300 K 15.521 16.719 18.318 16.327 16.519
4.12 4.22 4.31 4.33 4.33
16.120 17.023 17.119 16.629 16.720
3.81 3.71 3.81 3.81 3.91
14.821 16.019 16.118 15.627 15.719
15.3 15.3 15.3 15.3 15.3
T ) 360 K 9.210 9.716 8.314 9.616 9.418
2.01 2.12 2.31 2.32 2.42
9.411 9.817 8.515 9.718 9.720
2.21 2.51 2.51 2.61 2.51
8.710 9.216 7.814 9.116 8.918
9.0 9.0 9.0 9.0 9.0
a The LRCs and the total surface tension are given for each method. The subscripts give the accuracy of the last decimal(s), that is 14.821 means 14.8 ( 2.1. The experimental surface tensions71 (γEXP) are reported for comparison.
TABLE 6: Surface Tension Values (mN m-1) of n-Pentane Calculated at T ) 300 K and T ) 360 K from MC Simulations as a Function of the Lz Dimensiona γKB
γIK
Lx ) Ly
Lz
γLRC
γTOT
γLRC
35.0 35.0 35.0 35.0 35.0
235.0 251.0 266.0 282.0 297.0
5.01 4.91 4.31 4.01 4.51
16.519 16.818 16.417 15.524 16.823
4.31 4.61 4.51 4.61 4.81
35.0 35.0 35.0 35.0 35.0
235.0 258.0 275.0 293.0 311.0
2.91 3.01 3.21 3.11 3.21
9.410 9.514 9.114 9.918 10.120
2.71 2.11 3.41 3.21 3.01
γTA γTOT
γKBZ
γLRC
γTOT
γLRC
γTOT
γEXP
T ) 300 K 15.521 16.617 17.218 15.824 17.117
4.12 4.42 4.21 4.63 4.72
16.120 17.017 17.019 16.929 17.121
3.81 3.91 4.21 4.11 4.01
14.821 16.219 16.518 15.927 15.919
15.3 15.3 15.3 15.3 15.3
T ) 360 K 9.210 8.713 9.316 9.417 9.818
2.01 2.52 2.42 2.11 2.22
9.411 9.714 8.015 8.413 8.717
2.21 2.71 2.81 2.91 2.71
8.710 9.215 8.718 9.215 9.520
9.0 9.0 9.0 9.0 9.0
a The LRCs and the total surface tension are given for each method. The subscripts give the accuracy of the last decimal(s), that is, 16.219 means 16.2 ( 1.9. The experimental surface tensions71 (γEXP) are reported for comparison.
ULJ,m(r) is the intermolecular energy with r being the distance between the two centers of mass of molecules i and j. Ni
ULJ,m(r) )
Nj
[( ) ( ) ]
∑ ∑ 4εab a
b
σab r
12
-
σab r
6
(9)
The first part of the long-range contribution has an analytical form identical to the one associated with a homogeneous system but uses the local density F(zk) of the slab. The second part consists of a double integral that contains a series of density differences, which render this part cumbersome to calculate. As concerns the second part of the LRCs to the total configurational energy, it has been shown 35,49 that it represents only a minor contribution to the total long-range energy. 2.2. Computational Procedures. The simulation box was a rectangular parallelepipedic box of dimensions LxLyLz (Lx ) Ly) with N alkane molecules. The details of the geometry of the system are given in Table 3. The periodic boundary conditions were applied in the three directions. MC simulations were performed in the NVT ensemble. Each cycle consisted of N randomly selected moves with fixed probabilities: translation of the center of mass of a random molecule, rotation of a randomly selected molecule around its center of mass, and change of the internal conformation by using the configurational bias regrowth move and the pivot move.28,57 The pivot consisting of rotating a part of the molecule located beyond some atom around this atom was only used for branched alkanes. The
frequency of each type of move depends on the alkane considered (see Table 4). The initial configuration was built by placing N molecules on nodes of a fcc-centered orthorombic lattice included in a cubic box. The number N of molecules and the intial volume were chosen to build a dense phase. The nodes of the lattice as well as the orientation of the molecules were randomly chosen. MC simulations in the NpT ensemble were first performed on this bulk monophasic fluid configuration. The dimension of the resulting box was increased along the z-axis by placing two empty cells on both sides of the bulk liquid box. A typical MC run consisted of 200000 cycles for equilibration and 200000 cycles for the production phase. The thermodynamic properties were calculated every 10 cycles leading to the storage of 20000 configurations. The statistical errors for these properties were estimated using the Jackknife method.58 This method used four superblocks, which are formed by combining three blocks of 5000 configurations. One of the specificity of our MC methodology was the use of the LRCs to the configurational energy in the Metropolis scheme. The total LRC energy ULRC was updated after each move of molecular position and was added in the total energy of the system to be used in the Metropolis scheme. 3. Surface Tension In what follows, we present the different operational expressions of the surface tension with the corresponding
MC Simulations Using the AUA Model
J. Phys. Chem. B, Vol. 112, No. 44, 2008 13889
Figure 2. Surface tension values of the liquid-vapor interface of n-pentane at T ) 300 K and T ) 360 K calculated from different methods as indicated in the legend as a function of (a) the interfacial area and (b) the Lz dimension. pN(zk) - pT(zk) profiles calculated from the IK method in (c) two different interfacial areas and (d) different box dimensions along the z-axis. pN,LRC(zk) - pT,LRC(zk) profiles of the LRCs to the surface tension in (e) two interfacial areas and (f) different Lz dimensions. In parts c-f, the integrals of the pN - pT profiles are plotted in the right axis.
expressions of their LRCs. Let us consider a system of N molecules with two planar liquid-vapor surfaces lying in the x,y plane where z represents then the direction normal to the surface. 3.1. KB Relation. The molecular surface tension γKB was first introduced by KB50 and makes use of the molecular virial expression to give the following relationships:
1 γKB ) 2A
〈∑ ∑ ∑ ∑ N-1 i)1
Ni
Nj
rij · riajb - 3zij · ziajb dU(riajb) 2riajb driajb j)i+1 a)1 b)1 N
〉
(10) Specific LRCs were developed by Blockhuis et al.59 and are based upon the approximations that the radial distribu-
tion function is equal to unity for r greater than the cutoff radius and that the density profile can be fitted to a hyperbolic function. This LRC expression takes the following form:
π γKB,LRC ) (Fl - Fv)2 2
∫01 ds∫r+∞ dr coth( 2rs d ) c
dULJ,m 4 r × dr (3s3 - s) (11)
where d is an estimation of the thickness of the interface and s is a parameter defined by s ) (zi - zj)/rij. 3.2. IK Definition. The method of IK54 expresses the surface tension from the local components of the pressure tensor
13890 J. Phys. Chem. B, Vol. 112, No. 44, 2008
γIK )
1 2
∫-LL ⁄2⁄2 [pN(zk) - pT(zk)]dz z
Biscay et al.
(12)
z
where pN(zk) and pT(zk) are the normal and tangential components of the pressure tensor along the normal to the surface, respectively. The method of IK54 is based upon the notion of the force across a unit area. The pressure tensor is then written as a sum of a kinetic term and a potential term resulting from the intermolecular forces. Whereas the first term is well-defined, the potential term is subjected to arbitrariness because there is no unique way to determine which intermolecular forces contribute to the stress across dA. There are many ways of choosing the contour joining two interacting particles. IK54 have chosen as a contour the straight line between the two particles. Other choices are possible and result from the lack of uniqueness in the definition of the microscopic stress tensor. The components of the pressure tensor6,55,56 in the IK definition are expressed by
pRβ(zk) ) 〈F(zk)〉kBT I +
1 A
〈∑ ∑ N-1 N
(rij)R(Fij)β
i)1 j>i
( ) ( )〉
1 zk - zi × θ |zij| zij θ
zj - zk zij
(13)
where I is the unit tensor and T is the input temperature. R and β represent x, y, or z directions. θ(x) is the unit step function defined by θ(x) ) 0 when x < 0 and θ(x) ) 1 when x g 0. A is the surface area normal to the z-axis. The distance zij between two molecular centers of mass is divided into Ns slabs of thickness δz. Following IK, the molecules i and j give a local contribution to the pressure tensor in a given slab if the line joining the centers of mass of molecules i and j crosses, starts, or finishes in the slab. Each slab has 1/Ns of the total contribution from the i-j interaction. The normal component pN(zk) is equal to pzz(zk), whereas the tangential component is given by 1/2[pxx(zk) + pyy(zk)]. Fij in eq 13 is the intermolecular force between molecules i and j and is expressed as the sum of all of the site-site forces acting between these two molecules. Ni
Fij )
π (2) pT,LRC (zk) ) - F(zk) 2
∫r∞ dr∫-rr d∆z[F(z) - F(zk)] × c
dULJ,m(r) 2 [r - (∆z)2] (16) dr The first term of pN,LRC(zk) and pT,LRC(zk) is identical to that used in homogeneous molecular simulations by using a local density F(zk) instead of a scalar density F, whereas the second term takes into account the density differences between the slabs. From these LRC expressions, it is then possible to calculate the LRC parts relative to the surface tension. The operational expression of the LRC part of the surface tension49 within the IK formalism is then given by
Vs π γIK,LRC(zk) ) F(zk) 2 A
∫r ∫ ∞
dr
c
Ns
r
d∆z -r
∑ i)1
dULJ,m 2 [F(zk+i) - F(zk+i-1)] [r - 3(∆z)2] (17) dr The total LRC to the surface tension is obtained by summing up all of the contributions to the local values of each bin and dividing the result by 2. 3.3. TA Method. The TA method47 is based upon a thermodynamic route and expresses the surface tension as a change in the free energy for an infinitesimal change in the surface area in the constant NVT ensemble. This infinitesimal change in the area is performed throughout a perturbation process for which the perturbed system (stateA+∆A) is obtained from an infinitesimal change ∆A of the area A of the reference system. The box dimensions [L(A+∆A) ,L(A+∆A) ,Lz(A+∆A)] in the x y perturbed systems are changed using the following transforma-
Nj
∑ ∑ (fiajb)
a)1 b)1
Nj
Ni
) -
r
∑ ∑ riajb iajb
a)1 b)1
dU(riajb) driajb
(14)
The appropriate LRC to the IK definition of the normal and tangential components of the pressure tensor has been derived by Guo and Lu49 and is composed of two parts as expressed in eqs 15 and 16.
TABLE 7: Estimated Critical Temperature and Critical Density of the Hydrocarbons Considereda
(1) (2) pN,LRC(zk) ) pN,LRC (zk) + pN,LRC (zk)
)-
2π 2 F (zk) 3
πF(zk) F(zk)]
∫r∞ drr3 c
dULJ,m(r) dr
r
c
dULJ,m(r) (∆z)2 dr
Fc (kg m-3)
Tc (K)
∫r dr∫-r d∆z[F(z) ∞
Figure 3. Vapor-coexistence curve of n-alkanes calculated from MC simulations (open circles). The continuous lines represent the experimental coexisting densities.70 The experimental critical points are represented by filled symbols.
(15)
As concerns the tangential pressure, the first term is identical to the first term of the normal component, whereas the second term is expressed by
n-pentane n-hexane n-heptane n-decane 2,3-DMP 2,4-DMP a
AUA
exp.
%
AUA
exp.
%
461 503 533 610 530 509
470 508 540 618 537 520
-1.9 -1.0 -1.2 -1.2 -1.3 -2.1
247 235 234 234 250 234
232 233 233 228 254 239
6 -0.8 +0.4 +2.5 +1.5 +2.0
Computed values (AUA) are compared with experimental data.70,71
MC Simulations Using the AUA Model
J. Phys. Chem. B, Vol. 112, No. 44, 2008 13891
(A+∆A) (A+∆A) tions: LA+∆A ) L(A) ) L(A) ) x x √1+ξ,Ly y √1+ξ,Lz (A) Lz /(1 + ξ), where ξ,1. The area (A+∆A) of the perturbed (A) (A) (A) state thus equals to L(A) x Ly (1 + ξ) and ∆A is equal to Lx Ly ξ. These transformations conserve the volume of the box in the perturbed state. This means that it is possible to express the surface tension as a ratio of partition function between the perturbed and the reference states. U(A)(rN) and U(A+∆A)(r′N) are the configurational energies of the systems with an area A and a configurational space rN and an area A+∆A and a configurational space r′N, respectively. The transformations used in the perturbation process lead to the following equality dr′N ) drN. This expression is the key for writing the surface tension as the average of exp(-∆U/kT) in the constant NVT ensemble
γTA )
∂F ( ∂A )
N,V,T
) lim -
〈 (
ξf0
ln exp )
kBT × ∆A
[U(A+∆A)(r′N) - U(A)(rN)] kBT
)〉
0
∑ lim ξf0 k
〈 (
kBT [U(A+∆A)(zk, r′N) - U(A)(zk, rN)] ln exp ∆A kBT
)〉
k,A
(18) 〈...〉k,A indicates that the average is carried out over the reference state and the k slabs. U(A+∆A)(zk,r′N) and U(A)(zk,rN) are the configurational energies of the slab k in the perturbed and reference states. The ambiguity6 focuses on the part of the interaction energy to be included in the volume Vs of the slab. We adopt the definition of Ladd and Woodcok60 and choose to assign in the slab centered on zk two energy contributions: one contribution due to the energy between the molecules within the slab and a second contribution due to the energy of the molecules within the slab with those outside the slab. The energy of the slab at the position zk is defined as N
Uzk )
N
Ni
Nj
j*i
a
b
∑ ∑ ∑ ∑ Hk(zi)ULJ(riajb)
1 2 i)1
(19)
where Hk(zi) is a top-hat function with functional values of
{
Hk(zi) ) 1 for zk - δz < zi < zk + δz 2 2 Hk(zi) ) 0 otherwise
(20)
Using this definition, we respect the following condition
∫V dzkUz ) U k
(21)
where U is the total configurational energy of the simulation box and V is its volume. As concerns the LRC contributions to the surface tension calculated from the TA method, the total LRC contribution is expressed as a function of the total LRC energy ULRC (see eq 22). We have already shown34 that the difference between the (zk') and u(1),(A) u(1),(A+∆A) LRC LRC (zk) vanishes because the transformation conserves the volume. Thus, the total value of the tail correction of the surface tension within the TA formalism is given by
〈 (
(A+∆A) (A) kBT [ULRC - ULRC ] γTA,LRC ) lim ln exp ξf0 ∆Aξ kBT
)〉
0
Ns
)
∑ γLRC(zk)
(22)
k
The calculation of the surface tension was carried out in the direct (γTA,D) and reverse (γTA,I) directions. The calculation of the direct direction involves an increase of the surface area (A + ∆A), whereas a decrease of the surface area (A - ∆A) is performed in the reverse path. The ensemble average was carried out over the configurations of the reference system, whereas the configurations of the perturbed system did not participate to the Markov chain of states. Thermodynamic consistency requires that the surface tension in the direct and reverse directions must be equal in magnitude and in opposite sign. This is satisfied when the configuration space of the perturbed system matches the one of the reference system. This requirement implies use of an appropriate value of ξ. The value of ξ must satisfy two constraints:47 this value should be small enough to allow an accurate calculation of the surface tension from eq 18 and large enough to provide reasonable statistics for the Boltzmann factor. The surface tension value was averaged over the two directions as (γTA,D - γTA,I)/2. However, when the results differ between the two directions, there is no fundamental basis to justify that the errors cancel when averaging the direct and inverse values. One way to check the TA calculation is to compare the average surface tension with other values resulting from different approaches and to check the impact of the value of ξ on the calculated results. The influence of the value of ξ is illustrated in Figure 1 where the values of the surface tension of the n-pentane at 300 K are reported in the direct and reverse paths for ξ ranging from 10 -6 to 10 -2. We observe in Figure 1 that the value of the surface tension is little dependent on ξ between 10-5 and 10-3 with similar values in the backward and forward directions. The value of ξ ) 5 × 10-4 used both in this work and in most of the previous studies,34,45,47 represented by a dotted line in Figure 1, is shown to be a reasonable value for the calculation of the surface tension using the TA approach. 3.4. Local Expression of the Surface Tension from the Virial Route (KBZ). The original working expression of KB50 does not provide a profile of the surface tension as a function of the direction normal to the interface. The profile becomes a key element from a methodological viewpoint to check the validity of the calculation concerning the stabilization of the interfaces, the independence between the two interfaces, and the constancy of the pseudolocal surface tension γ(z) in the bulk regions. In a previous work, we established a local version48 of the surface tension based upon the KB expression. The working expression was obtained from the derivative of the potential with respect to the surface and is referred as the KBZ method in this paper.
γKBZ )
〈 〉
∂U ) ∂A 0
Ns
∑ k
〈 〉 ∂Uzk ∂A
0
(23)
Ns
)
∑ γKBZ(zk) k
After further elaboration, γKBZ(zk) can be expressed as
(24)
13892 J. Phys. Chem. B, Vol. 112, No. 44, 2008
γKBZ(zk) ) -
Biscay et al.
48ε
1 ij A(riajb, rij) × ∑ ∑ ∑ ∑ riajb riajb i∈k
a
j*i
b
[( ) ( ) ] σij riajb
12
-
1 σij 2 riajb
6
where the function A(riajb,rij) is defined as
[
A(riajb,rij) ) (riajb)x
(rij)x 4Lx2
+ (riajb)y
(rij)y 4Ly2
- (riajb)z
(rij)z 2LyLx
(25)
] (26)
The working local LRC expression was also established in this previous work48 and is then given by the following equation
γKBZ,LRC(zk) ) πF(zk)
(
Vs 2A
∫r∞ dr∫-rr d∆z[F(z) - F(zk)] × c
)(
)
∂ULJ,m(r) r - 3(∆z)2 ULJ,m(r) + r (27) r ∂r 2
4. Results and Discussion 4.1. Size Effects. In a preliminary step, we performed a number of system-size checks for the calculation of the surface tension. These checks are important because we focused our attention on the computational efficiency; therefore, we tried to use values of cutoff identical to those usually used in the Gibbs ensemble Monte Carlo (GEMC) simulations.61-63 Additionally, it was established in a previous paper35 that the dependence of the surface tension with the cutoff radius was canceled as appropriate LRCs to the surface tension were taken into account. In consequence, the value of the cutoff radius was
kept constant for all of the system-size checks. We calculated the surface tension for systems with N ranging from 500 to 900 molecules. The impact of the number of molecules is studied through both an increase of the interfacial area by keeping Lz constant and an increase of the Lz dimension by keeping the surface area constant. The values of the surface tension of the n-pentane are reported in Tables 5 and 6 and are shown in parts a and b of Figure 2. We do not observe a trend for an evolution of the surface tension over this range of system sizes. The changes of the surface tension with respect to the system size are within the fluctuations estimated by using the variation of the block averages. Such observations were already made on the liquid-vapor interface of carbon dioxide48,64,65 and water.66 However, some recent works36,45 report an oscillatory function of the surface tension of LJ fluids for small surface areas. Let us underline that the surface areas used in this work correspond to areas from which the oscillatory effect is significantly attenuated. To further investigate the system size dependence, the difference between the normal and the tangential pressure profiles for N ) 500 and 800 molecules as well as the surface tension profiles are reported in Figure 2c within the IK definition. The pN - pT profiles show two symmetric independent interfaces, and the γ(z) profile exhibits no contribution to the surface tension from the liquid and vapor bulk phases. The shape of these profiles is not dependent on the system size within the statistical fluctuations. Figure 2d presents the same profiles calculated with different box dimensions along the z-axis. As expected, the profiles present a liquid region that extends further when the number of molecules is increasing. However, the magnitude of the peaks in the (pN - pT) profiles is not dependent on the number of molecules leading to surface tension values
Figure 4. Surface tensions of the liquid-vapor interface of n-alkanes as a function of temperature calculated from different methods as indicated in the legend.
MC Simulations Using the AUA Model
J. Phys. Chem. B, Vol. 112, No. 44, 2008 13893
TABLE 8: Surface Tension Values (mN m-1) of n-Alkanes Calculated from MC Simulations Using Different Operational Expressionsa γKB T (K) γLRC
γIK
γTOT
γLRC
γTA
γTOT
γLRC
γKBZ
γTOT
γLRC
γTOT
γEXP
300 330 360 390 420
5.01 16.519 3.12 13.529 2.91 9.410 1.61 5.520 1.01 3.110
n-pentane 4.31 15.521 4.12 16.120 3.61 12.618 2.74 12.920 2.71 9.210 2.01 9.411 1.71 5.221 1.52 5.622 1.11 2.910 0.82 2.911
3.81 15.521 15.3 3.01 12.018 12.1 2.21 8.710 9.0 1.41 4.920 6.0 0.91 2.710 3.4
360 390 420 450 480
4.11 13.320 2.41 8.913 1.71 6.88 0.91 3.411 0.51 1.511
n-hexane 3.71 12.819 2.63 12.022 2.71 8.213 2.12 8.514 1.91 6.58 1.52 5.79 1.01 2.511 0.92 2.711 0.61 1.111 0.51 1.212
3.11 12.219 11.5 2.21 7.813 8.6 1.51 6.28 5.9 0.81 2.310 3.5 0.51 1.011 1.4
400 420 440 460 480
3.51 10.912 2.02 9.020 1.91 7.018 1.41 5.014 1.01 3.815
n-heptane 2.61 10.012 1.93 10.414 2.31 8.418 1.74 8.822 2.01 6.019 1.85 6.520 1.41 4.615 1.32 5.018 1.11 2.913 1.01 3.116
2.11 1.91 1.71 1.31 0.91
4.91 19.413 4.62 15.219 3.61 10.318 2.41 6.018 1.11 2.516
n-decane 18.511 4.84 19.014 14.318 4.44 14.821 9.217 3.55 9.920 5.117 2.24 5.619 2.117 0.92 2.117
350 400 450 500 550
4.61 4.31 3.21 2.01 0.91
9.612 10.1 8.118 8.3 5.717 6.6 4.413 5.0 2.714 3.5
4.61 18.412 18.7 4.31 14.220 14.4 3.31 9.319 10.4 1.91 4.818 6.6 0.91 2.014 3.3
a
The LRCs contribution and the total surface tension are given for each method. The subscripts give the accuracy of the last decimal(s), that is, 15.521 means 15.5 ( 2.1. The experimental71 surface tensions (γEXP) are reported for comparison.
TABLE 9: Calculated Surface Tensions (mN m-1) from MC Simulations Obtained by Using the AUA-427 and UA-TraPPE22 Models for the Three Linear Alkanes γKB γLRC
γTOT
UA AUA-4
5.41 5.01
UA AUA-4 UA AUA-4
γIK γLRC
γTOT
γTA γLRC
γKBZ
γTOT
γLRC
γTOT
16.810 16.518
n-pentane, T ) 300 K 5.01 16.411 4.73 16.614 4.33 15.521 4.12 16.120
4.21 3.83
15.611 15.521
4.32 4.11
12.022 13.320
n-hexane, T ) 360 K 3.91 11.621 3.24 12.124 3.71 12.819 2.63 13.022
3.31 3.11
10.921 12.219
3.91 3.51
10.323 10.912
n-heptane, T ) 400 K 4.61 11.023 3.52 11.124 2.61 10.012 1.93 10.414
3.81 2.11
10.222 9.612
in close agreement within the standard deviations. Parts e and f of Figures 2 present the pN - pT profiles of the LRCs to the pressure calculated from eqs 15 and 16 as well as the integral of the corresponding profile shown in the right axis. When the interfacial area is changed (Figure 2e), the change in the LRC to the surface tension is about 4% when N increases from 500 to 800. At fixed interfacial area, we observe in Figure 2f that the final value of γLRC is also little dependent on the system size. The fluctuations in γLRC are within 10% of the value of γ calculated with the smallest system size. This is the first time that the different profiles of γ and of its LRCs are given as a function of the system size. These profiles are very powerful to check that the two-phase simulations lead to two wellestablished interfaces that contribute in the same way to the surface tension. This preliminary study showed that the smallest system size (N ) 500) allowed to recover the main features of
a well-equilibrated two-phase system within a reasonable computational effort. 4.2. Linear Alkanes. The coexisting densities in the bulk phases are obtained by fitting the simulated density profiles to a hyperbolic tangent function of the form:
[
2(z - z0) 1 1 F(z) ) (Fl + Fv) - (Fl - Fv) tanh 2 2 d
]
(28)
In eq 28, Fl and Fv are the coexisting densities of the liquid and vapor phases, z0 is the position of the Gibbs dividing surface, and d is an approximate measure of the thickness of the interface. The density profile can also be fit using the error function model 67,68
[
2(z - z0) 1 1 F(z) ) (Fl + Fv) - (Fl - Fv)erf 2 2 √2d
]
(29)
where erf(z) ) (2/π)∫z0 exp (-t2)dt. The critical point69 (Tc and Fc) was calculated by fitting the simulated coexistence densities to the law of rectilinear diameters
Fl + Fv ) Fc + A(T - Tc) 2
(30)
and to the scaling law for the densities using an Ising type critical exponent β ) 0.325
Fl - Fv ) B(T - Tc)β
(31)
Figure 3 shows the coexisting densities of some linear alkanes calculated from the two-phase simulations using the AUA4 potential and the tanh fits. Fitting the coexisting densities to the erf function gives differences smaller than 0.1% for the liquid densities and 2% for the vapor densities. These changes in the coexisting densities are very small as compared to those induced by the statistical fluctuations calculated from the variations of the block averages.33,35 A comparison between the erf and the tanh fits will be shown in the case of the density profiles of branched alkanes. The comparison with experiments70,71 indicates that good predictions are obtained. The critical points of these n-alkanes are given in Table 7 and are shown in Figure 3. The results show that the error in the critical temperature does not exceed 10 K, leading to a relative maximum deviation of 2%. The critical densities are shown to be in reasonable agreement with experiments even if this property is subject to a larger uncertainty with a maximum deviation of 6% for n-C5. However, the critical densities are consistent with those calculated from both GEMC simulations using the AUA4 potential27 and the direct MC simulations using the UA TraPPE model.35 Figure 4 shows the results of liquid-vapor surface tensions obtained from different definitions and operational expressions for the four studied n-alkanes from n-C5 to n-C10. First, we observe that the surface tensions are correctly predicted using the AUA4 potential within the statistical errors. Table 8 reports the LRC contributions and the intrinsic part with respect to the method used. We observe that the surface tension calculated from the sum of the intrinsic part and LRCs is independent of the definition used once appropriate LRCs were added.34,48,49 We observe that the LRCs calculated using the TA approach are smaller than those calculated from the IK and KB methods. However, this difference is compensated by a larger intrinsic part of the surface tension with the TA approach. This was explained in a recent paper34 by the fact that the TA formalism uses the potential in its operational expression, whereas IK and KB methods used the derivative of the potential leading to slight
13894 J. Phys. Chem. B, Vol. 112, No. 44, 2008
Figure 5. (a) Vapor-coexistence curve of two heptane isomers (2,3DMP and 2,4-DMP) calculated from MC simulations (open circles). The continuous lines represent the experimental coexisting densities.71 The experimental critical points are represented in filled symbols. (b) Molecular density profiles along the z-axis of the liquid-vapor equilibrium at T ) 420 K. (c) Molecular density profile of the 2,3DMP at 420 K calculated from the MC simulation and from the tanh and erf fits.
differences in the surface tension due to the truncation procedures. We also observe from Table 8 that the LRCs to the surface tension represent up to 30% of the total surface tension with a cutoff value of 12 Å. This underlines the importance of the tail corrections and the need to associate an adequate expression of γLRC consistent with that used for the intrinsic γ part. Under these conditions, the surface tension is independent of the method used and of the value of cutoff. It is then established that the AUA4 force field is capable of reproducing the surface tensions of n-alkanes even if the parameters of this force field have not been optimized using values of surface tensions. This is the same observation that we made for the calculation of the surface tension of n-alkanes with the UA TraPPE models.33,35 Table 9 reports surface
Biscay et al.
Figure 6. (a) Surface tensions of the liquid-vapor interface of 2,3dimetylpentane and 2,4-DMP as a function of temperature calculated from different methods (same symbols as the ones used as in Figure 4). (b) Surface tensions of the two isomers calculated from the IK method as a function of temperature.
tensions of n-C5, n-C6, and n-C7 with respect to the UA TraPPE and AUA4 models. The first conclusion that can be drawn from the values of Table 9 is that the two models lead to similar surface tensions within the statistical uncertainties. It is also interesting to note that the LRCs to the surface tensions calculated with AUA4 are smaller than those calculated with UA TraPPE; thus, the intrinsic part of the surface tension is slightly increased using the AUA4 potential. However, the total surface tension is not affected by the model. This trend is observed for the different definitions IK, TA, KB, and KBZ used in this work. 4.3. Branched Alkanes. We focus now on the two-phase simulations of the 2,3-DMP and 2,4-DMP. These two hydrocarbon isomers differ only by the position of one methyl group. The liquid-vapor coexistence curve of these iso-alkanes is given in Figure 5a and shows that the vapor and liquid densities are
MC Simulations Using the AUA Model
J. Phys. Chem. B, Vol. 112, No. 44, 2008 13895
Figure 7. Configuration of the 2,4-DMP liquid-vapor interface at T ) 340 K.
TABLE 10: Surface Tension Values (mN m-1) of Branched Alkanes Calculated from MC Simulations Using Different Operational Expressionsa γKB T (K) γLRC
γTOT
γIK γLRC
γTOT
γTA γLRC
γTOT
γKBZ γLRC
γTOT
γEXP
300 340 380 420 460
5.71 18.526 3.61 15.912 2.61 12.021 1.61 8.919 0.91 5.017
2,3-DMP 4.51 19.325 3.74 18.329 4.01 15.611 2.91 16.014 3.11 11.719 2.13 11.922 1.81 8.018 1.73 7.319 1.11 5.38 0.92 5.49
3.71 16.625 19.3 3.41 15.011 15.5 2.61 11.219 11.8 1.51 7.718 8.3 0.91 5.18 5.0
300 340 380 420 460
3.71 17.634 3.21 14.329 2.41 10.623 1.51 7.129 0.51 4.915
2,4-DMP 4.61 16.433 3.13 16.835 3.61 13.528 2.84 14.231 2.71 9.722 1.93 10.024 1.71 6.027 1.35 6.332 0.61 3.814 0.53 4.017
3.81 15.633 17.5 3.01 13.928 13.7 2.21 9.222 10.2 1.41 5.727 6.8 0.51 3.714 3.7
a The LRCs contribution and the total surface tension are given for each method. The subscripts give the accuracy of the last decimal(s), that is 15.521 means 15.5 ( 2.1. The experimental surface tensions (γEXP) are reported for comparison.71
well-predicted using the direct MC simulations. The critical properties are given in Table 7. We also observe that the simulation is able to reproduce the differences in density between the two isomers over the range of temperatures studied. This was already shown from GEMC simulations72 using the same potential. The critical region is also well-described. The difference in the critical properties of the two molecules is also well-predicted with a maximum average deviation with experimental values of 2%. Figure 5b shows the density profiles along the direction normal to the surface for the two isomers. This profile calculated at T ) 420 K is shown to reproduce clearly the differences in the liquid densities between the two isomers. Part c of Figure 5 shows the density profile of the 2,3-DMP calculated from the MC simulations and the fits using either the tanh or the erf functions. As observed in this figure, the functional form (tanh or erf) chosen to fit the density profile does not affect the calculated coexisting densities within the statistical uncertainties. The ability to capture some small differences in the density of the two isomers is of prime importance if precise predictions of surface tensions values are expected. The surface tensions of the 2,3-DMP and 2,4-DMP are presented in Figure 6a,b for the KB, TA, IK and KBZ methods. The values of the total surface tensions and of the associated LRCs are listed in Table 10. Parts a and b of Figure 6 show that the surface tensions of these iso-alkanes calculated using the AUA4 model are found to be in good agreement with the DIPPR correlations.71 Additionally, Figure 6c plots the surface tension calculated using the IK method as a function of temperature for the 2,3-DMP and 2,4-DMP. Interestingly, we check that the calculation of the surface tension with the AUA4 potential allows us to reproduce the experimental differences of surface tension between the two isomers. This is an excellent result, which shows the performance of the AUA4 potential for
Figure 8. Surface tension profiles of the n-pentane at T ) 300 K calculated for (a) the intrinsic part and (b) the LRCs using the IK, TA, and KBZ approaches.
the prediction of the surface tension. A typical configuration of the liquid-vapor interface of the 2,4-DMP at 340 K is given in Figure 7. 4.4. Discussion. The calculation of the surface tension can be carried out using different definitions. As a function of the definition used, either the potential or the derivative of the potential is involved in the calculation. Additionally, some methods provide a profile of the surface tension as a function of the direction normal to the surface. Because the calculation of the surface tension is associated with relatively large statistical fluctuations, it is interesting to compare the different methods between them. We propose here to show the different profiles of the surface tension calculated from the IK, TA, and KBZ methods. The profiles of the intrinsic part and of the LRCs are given in Figure 8a,b, respectively. The values calculated from the profiles using the TA and KBZ approaches are similar to those calculated from the original expressions of TA and KB. Figure 8a confirms that the different IK, TA, and KBZ profiles match very well within the errors bars. These profiles also exhibit two symmetric and independent interfaces with no contribution due to the bulk liquid and vapor phases. The profiles
13896 J. Phys. Chem. B, Vol. 112, No. 44, 2008 of the LRCs to the surface tension are also shown in Figure 8b and establish the stability of the two interfaces. 5. Conclusions A certain number of works have contributed to show that the shift of force center used in AUA potentials allows significant improvement of the potential transferability without an increase of the computational time. Because the parameters of the AUA4 potential were not optimized on surface tensions, the calculation of this property using this force field represents a challenging test of transferability for this potential. To our knowledge, this was the first time that the calculation of surface tensions was performed on alkanes using the AUA4 potential. Because of some relatively important statistical fluctuations of the surface tension, different operational expressions have been applied to calculate the liquid-vapor interfacial properties of linear and branched alkanes. The fundamental checks of the system sizes have been carried as a function of the interfacial area and of the length of the box dimension along the normal to the surface. The AUA4 potential has allowed efficient predictions of the surface tension of n-alkanes. The surface tensions calculated from AUA4 are shown to be similar to those calculated using the UA TraPPE model within the standard deviations. The AUA4 model produces larger intrinsic surface tension values and smaller LRC contributions with respect to the UA TraPPE model. The agreement with experimental values is also found to be excellent in the case of branched alkanes with the 2,3DMP and 2,4-DMP. Additionally, our calculations of the surface tension allow the prediction of subtle experimental differences between the two isomers. The coexisting densities calculated from the direct MC twophase simulations are consistent with those calculated from the GEMC method using the same potential. The critical temperature and densities are also well-reproduced. This work establishes the AUA4 potential as an efficient model for the calculation of surface tension of linear and branched alkanes. Acknowledgment. We acknowledge the support of the French Agence Nationale de la Recherche (ANR), under Grant SUSHI (ANR-07-BLAN-0268) from “SimUlation de Syste`mes He´te´roge`nes et d’Interfaces”. References and Notes (1) Adamson, A. W. Phys. Chem. Surf; Interscience Publishers: New York, 1960. (2) Guggenheim, E. A. J. Chem. Phys. 1945, 13, 253. (3) Zuo, Y.-X.; Stenby, E. H. Can. J. Chem. Eng. 1997, 75, 1130. (4) van der Waals, J. D. Z. Phys. Chem. 1894, 13, 657. (5) Rowlinson, J. S. J. Stat. Phys. 1979, 20, 197. (6) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Clarendon Press: Oxford, 1982. (7) Cahn, J. W.; Hilliard, J. E. J. Chem. Phys. 1958, 28, 258. (8) Carey, B. S.; Scriven, L. E.; Davies, H. T. AIChE J. 1980, 26, 705. (9) Abbas, S.; Nordholm, S. J. Colloid Interface Sci. 1994, 166, 481. (10) Cornelisse, P. M. W.; Peters, C. J.; de Swaan Arons, J. J. Chem. Phys. 1997, 106, 9820. (11) Kahl, H.; Enders, S. Fluid Phase Equilib. 2000, 172, 27. (12) Kahl, H.; Enders, S. Phys. Chem. Chem. Phys. 2002, 4, 931. (13) Duque, D.; Pa`mies, J. C.; Vega, L. F. J. Chem. Phys. 2004, 121, 11395. (14) McCabe, C.; Kiselev, S. B. Ind. Eng. Chem. Res. 2004, 43, 2839. (15) Mejı´a, A.; Pa`mies, J. C.; Duque, D.; Segura, H.; Vega, L. F. J. Chem. Phys. 2005, 123, 034505. (16) Mejı´a, A.; Segura, H.; Wisniak, J.; Polishuk, I. Phys. Chem. Liq. 2006, 44, 45. (17) Gloor, G. J.; Jackson, G.; Blas, F. J.; del Rio, E. M.; de Miguel, E. J. Phys. Chem. C 2007, 111, 15513.
Biscay et al. (18) Gloor, G. J.; Jackson, G.; Blas, F. J.; del Rio, E. M.; de Miguel, E. J. Chem. Phys. 2004, 121, 12740. (19) Toxvaerd, S. J. Chem. Phys. 1990, 93, 4290. (20) Toxvaerd, S. J. Chem. Phys. 1997, 107, 5197. (21) Smit, B.; Karaborni, S.; Siepmann, I. J. J. Chem. Phys. 1995, 102, 2126. (22) Martin, M. G.; Siepmann, I. J. J. Phys. Chem. B 1998, 102, 2569. (23) Martin, M. G.; Siepmann, I. J. J. Phys. Chem. B 1999, 103, 4508. (24) Nath, S. K.; Escobedo, F. A.; de Pablo, J. J. J. Chem. Phys. 1998, 108, 9905. (25) Jorgensen, W. L.; Madura, J. D. J. Am. Chem. Soc. 1984, 106, 6638. (26) Jorgensen, W.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225. (27) Ungerer, P.; Beauvais, C.; Delhommelle, J.; Boutin, A.; Rousseau, B.; Fuchs, A. J. Chem. Phys. 2000, 112, 5499. (28) Bourasseau, E.; Ungerer, P.; Boutin, A.; Fuchs, A. H. Mol. Simul. 2002, 28, 317. (29) Bourasseau, E.; Ungerer, P.; Boutin, A. J. Phys. Chem. B 2002, 106, 5483. (30) Ungerer, P.; Tavitian, B.; Boutin, A. Applications of Molecular Simulation in the Oil and Gas Industry; Editions Technip: Paris, France, 2005. (31) Nieto-Draghi, C.; Ungerer, P.; Rousseau, B. J. Chem. Phys. 2006, 125, 044517. (32) Trokhymchuk, A.; Alejandre, J. J. Chem. Phys. 1999, 111, 8510. (33) Goujon, F.; Malfreyt, P.; Simon, J.; Boutin, A.; Rousseau, B.; Fuchs, A. H. J. Chem. Phys. 2004, 121, 12559. (34) Ibergay, C.; Ghoufi, A.; Goujon, F.; Ungerer, P.; Boutin, A.; Rousseau, B.; Malfreyt, P. Phys. ReV. E 2007, 75, 051602. (35) Goujon, F.; Malfreyt, P.; Boutin, A.; Fuchs, A. J. Chem. Phys. 2002, 116, 8106. (36) Orea, P.; Lopez-Lemus, J.; Alejandre, J. J. Chem. Phys. 2005, 123, 114702. (37) Binder, K. Phys. ReV. A 1982, 25, 1699. (38) Errington, J. R. Phys. ReV. E 2003, 67, 012102. (39) Errington, J. R. J. Chem. Phys. 2003, 118, 9915. (40) Potoff, J. J.; Panagiotopoulos, A. Z. J. Chem. Phys. 2000, 112, 6411. (41) Singh, J. K.; Kofke, D. A.; Errington, J. R. J. Chem. Phys. 2003, 119, 3405. (42) Singh, J. K.; Errington, J. R. J. Phys. Chem. B 2006, 110, 1369. (43) Shen, V. K.; Mountain, R. D.; Errington, J. R. J. Phys. Chem. B 2007, 111, 6198. (44) Shen, V. K.; Errington, J. R. J. Chem. Phys. 2006, 124, 024721. (45) Errignton, J. R.; Kofke, D. A. J. Chem. Phys. 2007, 127, 174709. (46) Lopez-Lemus, J.; Alejandre, J. Mol. Phys. 2002, 100, 2983. (47) Gloor, G. J.; Jackson, G.; Blas, F. J.; de Miguel, E. J. Chem. Phys. 2005, 123, 134703. (48) Ghoufi, A.; Goujon, F.; Lachet, V.; Malfreyt, P. Phys. ReV. E 2008, 77, 031601. (49) Guo, M.; Lu, B. C. Y. J. Chem. Phys. 1997, 106, 3688. (50) Kirkwood, J. G.; Buff, F. P. J. Chem. Phys. 1949, 17, 338. (51) Buff, F. B. Z. Elektrochem. 1952, 56, 311. (52) McLellan, A. G. Proc. R. Soc. A 1952, 213, 274. (53) McLellan, A. G. Proc. R. Soc. A 1953, 217, 92. (54) Irving, J. H.; Kirkwood, J. G. J. Chem. Phys. 1950, 18, 817. (55) Walton, J. P. R. B.; Tildesley, D. J.; Rowlinson, J. S.; Henderson, J. R. Mol. Phys. 1983, 48, 1357. (56) Walton, J. P. R. B.; Tildesley, D. J.; Rowlinson, J. S. Mol. Phys. 1986, 58, 1013. (57) Dodd, L. R.; Boone, T. D.; Theodorou, D. N. Mol. Phys. 1993, 78, 961. (58) Efron, B. The Jackknife the Boostrap and Other Resampling Plans; Design Institute for Physical Properties, SIAM: Philadelphia, 1982. (59) Blokhuis, E. M.; Bedeaux, D.; Holcomb, C. D.; Zollweg, J. A. Mol. Phys. 1995, 85, 665. (60) Ladd, A. J. C.; Woodcok, L. V. Mol. Phys. 1978, 36, 611. (61) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813. (62) Panagiotopoulos, A. Z.; Quirke, N.; Stapleton, M.; Tildesley, D. J. Mol. Phys. 1988, 63, 527. (63) Panagiotopoulos, A. Z. Mol. Simul. 1992, 9, 1. (64) Ghoufi, A.; Goujon, F.; Lachet, V.; Malfreyt, P. J. Chem. Phys. 2008, 128, 154716. (65) Ghoufi, A.; Goujon, F.; Lachet, V.; Malfreyt, P. J. Chem. Phys. 2008, 128, 154718. (66) Alejandre, J.; Tildesley, D. J.; Chapela, G. A. J. Chem. Phys. 1995, 102, 4574. (67) Sides, S. W.; Grest, G. S.; Lacasse, M. D. Phys. ReV. E 1999, 60, 6708. (68) Mitrinovic, D. M.; Tikhonov, A. M.; Li, M.; Huang, Z.; Schlossman, M. L. Phys. ReV. Lett. 2000, 85, 582. (69) Rowlinson, J.; Swinton, F. Liquid and Liquid Mixtures, 3rd ed.; Butterworths: London, 1982.
MC Simulations Using the AUA Model (70) Gmehling, J. CODATA Bull. 1985, 58, 56. (71) Rowley, R. L.; Wilding, W. V.; Oscarson, J. L.; Tundel, N. A.; Marshall, T. L.; Daubert, T. E.; Danner, R. P. DIPPR Data Compilation of Pure Compound Properties; Design Institute for Physical Properties, AIChE: New York, 2002.
J. Phys. Chem. B, Vol. 112, No. 44, 2008 13897 (72) Ahunbay, M. G.; Kranias, S.; Lachet, V.; Ungerer, P. Fluid Phase Equilib. 2005, 228, 311 .
JP806127J
