178
J. Phys. Chem. B 2009, 113, 178–187
All-Atom Force Field for the Prediction of Vapor-Liquid Equilibria and Interfacial Properties of HFA134a Robson P. S. Peguin, Ganesh Kamath, Jeffrey J. Potoff, and Sandro R. P. da Rocha* Department of Chemical Engineering and Materials Science, Wayne State UniVersity, 5050 Anthony Wayne DriVe, Detroit, Michigan 48202 ReceiVed: July 14, 2008; ReVised Manuscript ReceiVed: October 23, 2008
A new all-atom force field capable of accurately predicting the bulk and interfacial properties of 1,1,1,2tetrafluoroethane (HFA134a) is reported. Parameterization of several force fields with different initial charge configurations from ab initio calculations was performed using the histogram reweighting method and Monte Carlo simulations in the grand canonical ensemble. The 12-6 Lennard-Jones well depth and diameter for the different HFA134a models were determined by fitting the simulation results to pure-component vaporequilibrium data. Initial screening of the force fields was achieved by comparing the calculated and experimental bulk properties. The surface tension of pure HFA134a served as an additional screening property to help discriminate an optimum model. The proposed model reproduces the experimental saturated liquid and vapor densities, and the vapor pressure for HFA134a within average errors of 0.7%, 4.4%, and 3.1%, respectively. Critical density, temperature, vapor pressure, normal boiling point, and heat of vaporization at 298 K are also in good agreement with experimental data with errors of 0.2%, 0.1%, 6.2%, 0%, 2.2%, respectively. The calculated surface tension is found to be within the experimental range of 7.7-8.1 mN · m-1. The dipole moment of the different models was found to significantly affect the prediction of the vapor pressure and surface tension. The ability of the HFA134a models in predicting the interfacial tension against water is also discussed. The results presented here are relevant in the development of technologies where the more environmentally friendly HFA134a is utilized as a substitute to the ozone depleting chlorofluorocarbon propellants. 1. Introduction Chlorofluorocarbons (CFCs) were widely used in the past due to their nontoxic and nonflammable nature and their high thermodynamic efficiency.1 However, with the banning of ozone-depleting substancessas dictated by the Montreal ProtocolsCFCs have been almost completely phased out.2 The hydrofluoroalkanes (HFAs) have emerged as environmentally friendly alternatives to CFCs. While HFAs are greenhouse gases, they have zero ozone-depleting potential.3 HFAs have been shown to be fairly effective as lubricants, water repellant and stain repellant products, blowing agents in the production of insulating foam, chemical reagents, refrigerants, and industrial solvents.4 Of particular interest to us is the use of HFAs as propellants in pressurized metered-dose inhaler (pMDI) formulations.5 HFA134a (1,1,1,2-tetrafluoroethane) and HFA227 (1,1,1,2,3,3,3heptafluoropropane) are nontoxic,6 noncarcinogenic,7 and are rapidly absorbed and eliminated from the body.8 pMDIs are the most widely used oral inhalation therapy devices for treating pulmonary disorders, including bronchial asthma and obstructive pulmonary disease.9 pMDIs have also been recognized as candidate devices for the delivery of therapeutic biomolecules, including peptides, oligonucleotides, and proteins to and through the lungs.10 However, there have been innumerous reformulation issues associated with the replacement of CFCs with HFAs in pMDI formulations. While novel HFA-based pMDIs have been recently reported in the literature,11-16 more extensive development in the area has been limited by challenges in studying * To whom correspondence should be addressed. Tel.: +1-313-577-4669. Fax: +1-313-578-5820. E-mail:
[email protected].
compressible systems in situ (under pressure),17,18 and by the lack of fundamental understanding of the solvent environment in HFAs.11,19,20,17 Water-in-HFA microemulsions17,20-24 have also been suggested as potential carriers for the noninvasive delivery of therapeutic molecules to and through the lungs.13,25 An important first step in the design of aqueous reverse aggregates in HFAs is the evaluation of the properties of the bare HFA|water (HFA|W) interface,20,26 which to a large extent dictates the activity of the amphiphiles adsorbed at the interface.27 However, there are several limitations in assessing the microscopic properties of interfacial systems experimentally. Investigations are even more challenging in the case of systems under pressure.27,28 On the other hand, molecular modeling tools can provide insight not only into the bulk thermodynamic propertiess which can be directly checked against the experimental valuessbut also into the microstructure of the interface. Such computational studies may, therefore, serve as a guide for the design of interfacially active species capable of forming and stabilizing reverse aqueous aggregates in HFAs. A reliable description of the thermodynamic and microstructural properties of interfaces by atomistic simulations can only be achieved, however, if accurate molecular models are available. There have been several previous studies aimed at developing potential models, also referred to as force fields, for HFA134a and other fluorocarbons. Effective potentials for HFA134a and pentafluoroethane using the Halgren’s Buff 14-7 and Coulombic functions have been proposed.29,30 A molecular mechanics study of hydrofluorocarbon (up to C-6) using MM3 force field showed good agreement with experimental results
10.1021/jp806213w CCC: $40.75 2009 American Chemical Society Published on Web 12/16/2008
VLE and Interfacial Properties of HFA134a only when one or two fluorine atoms were present in the molecule.31 More recently, an all-atom model for HFA134a based on ab initio quantum mechanical calculations was developed.32 The proposed model uses the 9-6 Lennard-Jones (LJ) potential to represent the van der Waals (vdW) forces. Also, phase equilibria of mixtures containing HFA134a was predicted using statistical associating fluid theory (SAFT) with a foursite model and square-well potential.33,34 Mixtures of hydrochlorofluorocarbon/hydrofluorocarbon were described by the Stockmayer potential.35 A dipolar two-center Lennard-Jones model was also used to describe vapor-liquid equilibria of pure components and mixtures.36 Significant effort has occurred in the development of force fields for other fluorocarbons as well.37-43 However, no all-atom potential model based on a 12-6 LJ potential has been proposed for HFAs. Since the 12-6 LJ is the most frequently used potential in computer simulations to describe intermolecular interactions of water and other substances, it becomes relevant to develop a force field for HFAs that allow us to investigate interfaces with 12-6 LJ fluids.44 In our previous work,20,26 we have reported molecular dynamics (MD) simulations of the HFA|air and HFA|W interfaces. Those simulations were performed using an adaptation of a model proposed in the literature that used the Halgren’s Buff 14-7 potential.30 The model performed well in terms of predicting the surface tension of HFA134a and interfacial tension of the HFA134a-water system. However, as it will be shown here, the adapted model is not capable of accurately describing the pure component properties of HFA134a over a broad range of conditions. The ability of the adapted model to predict the properties of mixtures with HFA134a without the use of additional adjustable parameters for the unlike-pair interactions, as would be necessary when studying interfacially active species at the interface (e.g., hydrogenated surfactants), is also questionable. Different from empirically based approaches, molecular-based simulation procedures may be used to predict the phase equilibria in absence of experimental data. A variety of simulation methods,45,46 including Gibbs ensemble Monte Carlo,47 Gibbs-Duhem integration,48 pseudoensembles, and constantpressure-constant-temperature (NPT) + test particle methods45 can be used in the model development process. Usually, the histogram reweighting technique can predict the phase behavior of fluids and critical points within much smaller uncertainties.49 However, it requires more computational time than other methods since multiple runs are needed to cover the range of temperatures and densities relevant for determining a coexistence curve.46,50 In many cases, only adjustments in energy and size parameters are necessary to accurately describe the vapor-liquid coexistence curve. However, more realistic models need also to reproduce the vapor pressures. Since the acentric factor is a function of both molecule shape/size and dipole moment, the partial charges can significantly affect the prediction of the vapor pressures,51-53 and therefore, they need to be considered as one of variables during the optimization process.53,54 More recent studies have also demonstrated the relevance of including interfacial55 and transport proprieties55-57 as part of the parameter optimization algorithm. The goal of our work is to develop a new all-atom force field with a 12-6 LJ potential capable of accurately predicting both the bulk thermodynamic properties and surface tension for pure HFA134a. For the most promising HFA134a models, we also evaluated the interfacial tension against a commonly used water model. The remainder of this article is organized as follows. In the next section we provide an outline of the procedure applied
J. Phys. Chem. B, Vol. 113, No. 1, 2009 179 TABLE 1: Intramolecular Potential Energy Parameters for HFA134ac
a Parameters are assigned from the OPLS-AA force field.58 Parameters obtained by potential energy scans for a series of dihedral angles from ab initio calculations using Gaussian 03.60 c The inset in the table is the ball and stick representation of the HFA134a molecule (CF3CH2F). kB is the Boltzmann constant. b
to optimize the force field and the computational details used in the molecular dynamics (MD) simulations of the interfacial systems. The predictions of the coexistence and interfacial properties, and the comparison with the experimental observations are reported in the Results and Discussions section. Finally, the concluding remarks and future work are presented. 2. Force Field Development The intramolecular energy of the system includes a bond angle bending term (Ebend) and a torsion term (Etorsion). The length between bonded atoms is constant; i.e., they are considered rigid. The intermolecular energy is represented by the sum of the van der Waals (EvdW) and electrostatic (Eelect) interactions. The total energy of the system is given by the sum of these contributions, as shown in eq 1:
E ) Ebend + Etorsion + EvdW + Eelect
(1)
The bond lengths, bond angles, and bending constants used in this work are listed in Table 1. They are the same as those usedintheOPLS-AA(optimizedpotentialsforliquidsimulations-allatoms) force field.58 A harmonic potential was used to control bond-angle bending:
Ebend )
kθ (θ - θo)2 2
(2)
where θ is the measured bond angle, θo is the equilibrium bond angle, and kθ is the force constant. Torsion potentials were described by the following functional form:
Etorsion ) C0C1[1 + cos(φ)] + C2[1 - cos(2φ)] + C3[1 + cos(3φ)] (3) where φ and Ci are the dihedral angle and the Fourier coefficients in the cosine series function, respectively. The coefficients for the dihedral potentials for HCCF and FCCF (Ci) were determined from relaxed potential energy scans for a series of dihedral angles, similarly to the procedure employed in
180 J. Phys. Chem. B, Vol. 113, No. 1, 2009
Peguin et al.
TABLE 2: 12-6 Lennard-Jones Parameters, Partial Charges (q), and the Calculated Dipole Moment (µ) and Acentric Factor (ω) of the HFA134a Models Proposed Herec force field previous work20 1 2 3 4
qC at CF3 |e|
qC at CH2 |e|
qH |e|
εF/kB (K)
σF (Å)
qF |e|
µ (D)a
+0.9000
+0.2400
+0.0300
30.7
3.50
-0.3000
2.496
+0.5750 +0.5340 +0.6030 +0.5470
+0.0430 +0.0020 +0.0450 +0.0410
+0.1090 +0.1150 +0.1140 +0.1040
22.0 24.5 20.0 24.0
2.92 2.92 2.91 2.92
-0.2090 -0.1915 -0.2190 -0.1990
2.533 2.375 2.650 2.415
ωb
0.392 0.366 0.443 0.358
a The experimentally estimated ω for HFA134a is 0.324-0.348.65,95 The experimentally measured dipole moment of HFA134a is 2.06.97 The dipole moment calculated in this work using the visual molecular dynamics (VMD) program64 is 2.12 D. b Calculated for a bulk HFA134a molecule using eq 11. c Previous work (old): εC/kB ) 55.05 K, σC ) 3.82 Å, εH/kB ) 7.55 K, σH ) 2.58 Å. This work (model nos. 1-4): εC/kB ) 47 K, σC ) 3.6 Å, εH/kB ) 10 K, σH ) 2.5 Å. kB is the Boltzmann constant.
previous works.41,59 The angles were held fixed while the remaining internal coordinates of the molecule were optimized toward the energy minima. The Ci values are listed in Table 1. Ab initio calculations were performed using Gaussian 03.60 Relaxed potential energy scans on the necessary dihedrals were performed at the Hartree-Fock (HF) level of theory and 6-31 g+(d,p) basis set. The curve fitting toolbox (cftool) of MATLAB 7.0 was used to determine the Ci values. The vdW interactions were described by the pairwise-additive 12-6 LJ potential. Electrostatic interactions were represented by the Coulombic potential. The total nonbonded potential is given in eq 4:
Eij )
∑ ij
{ [( ) ( ) ] 4εij
σij rij
12
-
σij rij
6
+
1 qiqj 4πεo rij
}
(4)
where qi and qj are the charges centered on the individual atoms, rij is the distance between the sites i and j, and εo is the electrical permittivity of space. σij and εij are the size and energy parameters for the LJ potential, respectively. The parameters for unlike interactions were computed using standard LorentzBerthelot combining rules.61,62 The ε and σ parameters for carbon (C) and hydrogen (H) atoms were transferred from the ethane model developed previously, using the same parametrization procedure as that for HFA134a described below.63 These parameters are shown in Table 2. 2.1. Partial Charges. The procedure for obtaining the partial charges of individual atoms in HFA134a, as listed in Table 2, consisted of first fitting the electrostatic potential from geometries optimized at a given level of theory and basis sets using the CHELPG (charges from electrostatic potentials using a grid based method) subroutine in Gaussian 03.57,60 For simplicity, however, a single (same) charge was taken for all F atoms in the molecule. The charge for F was set to be the average charge of all F atoms obtained from the CHELPG analysis. For the first model (force field, FF), here called FF no. 1, the average charge value was obtained from CHELPG charges at the MP2 level of theory and 6-31+g(d,p) basis set of the HFA134a geometry optimized at MP2/6-31+g(d,p). For FF no. 2, the charges were derived at the MP2 level of theory and aug-ccpVDZ basis sets of the geometry optimized at MP2/6-31+g(d,p). Two other FFs were developed, with charges arbitrarily set to (5% those of FF no. 1. The model with 5% more negative F charge is named FF no. 3, and that with less negative charge, FF no. 4. For a given level of theory and basis set, the charges obtained from the CHELPG analysis for all H atoms were the same, thus no average was necessary. The charges for the different FFs are shown in Table 2. The coordinates of the optimized geometries for HFA134a, and the partial charges
obtained in the two different basis sets are listed in Tables S1 and S2, respectively, in the Supporting Information. The gas phase dipole moment of the HFA134a molecule for the different force fields was calculated using the visual molecular dynamics (VMD) program.64 Cartesian coordinates input to VMD were obtained from the MD simulations described below, using molecules located in the bulk region of the simulation box. The dipole moment for the HFA134a molecule optimized (geometry and CHELPG charges) in Gaussian 03 at MP2/aug-cc-pVDZ//MP2/6-31+g(d,p) (2.12 D) is also reported in Table 2. 2.2. Force Field Parametrization. Given that ε and σ for C and H atoms in HFA134a were transferred directly from the ethane model developed previously and the charges were set as discussed above, we now have ε and σ for fluorine (εF and σF) as adjustable parameters. The LJ parameters for F were initially assigned as those from the OPLS-AA force field.58 The parameters were then optimized until the calculations (discussed below) reproduced the experimentally measured phase diagram, vapor pressure, and heat of vaporization of HFA134a.65 Surface tension results were then used to further screen the developed models. 2.3. Grand Canonical Monte Carlo Simulations. Grand canonical Monte Carlo simulations (GCMC) were carried out to determine the vapor and liquid equilibrium coexistence curves, vapor pressure, and enthalpy of vaporization of the proposed FFs for HFA134a. The histogram reweighting method was used to construct the pure component phase diagram from probability distributions extracted from the simulations.49,66-68 The insertion of molecules in the GCMC simulations was enhanced through multiple first bead insertions69 and the application of the coupled-decoupled configurational-bias Monte Carlo method.70 The fractions of the various moves for each simulation were set to 15% for particle displacements, 15% for rotations, 10% for configurational-bias regrowths, and 50% for insertions and deletions. Two types of simulations were carried out. First, a cubic box of L ) 20 Å was used to obtain an optimum set of LJ parameters for F. The discussions comparing all optimized models for HFA134a reported here are based on the models obtained with this system size. Additional simulations for a larger system size (L ) 25 Å) were performed in order to ascertain the effect of system size on the parameter set.49 No finite-size effects near the critical point were observed. LJ interactions were truncated at L ) 10 Å, and standard longrange corrections were applied.71,72 An Ewald sum with tinfoil boundary conditions (kL ) 5 and kmax ) 5) was used to calculate the long-range electrostatic interactions. Simulations were equilibrated for 1 million Monte Carlo steps (MCS) before statistics were recorded. Production runs were 25 million MCS. Over the course of each simulation, the number of molecules
VLE and Interfacial Properties of HFA134a
J. Phys. Chem. B, Vol. 113, No. 1, 2009 181
(N) and energy (E) were stored in the form of a list, which was updated every 250 MCS. The necessary probability distributions were extracted from this list after the completion of the simulation. 2.4. Estimation of the Equilibrium Properties. Critical parameters for the optimized FFs were estimated by fitting the saturated liquid and vapor densities over the temperature range to the density scaling law for critical temperatures73
Fliquid - Fvapor ) B(T - Tc)β
(5)
and the law of rectilinear diameters74
Fliquid + Fvapor ) Fc + A(T - Tc) 2
(6)
where A and B are constants and β ) 0.325 is the critical exponent for Ising-type fluids in three dimensions.75 Tc and Fc are the critical temperature and density, respectively. The boiling point (Tb) was calculated via the Clausius-Clapeyron equation. Tb was obtained through linear extrapolation of that plot to 1.01 bar. The heat of vaporization (∆Hv) was calculated for each parameter set over the temperature range of 200 < T (K) < 360. Since the enthalpy is defined as H ) E + pV, ∆Hv can be calculated directly from the following:
∆Hv ) Evapor - Eliquid + p(Vvapor - Vliquid)
(7)
∆Hv is computed from the average internal energy per molecule (E), the average saturation pressure (p), and the average volume per molecule (V). The energies and molar volumes of each phase were determined from the same histogram data collected for the vapor-liquid equilibrium calculations. During the discussion of the results, the experimental and simulated values are compared and the error is defined as
error(%) ) 100
|
Micalc - Miexpt Miexpt
|
fields for molecules with no or minimum experimental data available79 and in biological systems.80 3.2. Simulation Details. The HFA134a|A interface was modeled by placing 156 molecules in the center of a rectangular box (28 × 28 × 100 Å), giving rise to two interfaces. The box size was chosen in order to represent the correct bulk density of HFA134a of 1210 kg · m-3, at 298 K and 1.5 MPa.65 The initial configuration was equilibrated in the canonical (constant number of particles, volume, and temperaturesNVT) ensemble, at 298 K, for at least 1 ns. Production runs were carried out in the microcanonical (constant number of particles, volume, and energysNVE) ensemble for 4 ns. Statistics were taken in the last 3 ns of the NVE run. The HFA134a|W interface was constructed by combining two cubic boxes of 21952 Å3 each one containing 156 HFA134a molecules and the other 735 water molecules. In this case, due to the periodic boundary conditions, the system had two HFA134a|W interfaces. The initial linear dimension in the z direction (Lz) of the HFA134a|W simulation box was chosen so that the system represented the correct bulk densities of water (1000 kg · m-3) and HFA134a (1210 kg · m-3) at 298 K and 1.5 MPa.65 Periodic boundary conditions were applied in all directions. Equilibration and production runs of the HFA134a|W interface were performed at 298 K for at least 1000 ps in the NVT and NVE ensembles, respectively. Statistics were taken in the last 500 ps of the NVE run. Standard deviations of the ensemble averages were computed by breaking the statistics into blocks of 25 ps. The simulations were carried out in parallel, using DL_POLY version 2.13.81 Equilibration was monitored through the configurational energy. The procedure adopted here is similar to previous works on both fluid|air and fluid|fluid interfaces.20,82,83 The Ewald sum was used for the longrange electrostatic potential. The cutoff for vdW interactions and the real part of the Ewald sum was 14 Å. Ewald precision was set as 10-6. A time step of 1 fs was adopted. Bond lengths were constrained using the SHAKE algorithm. All simulations in the NVT ensemble were performed using a Nose´-Hoover thermostat with relaxation constant of 0.2 ps. 3.3. Interfacial Properties. Both surface (σHFA134a) and interfacial (γHFA134a|W) tensions were calculated from the trajectories of the MD simulations discussed above using the pressure tensor equation:84
(8)
where Miexpt is the experimental value of a generic property M and Micalc is the corresponding calculated value. 3. HFA134a|Air and HFA134a|Water Interface 3.1. Molecular Models and Potentials. MD computer simulations of the HFA134a|air (HFA134a|A) and HFA134a|water (HFA134a|W) interfaces were performed using the potential models for HFA134a described here (FF nos. 1-4), and the single point charge-extended (SPCE) model for water.76 The SPCE model was chosen because it satisfactorily estimates the surface tension77 and the energetic and structural properties of water.76,77 The SPCE model was previously used in our simulations of the HFA134a|W interface.20 In the surface/ interfacial tension calculations, both HFA134a and water were treated as rigid bodies; i.e., there are no intramolecular degrees of freedom.27,78 Internal degrees of freedom of such small molecules play a small role in the prediction of thermophysical and equilibrium interfacial properties, which is the general objective of the proposed force field. However, the internal degrees of freedom become relevant in the development of force
tension )
(
)
Pxx + Pyy 1 Pzz Lz 2 2
(9)
where PRR (R ) x, y, or z) is the RR element of the pressure tensor. The factor of 1/2 in the equation arises from the two liquid|liquid (or air) interfaces in the system. An estimate of the width (w) of the HFA134a|W interface was obtained by fitting the water density profile to a hyperbolic tangent of the form
Fwater
[
( )]
z - zwo 1 o ) Fwater 1 + tanh 2 w
(10)
o is the bulk density and zow is the position of the where Fwater Gibb’s dividing surface. The “10-90” thickness,85 t, of a hyperbolic tangent is mathematically related to w by t ) 2.197w.86 The thickness is defined as the distance along the interface over which the density changes from a value of 10-90% of the water bulk density.83 The same equation was o o and FHFA134a also used to obtain the bulk densities (Fwater ) in the HFA|W system.
182 J. Phys. Chem. B, Vol. 113, No. 1, 2009
Peguin et al.
Figure 1. Vapor-liquid coexistence curves for HFA134a: experimental65 (s); models: FF no. 1 (4), FF no. 2 (3), FF no. 3 (0), and FF no. 4 (O). The error bars are smaller than the symbols.
Figure 2. Clausius-Clapeyron plot for HFA134a: experimental65 (s); models: FF no. 1 (4), FF no. 2 (3), FF no. 3 (0), and FF no. 4 (O). The error bars are smaller than the symbols.
4. Results and Discussion 4.1. Vapor-Liquid Equilibria. A comparison between the experimental65 and predicted vapor-liquid coexistence curves for HFA134a is shown in Figure 1. The absolute values of the vapor and liquid densities are reported in Tables S3 and S4, respectively, in the Supporting Information. All FFs are shown to satisfactorily represent the liquid densities, with a maximum error of 1.1%. The averaged unsigned deviations were: 0.3% for FF no. 3, 0.4% for FF nos. 1 and no. 4, and 0.7% for FF no. 2. These errors are similar to the experimental uncertainty reported in the measured density for HFA134a (in the critical region) of (0.50%.87 An error in liquid density of 2.4% was found for a previously reported model for HFA134a with the 9-6 LJ potential.32 The error in the saturated liquid densities reported for C2F6 and C2F4 models using the two-center LJ plus point quadrupole pair potential (2CLJQ), parametrized by NPT + test particle method, was FF no. 1 > FF no. 2 and no. 4ssee Table 2) as shown in Figure 2. It is expected that dispersive contributions in the intermolecular interactions between hydrofluoroalkane dimers will decrease with an increase in the dipole moment,94 thus reducing the vapor pressure of the system. While FF nos. 2 and 4 can be clearly identified as better models (in predicting vapor-liquid equilibria and vapor pressures) compared to FF nos. 1 and 3, the results so far do not allow us to discriminate between FF nos. 2 and 4. In this case, it is necessary to take in consideration how the other bulk and interfacial properties are being predicted. 4.3. Acentric Factor. The acentric factor (ω) is a parameter that measures the complexity of a molecule with respect to both its geometry and polarity,51,52 and it is directly correlated to the saturation pressure.53 The acentric factor is defined as
VLE and Interfacial Properties of HFA134a
ω ) -log10 Psat r - 1
J. Phys. Chem. B, Vol. 113, No. 1, 2009 183
(11)
where Prsat is evaluated at the reduced temperature Tr ) T/Tc ) 0.7. The ω values calculated in this work using eq 11 are reported in Table 2. The accuracy of the determination of ω depends on the quality of estimation of the Pc and Tc. The experimentally estimated ω for HFA134a is between 0.32495 and 0.348.65 Deviations in ω against the experimental value of 0.348 are: FF no. 2 (11.1%) < FF no. 4 (12.5%) < FF no. 1 (19.8%) < FF no. 3 (39.8%). Even though ω is not directly related with µ, the predicted Tc and Pc (Table 2) used in the calculation of ω are affected by µ. Figure 3 shows that the acentric factor decreases as the dipole moment and partial charges (see Table 2) decrease. The ω value is a measure of the slope of the vapor pressure curve. As the estimated ω gets larger than the experimental values, the vapor pressure curves from the models drift away from the experimental vapor pressures as seen in Figure 2. The ability of model nos. 2 and 4 in estimating ω again compare well with each other. It has been previously reported that a decrease in dipole moment (as it gets closer to the experimental value) increases the vapor pressure and decreases the acentric factor.30 4.4. Enthalpy of Vaporization. The ∆Hv was calculated for each FF, over the temperature range 200 < T (K) < 360 and are shown in Figure 4. The absolute values of ∆Hv are reported in Table S6 in the Supporting Information. All FFs showed satisfactory estimation of the ∆Hv at high temperatures. However, significant deviations from the experimental values were observed at lower temperatures. The averaged unsigned deviations ranged from 1.9% (FF no. 2) < 2.2% (FF no. 4) < 3.5% (FF no. 1) < 4.2% (FF no. 3). As reported in Table 3, FF nos. 2 and 4 best describe the ∆Hv at 298 K, with errors of 2.2%. The maximum observed error in the estimation is 5.0% for FF no. 3. Such deviations observed for FF nos. 2 and 4 are slightly higher than the errors reported in the calculations for C2F4 and C2F6 of within 1%,40 and similar to that reported for the simulation of CHF3.39 Average ∆Hv valueswereobtainedfromtheslopeoftheClausius-Clapeyron plot (Figure 2) and compared to the average experimental value (22.3 kJ · mol-1) obtained within the same temperature range. Model no. 2 is in perfect agreement with experiment, while a small deviation (0.2 kJ · mol-1) is seen for model no. 4. The average ∆Hv for model nos. 1 and 3 were 23.0 and 23.5 kJ · mol-1, respectively. It is important to note that the predictions of the ∆Hv mirror the estimations of the vapor pressures because the two of them are dependent on the dipole moment.49 One can also see that an increase in dipole moment (going away from the experimental value) accentuated the overprediction of the ∆Hv, in the same way that larger dipole moments increased the underprediction in vapor pressure. Similar effect was observed in the parametrization of other molecules.53 4.5. Critical Properties. The terms Fc and Tc for the HFA134 FFs, calculated with eqs 5 and 6, respectively, are compared against the experimental values as shown in Table 3. Tc and Fc agree very well with the experimental values. Deviations in Tc are 0.1% (FF nos. 2 and 4) and 0.2% (FF nos. 1 and 3), and those in Fc are 0.1% (FF no. 4) < 0.2% (FF no. 2) < 1.0% (FF no. 3) < 1.1% (FF no. 1). All deviations in Tc are larger than the reported experimental error of 0.02%.96 While FF nos. 2 and 4 show deviations in Fc smaller than the experimental error of 0.58%,96 larger deviations are seen for FF nos. 1 and 3. All deviations in Fc are lower than the values reported for the optimized parameters for C2F6 and C2F4 (1-2%).40
Figure 3. Acentric factor versus dipole moment: experimental (9);65,95 (b) all four models of HFA134a.
Figure 4. Heat of vaporization for HFA134a: experimental65 (s); models FF no. 1 (4), FF no. 2 (3), FF no. 3 (0), and FF no. 4 (O). Error bars are smaller than the symbols.
The Clausius-Clapeyron plot was used to determine the critical pressure (Pc) and normal boiling point (Tb) for the different FFs. Tb was obtained through linear extrapolation of the plot to 1.01 bar. Pc was determined by linear extrapolation of the vapor pressure to the critical temperature. The calculated properties and the errors are reported in Table 3. All force fields predicted a normal boiling point that is in close agreement with the experiment, with a maximum error of 0.8% (observed for FF no. 3). Larger errors were observed in the prediction of the Pc. FF nos. 2 and 4 showed the smallest deviations, with the highest deviation of 13.2% observed for FF no. 3. The errors reported in Pc (FF nos. 2 and 4) for HFA134a models developed in this work were slightly larger than the 4% reported by a series of molecules parametrized by NPT + test particle method.40 The deviations observed here are larger than the experimental deviation of 0.15%.96 As observed before for HFA134a using 2CLJD potential, an increase in the dipole moment raises the critical pressure.90 Here, the FF no. 2 best estimates the critical pressure because its dipole moment of 2.375 D is the closest (from all other force fields) to the experimental value of 2.06 D.97 On the other hand, FF no. 3 has the largest dipole moment (2.650 D) and poorly predicts the critical pressure. 4.6. Surface Tension. As an additional thermodynamic variable during the parametrization procedure, surface tension calculations were performed with the optimized force fields for HFA134a. The results are shown in Table 3. We observe that FF no. 2, with a surface tension of 7.9 ( 2.9 mN · m-1 does the best job in predicting the experimental surface tension of HFA134a, which is reported to be between 7.7 and 8.1 mN · m-1.20,98-101 The performance of FF no. 2 is closely followed by FF no. 1, with a surface tension of 8.4 ( 4.0 mN · m-1. Larger deviations are seen for FF nos. 3 and 4. When compared with the other FFs tested, we can see that the F atoms of FF no. 2 are the least negative and that the
184 J. Phys. Chem. B, Vol. 113, No. 1, 2009
Peguin et al.
TABLE 3: Experimental and Predicted Critical Density (Gc), Temperature (Tc), and Pressure (Pc), Normal Boiling Point (Tb), heat of vaporization at 298 K (∆Hv298K), and Surface Tension (σHFA134a) for HFA134ac force field
Fc (kg · m-3)
Tc (K)
Pc (bar)
Tb (K)
∆H298K (kJ · mol-1) v
expt 1 2 3 4
507.0a 501.3 (1.1) 505.8 (0.2) 501.7 (1.0) 506.5 (0.1)
374.3a 375.2 (0.2) 374.5 (0.1) 375.1 (0.2) 374.1 (0.1)
40.6a 44.2 (8.8) 43.1 (6.2) 46.0 (13.2) 43.6 (7.3)
246.9a 247.5 (0.2) 247.0 (0) 249.0 (0.8) 245.7 (0.5)
18.1a 18.8 (3.9) 18.5 (2.2) 19.0 (5.0) 18.5 (2.2)
a
σHFA134a (mN · m-1) 7.7-8.1b 8.4 ( 4.0 7.9 ( 2.9 11.0 ( 4.0 11.4 ( 3.8
Reference 65. b References 20 and 98-101. c Errors are in parenthesis (%).
dipole moment of this model is the closest to the experimental gas phase value of 2.06 D.97 The magnitude of the LJ parameters can also directly affect the calculated surface tension.55 However, in our case, we observed that an increase/decrease in εF (FF nos. 3 < 1 < 4 < 2) does not directly correlate with the surface tension value. It is worth noticing that, different from the MD simulations used to obtain the surface tension results, long tail corrections were employed in the MC calculations during the optimization of the force field parameters. At the same time that the tensions, which are determined from the pressure tensor (eq 9), directly depend on the intermolecular forces, such interactions are affected by the optimized parameters of the force field and the length of the cutoff assigned in the computer simulations.102 One exception is the surface tension obtained by MD with the lattice summation, which is independent of cutoff, but more computationally intensive.103 Since the magnitude of such corrections decreases with the increasing cutoff distance, it is expected that the use of a cutoff of 14 Å (maximum cutoff possible for the system size used in this work ) L/2 and the cutoff used in the MD work) in the MD simulations would help minimize the impact of the absence of the long tail corrections in the estimation of the interfacial properties.102 In order to probe the effect of the cutoff distance, the σHFA134a for FF no. 1 was also evaluated at cutoff of 12 Å. We noticed a slight decrease in σHFA134a from 8.4 to 7.9 mN · m-1, thus indicating that there could still be room for improvement in the calculated surface tension values. In conclusion, while all FFs were able to satisfactorily predict the VLE, saturation pressure, enthalpy of vaporization, and critical properties of HFA134a, FF nos. 2 and 4 showed the smallest deviations. Surface tension results allowed us to further screen the FFs, with the FF no. 2 best representing the experimental bulk and surface properties of HFA134a. 4.7. Comparison with Other HFA134a Models. It is instructive to compare the vapor-liquid equilibria and vapor pressures predicted by our most promising model (FF no. 2) and the force fields developed using 2CLJD model36 and the statistical associating fluid theory (SAFT).34 Although a very good prediction of the coexistence curve can be obtained using the 2CLJD model developed by Gibbs-Duhem integration,36 the vapor pressures are significantly overpredicted at lower temperatures. On the other hand, the SAFT model provides a good description of the vapor pressures,34 but the critical temperature is overpredicted, and the liquid densities at intermediate temperatures are slightly underpredicted. Different from those approaches, the most promising force field presented in this work is able to accurately predict both the coexistence curve and the vapor pressures. We also compare the results from FF no. 2 to the model previously used in the calculation of the surface tension of HFA134a.20 While the surface tension for that model (9.1 mN · m-1) is not in quite as good agreement as that of FF no. 2, it is still very close to the experimental value. The vapor-liquid
Figure 5. Comparison of the vapor-liquid coexistence curves for HFA134a obtained in this work versus (FF no. 2, 0) previous work (O).20 The experimental (s) values are also plotted.65 Error bars are smaller than the symbols.
TABLE 4: Experimental and Calculated Tension of the HFA134a|Water Interface (γHFA134a|W) and Bulk Densities (Go) Obtained with Several Optimized Force Fields and the SPCE Water Model at 298 K o o force field γHFA134a|W (mN · m-1) Fwater (kg · m-3) FHFA134a (kg · m-3)
expt 1 2 3 4 a
33.5 ( 0.3a 24.8 ( 7.8 26.3 ( 6.6 24.2 ( 4.7 25.2 ( 6.7
1000b 991.0 ( 0.4 988.2 ( 0.4 987.1 ( 0.3 988.7 ( 0.3
1202b 1210.7 ( 1.4 1210.4 ( 1.4 1210.7 ( 1.4 1211.6 ( 1.6
Reference 20. b Reference 65.
equilibria results for both models are shown in Figure 5. It is observed that, despite the fact that the force field used in that earlier study20 (adapted from Buff 14-7 potential)30 satisfactorily described the surface tension of HFA134a, it poorly describes its coexistence properties. 4.8. Tension of the HFA134a|Water Interface. As a further test of our models, we determined the interfacial tension of the HFA134a|W interface. The ability of the FF obtained for HFA134a to accurately describe the thermodynamic and microstructural properties of the neat HFA|W interface is important because that system will serve as baseline during the evaluation of surfactant activity at the HFA134a|W interface. Highly HFAphilic surfactants are expected to decrease the interfacial tension tolevelsthatfavortheformationofwater-in-HFAmicroemulsion.104,105 The results obtained for all HFA134a FFs and the SPCE water model at 298 K are summarized in Table 4. For comparison, we also include the experimental values for γHFA134a|W20 and bulk density (Fo) of both phases.65 Although the surface tension can be used as an additional property to optimize the force field parameters, the same is not expected for γHFA134a|W since the FF used for water will also affect the final result and further optimization of either force field (water or HFA134a) may lead to discrepancies in the pure component properties. A significant deviation in γHFA134a|W is observed irrespective of the FF used. However, FF no. 2 again seems to outperform
VLE and Interfacial Properties of HFA134a the other models, with a γHFA134a|W of 26.3 mN · m-1. We also notice that a decrease in dipole moment of the HFA134a FF (toward the experimental value) moves γHFA134a|W closer to the experimental interfacial tension of 33.5 mN · m-1. This is reasonable considering that an increase in HFA-water (dipoledipole) interactions is expected to decrease the tension. It can thus be expected that an HFA134a FF with a dipole closer to the experimental value (2.06 D)97 can help to further improve the prediction of the γHFA134a|W. However, it remains to be seen whether the agreement with the pure component properties can still be preserved in that case since the surface tension of pure HFA134a is being accurately predicted by FF no. 2. However, this effect (dipole moment of the HFA134a models) is not expected to be the dominant one in terms of the interfacial tension with the water. The predicted surface tension of water using the SPCE model under similar conditions is 63.6 mN · m-1,77 which is about 8.1 mN · m-1 away from the experimental value. This is approximately the difference between the calculated and experimental interfacial tension observed for the HFA134a (FF no. 2)|water (SPCE) system studied here. To verify the effect of the water model in the estimation of FFno.2 FFno.2 ), the γHFA134a|water the γHFA134a|W using the FF no. 2 (γHFA134a|water was calculated using the following equation:106 FFno.2 FFno.2 d d γHFA134a|water ) σHFA134a + σwater - 2√σHFA134a σwater -
p p 2√σHFA134a σwater (12)
where the superscripts d and p represent the dispersive and polar components of the surface tension, respectively. Using the experimental value of σwater ) 71.9 mN · m-1 at 298 K and its components σdwater ) 21.5 mN · m-1 and σpwater ) 50.4 mN · m-1,107 to assess the contribution of FF no. 2 to the calculated tension value, the computer estimated value of σHFA134a ) 7.90 mN · m-1 FFno.2 at 298 K was used for FF no. 2 (σHFA134a ). It was assumed that the dispersive and polar components of the FF contribute equally d p to the total surface energy; i.e., σHFA134a ) σHFA134a ) σHFA134a/2 -1 ) 3.95 mN · m . Substituting those values in eq 12, a -1 γFFno.2 can be calculated. It is interesting HFA134a|water ) 33.15 mN · m to notice that this value is in very close agreement to the experimental value of 33.50 mN · m-1. This suggests that, provided an accurate model for the surface tension of water is obtained, the interfacial tension should be correctly described with the FF no. 2 developed here. 4.9. Microstructure of the HFA134a|Water Interface. Figure 6a is a snapshot from an equilibrium MD configuration of the HFA134a|W interface. Two copies of the same simulation box are shown next to each other. The equilibrium density profiles for HFA134a (FF no. 2) and water (SPCE) as a function of the z coordinate are shown in Figure 6b. The formation of a stable interface is the first indicator of the appropriateness of the models. Protrusions of water into HFA134a and vice versa can be qualitatively observed. The density profile of water is smooth and fast decaying, while that of HFA134a exhibits some oscillations. The characteristics of the interfacial region of the HFA134a|W system observed here are typical of organic|water systems, showing a sharp but smooth transition from one liquid phase to the other.108,109 The calculated interfacial thickness (t) for the HFA134a|W interface obtained for FF no. 2 is 4.9 Å. This values is somewhat smaller than that observed in our previous work of 5.2 Å.20 The bulk densities of HFA134a and water are reported in Table 4 and were obtained by fitting the density profile to eq
J. Phys. Chem. B, Vol. 113, No. 1, 2009 185
Figure 6. (a) Snapshot of an equilibrium configuration of the HFA134a(FF no. 2)|water(SPCE) simulation box. Two copies of the same snapshot are plotted next to each other. (b) z-Dependent density profile for the HFA134a|W interface. Error bars for density are smaller than the symbols.
10. The experimental bulk density of water at 298 K and saturation conditions is 1000 kg.m-3.65 The experimental density for HFA134a, at the same conditions, is 1202 kg.ml-3.65 While all simulations slightly underpredict the water density, they overpredict (to a smaller degree) the HFA density. The deviations in density are as follows: for SPCE water in the presence of FF no. 1 (0.9%) < FF no. 4 (1.1%) < FF no. 2 (1.2%) < FF no. 3 (1.3%) and for HFA134a 0.7% (FF nos. 1, 2, and 3) < 0.8% (FF no. 4). 5. Conclusions In this work, new all-atom force fields for HFA134a that use the 12-6 LJ-Coulombic potential to describe the nonbonded interactions have been investigated. Monte Carlo simulations and histogram reweighting technique were used during the parametrization procedure and in the calculation of the phase coexistence properties for the HFA134a models. Surface tension was used as an additional variable in the screening of the potential models. The best performing model, FF no. 2, was able to satisfactorily reproduce the liquid (average deviation of 0.7% compared with experimental data) and vapor densities (average deviation of 4.4%); vapor pressure (average deviation of 3.1%); critical density (deviation of 0.2%), critical temperature (deviation of 0.1%), and critical pressure (deviation of 6.2%); boiling temperature (exact); heat of vaporization at 298 K (deviation of 2.2%); and surface tension (within the experimental range of 7.7-8.1 mN · m-1). Overall, the model proposed here is shown to be a better descriptor of the experimental properties of HFA134a compared to those previously discussed in the literature.20,34,36,87 It also represents the first all-atom 12-6 LJ model for HFA134a. The potential model was further tested by investigating the interfacial tension of HFA134a (FF no. 2) against water (SPCE). The interfacial tension was seen to deviate 21.5% from the experimental value. The deviation was attributed to the inability of the SPCE water model to accurately describe the surface tension of water77 and to a lesser extent to the fact that the dipole for FF no. 2 was slightly larger than the experimental gas phase value, thus enhancing the water-HFA interactions (decreasing tension). Acknowledgment. Financial support for this work has been provided by the National Science Foundation, under CBET0553537 (S.R.P.R.) and CBET-0522005 (J.J.P.) grants, and the office of the OVPR at WSU for a Nano@WSU grant (S.R.P.R.). We would like to acknowledge a graduate research assistantship for R.P.S.P from IMR at WSU and GRID/WSU for computer time.
186 J. Phys. Chem. B, Vol. 113, No. 1, 2009 Supporting Information Available: Cartesian coordinates and partial charges for the optimized HFA134a at an MP2 level of theory; predicted liquid and gas densities, vapor pressure, heat of vaporization, and their deviations from the experimental values. This information is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Opportunities to Minimise Emissions of Hydrofluorocarbons (HFCs) from the European Union; March Consulting Group: UK 1998; p 119. (2) Montreal Protocol on Substances that Deplete the Ozone Layer, United Nations EnVironment Programme; UNEP: Nairobi, Kenya, 2000; p 52. (3) Climate Change - Emissions of Industrial Greenhouse Gases (HFCs, PFCs and Sulphur Hexafluoride), Environmental Protection Agency: Wexford, Ireland, 2000; p 81. (4) Dieckmann, J. T.; Magid, H. Global ComparatiVe Analysis of HFC and AlternatiVe Technologies for Refrigeration, Air Conditioning, Foam, SolVent, Aerosol Propellant, and Fire Protection Applications; Arthur D. Little, Inc.: Cambridge, MA, 2002; p 150. (5) Terzano, C. Pulm. Pharmacol. Ther. 2001, 14351, 366. (6) Graepel, P.; Alexander, D. J. J. Aerosol Med. 1991, 4 (3), 193– 200. (7) Aiache, J. M. Ann. Pharm. Fr. 1997, 55 (2), 62–68. (8) Harrison, L. I.; Soria, I.; Cline, A. C.; Ekholm, B. P. J. Pharm. Pharmacol. 1999, 51 (11), 1235–1240. (9) Hofford, J. M. J. Fam. Practice 1992, 34 (4), 485–492. (10) Keller, M. Int. J. Pharm. 1999, 186 (1), 81–90. (11) Peguin, R. P. S.; Wu, L.; da Rocha, S. R. P. Langmuir 2007, 23 (16), 8291–8294. (12) Jones, S. A.; Martin, G. P.; Brown, M. B. J. Controlled Release 2006, 115 (1), 1–8. (13) Liao, Y.-H.; Brown, M. B.; Jones, S. A.; Nazir, T.; Martin, G. P. Int. J. Pharm. 2005, 304 (1-2), 29–39. (14) Stefely, J. S. Drug DeliV. Technol. 2002, 2 (662), 64–69. (15) Wu, L.; Al-haydari, M.; da Rocha, S. R. P. Eur. J. Pharm. Sci. 2008, 33 (2), 146–158. (16) Wu, L.; Bharatwaj, B.; Panyam, J.; da Rocha, S. R. P. Pharm. Res. 2008, 25 (2), 289–301. (17) Selvam, P.; Peguin, R. P. S.; Chokshi, U.; da Rocha, S. R. P. Langmuir 2006, 22 (21), 8675–8683. (18) Selvam, P.; Chokshi, U.; Gouch, A.; Wu, L.; Porcar, L.; da Rocha, S. R. P. Soft Matter 2008, 4357, 366. (19) Wu, L.; Peguin, R. P. S.; da Rocha, S. P. R. J. Phys. Chem. B 2007, 111 (28), 8096–8104. (20) Peguin, R. P. S.; Selvam, P.; da Rocha, S. R. P. Langmuir 2006, 22 (21), 8826–8830. (21) Patel, N.; Marlow, M.; Lawrence, M. J. J. Colloid Interface Sci. 2003, 258 (2), 345–353. (22) Patel, N.; Marlow, M.; Lawrence, M. J. J. Colloid Interface Sci. 2003, 258 (2), 354–362. (23) Steytler, D. C.; Thorpe, M.; Eastoe, J.; Dupont, A.; Heenan, R. K. Langmuir 2003, 19 (21), 8715–8720. (24) Jackson, K.; Fulton, J. L. Langmuir 1996, 12 (22), 5289–5295. (25) Williams, R. O.; Liu, J. Eur. J. Pharm. Sci. 1999, 7 (2), 137– 144. (26) Wu, L.; Peguin, R. P. S.; Selvam, P.; Chokshi, U.; da Rocha, S. P. R. Molecular Scale Behavior in Alternative Propellant-Based Inhaler Formulations. In Inhalation Aerosols - Physical and Biological Basis for Therapy; 2nd ed.; Hickey, A. J., Ed.; Informa Healthcare USA: New York, 2007; pp 373-397. (27) da Rocha, S. R. P.; Johnston, K. P.; Rossky, P. J. J. Phys. Chem. B 2002, 106 (51), 13250–13261. (28) da Rocha, S. R. P.; Johnston, K. P.; Westacott, R. E.; Rossky, P. J. J. Phys. Chem. B 2001, 105 (48), 12092–12104. (29) Lisal, M.; Budinsky, R.; Vacek, V.; Aim, K Int. J. Thermophys. 1999, 20 (1), 163–174. (30) Lisal, M.; Vacek, V. Fluid Phase Equilib. 1997, 127 (1-2), 83– 101. (31) Chen, K. H.; Walker, G. A.; Allinger, N. L. THEOCHEM 1999, 49087, 107. (32) Fermeglia, M.; Ferrone, M.; Pricl, S. Fluid Phase Equilib. 2003, 210 (1), 105–116. (33) Galindo, A.; Gil-Villegas, A.; Whitehead, P. J.; Jackson, G.; Burgess, A. N. J. Phys. Chem. B 1998, 102 (39), 7632–7639. (34) Galindo, A.; Whitehead, P. J.; Jackson, G.; Burgess, A. N. J. Phys. Chem. B 1997, 101 (11), 2082–2091. (35) Gao, G. T.; Wang, W.; Zeng, X. C. Fluid Phase Equilib. 1999, 158, 16069–78.
Peguin et al. (36) Budinsky, R.; Vacek, V.; Lisal, M. Fluid Phase Equilib. 2004, 222, 223213–220. (37) Bonifacio, R. P.; Padua, A. A. H.; Gomes, M. F. C. J. Phys. Chem. B 2001, 105 (35), 8403–8409. (38) Palmer, B. J.; Anchell, J. L. J. Phys. Chem. 1995, 99 (32), 12239– 12248. (39) Potter, S. C.; Tildesley, D. J.; Burgess, A. N.; Rogers, S. C. Mol. Phys. 1997, 92 (5), 825–833. (40) Vrabec, J.; Stoll, J.; Hasse, H. J. Phys. Chem. B 2001, 105 (48), 12126–12133. (41) Watkins, E. K.; Jorgensen, W. L. J. Phys. Chem. A 2001, 105 (16), 4118–4125. (42) Yamamoto, R.; Kitao, O.; Nakanishi, K. Mol. Simul. 1994, 12 (3-6), 383–391. (43) Yamamoto, R.; Kitao, O.; Nakanishi, K. Fluid Phase Equilib. 1995, 104349, 361. (44) Benjamin, I. Annu. ReV. Phys. Chem. 1997, 48407, 451. (45) Lotfi, A.; Vrabec, J.; Fischer, J. Mol. Phys. 1992, 76 (6), 1319– 1333. (46) Panagiotopoulos, A. Z. J. Phys.: Condens. Matter 2000, 12 (3), R25–R52. (47) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61 (4), 813–826. (48) Kofke, D. A. J. Chem. Phys. 1993, 98 (5), 4149–4162. (49) Potoff, J. J.; Panagiotopoulos, A. Z. J. Chem. Phys. 1998, 109 (24), 10914–10920. (50) Potoff, J. J.; Siepmann, J. I. AlChE J. 2001, 47 (7), 1676–1682. (51) Pitzer, K. S. J. Am. Chem. Soc. 1955, 773427, 3433. (52) Pitzer, K. S.; Lippmann, D. Z.; Curl, R. F., Jr.; Huggins, C. M.; Petersen, D. E. J. Am. Chem. Soc. 1955, 77, 3433–3440. (53) Ketko, M. B. H.; Potoff, J. J. Mol. Simul. 2007, 33 (9-10), 769– 776. (54) Kamath, G.; Lubna, N.; Potoff, J. J J. Chem. Phys. 2005, 123, 124505/124501–124505/124507. (55) Olivet, A.; Vega, L. F J. Chem. Phys. 2007, 126 (14), 144502/ 144501–144502/144511. (56) Zhang, H.; Ely, J. F. Fluid Phase Equilib. 2004, 217 (1), 111– 118. (57) Maginn, E. J. Acc. Chem. Res. 2007, 40 (11), 1200–1207. (58) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118 (45), 11225–11236. (59) Padua, A. A. H. J. Phys. Chem. A 2002, 106 (43), 10116–10123. (60) Frisch, M. J. T., G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision A, 1st ed.; Gaussian, Inc.: Pittsburgh, 2003. (61) Berthelot, D. C.R. Hebdomadaires Seances Acad. Sci. 1898, 1261703. (62) Lorentz, H. A. Ann. Phys. 1881, 12127. (63) Potoff, J. J. internal communication;Wayne State University, 2007. (64) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graph. 1996, 14 (1), 33–38; plates 27-28. (65) Lemmon, E. W.; McLinden, M. O.; Huber, M. L. REFPROP Reference Fluid Thermodynamic and Transport Properties, 7.0 ed.; National Institute of Standards and Technology: Gaithersburg, MD, 2002. (66) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1988, 61 (23), 2635–2638. (67) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1989, 63 (12), 1195–1198. (68) Potoff, J. J.; Errington, J. R.; Panagiotopoulos, A. Z. Mol. Phys. 1999, 97 (10), 1073–1083. (69) Esselink, K.; Loyens, L. D. J. C.; Smit, B. Phys. ReV. E: Stat. Phys., Plasmas, Fluids 1995, 51 (2), 1560–1568. (70) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1999, 103 (21), 4508–4517. (71) McDonald, I. R. Mol. Phys. 1972, 23 (1), 41–58. (72) Wood, W. W. J. Chem. Phys. 1968, 48, (1), 415-434. (73) Rowlinson, J. S.; Swinton, F. L. Liquids and Liquid Mixtures, 3rd ed.; Butterworth: London, 1982. (74) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Clarendon Press: Oxford, 1982.
VLE and Interfacial Properties of HFA134a (75) Privman, V. Encyclopedia of Applied Physics; Wiley-VCH: Berlin, 1998; Vol. 23. (76) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91 (24), 6269–6271. (77) Vega, C.; de Miguel, E. J. Chem. Phys. 2007, 126 (15), 154707/ 154701–154707/154710. (78) Dang, L. X. J. Chem. Phys. 1999, 110 (20), 10113–10122. (79) Maple, J. R.; Hwang, M. J.; Stockfisch, T. P.; Dinur, U.; Waldman, M.; Ewig, C. S.; Hagler, A. T. J. Comput. Chem. 1994, 15 (2), 162–182. (80) Kitson, D. H.; Avbelj, F.; Moult, J.; Nguyen, D. T.; Mertz, J. E.; Hadzi, D.; Hagler, A. T. Proc. Natl. Acad. Sci. U.S.A. 1993, 90 (19), 8920– 8924. (81) Forester, T. R.; Smith, W. DL-POLY Package of Molecular Simulations; 2.13 ed.; CCLRC, Daresbury Laboratory: Daresbury, Warrington, England, 2001. (82) Dominguez, H.; Berkowitz, M. L. J. Phys. Chem. B 2000, 104 (22), 5302–5308. (83) Rivera, J. L.; McCabe, C.; Cummings, P. T. Phys. ReV. E: Stat., Nonlin. Soft Matter Phys. 2003, 67 (1-1), 011603/011601–011603/011610. (84) Croxton, C. A. Statistical Mechanics of the Liquid Surface; John Wiley & Sons: Chichester: New South Wales, Australia, 1980. (85) Sides, S. W.; Grest, G. S.; Lacasse, M.-D. Phys. ReV. E: Stat. Phys., Plasmas, Fluids 1999, 60 (6-A), 6708–6713. (86) Chapela, G. A.; Saville, G.; Thompson, S. M.; Rowlinson, J. S. J. Chem. Soc., Faraday Trans. 1977, 73 (8), 1133–1144. (87) Kabata, Y.; Tanikawa, S.; Uematsu, M.; Watanabe, K Int. J. Thermophys. 1989, 10 (3), 605–615. (88) Kamath, G.; Georgiev, G.; Potoff, J. J. J. Phys. Chem. B 2005, 109 (41), 19463–19473. (89) Kamath, G.; Cao, F.; Potoff, J. J. J. Phys. Chem. B 2004, 108 (37), 14130–14136. (90) Lisal, M.; Aim, K.; Fischer, J. Mol. Simul. 2000, 23 (6), 363– 388. (91) Leeuwen;, D. J. v.; Hermens, J. L. M. Risk Assessment of Chemicals; Kluwer Academic: Boston, 1995; p 374.
J. Phys. Chem. B, Vol. 113, No. 1, 2009 187 (92) Jaksland, C. A.; Gani, R.; Lien, K. M. Chem. Eng. Sci. 1995, 50 (3), 511–530. (93) Tzou, T.-Z.; Pachuta, R. R.; Coy, R. B.; Schultz, R. K. J. Pharm. Sci. 1997, 86 (12), 1352–1357. (94) Tsuzuki, S.; Uchimaru, T.; Mikami, M.; Urata, S. J. Phys. Chem. A 2003, 107 (39), 7962–7968. (95) Sheikh, S.; Papari, M. M.; Boushehri, A. Ind. Eng. Chem. Res. 2002, 41 (13), 3274–3281. (96) Marsh, K. N.; Abramson, A.; Ambrose, D.; Morton, D. W.; Nikitin, E.; Tsonopoulos, C.; Young, C. L. J. Chem. Eng. Data 2007, 52 (5), 1509– 1538. (97) HFA Propellants for Medical Use; Solvay Fluor und Derivate GmbH: Solvay, Hannover, 2003; p 43. (98) Chae, H. B.; Schmidt, J. W.; Moldover, M. R. J. Chem. Eng. Data 1990, 35 (1), 6–8. (99) Froba, A. P.; Will, S.; Leipertz, A Int. J. Thermophys. 2000, 21 (6), 1225–1253. (100) Heide, R Int. J. Refrig. 1997, 20 (7), 496–503. (101) Higashi, Y.; Shibata, T.; Okada, M. J. Chem. Eng. Data 1997, 42 (3), 438–440. (102) Shen, V. K.; Mountain, R. D.; Errington, J. R. J. Phys. Chem. B 2007, 111 (22), 6198–6207. (103) Lopez-Lemus, J.; Alejandre, J. Mol. Phys. 2002, 100 (18), 2983– 2992. (104) da Rocha, S. R. P.; Harrison, K. L.; Johnston, K. P. Langmuir 1999, 15 (2), 419–428. (105) da Rocha, S. R. P.; Johnston, K. P. Langmuir 2000, 16 (8), 3690– 3695. (106) Good, R. J.; Girifalco, L. A. J. Phys. Chem. 1960, 64561, 565. (107) Binks, B. P. Langmuir 1993, 925, 28. (108) Benjamin, L. J. Chem. Phys. 1992, 97 (2), 1432–1445. (109) Chang, T.-M.; Dang, L. X. J. Chem. Phys. 1996, 104 (17), 6772–6783.
JP806213W
