J. Phys. Chem. C 2008, 112, 5029-5035
5029
Surface Tension and Surface Orientation of Perfluorinated Alkanes Mesfin Tsige* Department of Physics, Southern Illinois UniVersity, Mail Code 4401, Carbondale, Illinois 62901
Gary S. Grest Sandia National Laboratories, Albuquerque, New Mexico 87185 ReceiVed: NoVember 7, 2007; In Final Form: January 18, 2008
Molecular dynamics simulations using two all-atom force fields, the OPLS force field of Jorgensen et al. and the recently developed force field of Borodin et al. (exp-6), have been used to study the equilibrium liquidvapor interface properties of perfluorinated alkanes. Both force fields predict similar liquid densities and are in excellent agreement with existing experimental data for short chains. The OPLS force field offers better agreement with experimental data for surface tension at lower molecular weights than the exp-6 force field, which overpredicts it by as much as 9%. However, for longer chains the OPLS force field predicts higher surface tension than the exp-6 force field. Also, the OPLS substantially overpredicts the melting temperature confirming a recent suggestion that this force field may not be transferable to longer perfluorinated alkane chains. The orientation of the chains at the surface is found to be strongly dependent on both chain length and temperature. On average, the chain segments tend to orient perpendicular to the surface with the -CF3 end segments at the surface.
I. Introduction Perfluorinated alkanes with a molecular formula CnF2n+2 have outstanding properties, very different from those of normal hydrocarbons that make them ideal candidates for a variety of applicationsrangingfromlubricantstoartificialbloodsubstitutes.1-16 Their chemical structure and the weak intermolecular interactions are mainly responsible for their unique properties.17 Because of their low surface tension and refractive index, high viscosity, and gas solubility, most of the technological applications of perfluorinated alkanes occur at the interface. It is common to introduce perfluorinated segments into polymers in order to cause significant changes in their surface properties. Therefore, predicting the surface properties of fluorinated materials before their synthesis can provide the means to tailor their surface properties for specific applications. Understanding the structure and thermodynamic properties of the free surface of fluorocarbons is thus an area of fundamental and current interest with numerous technological implications.18 The properties of the free surface are in general governed by surface tension, which in turn is related to the difference between the intermolecular interaction in the bulk and at the surface.19 Thus, a theory that exactly describes intermolecular interactions is essential to predict the properties of the free surface accurately. To date, few theoretical studies of the free surface of fluorocarbons exist despite their great promise for nanotechnological and biomedical applications.18 From both theoretical and experimental perspectives, our knowledge of the surface of fluorocarbons is very limited. Simulation studies of fluorocarbon systems have been less well-developed than their hydrocarbon counterparts mainly because of the lack of accurate force fields in the past that are capable of describing the complex conformational nature of * Corresponding author. E-mail:
[email protected].
these systems. However, reliable ab inito-based force fields have been proposed recently for perfluorinated alkanes.6,10,20 The properties of fluorocarbons are principally controlled by their chemical structure and intermolecular interactions. Consequently, a force field that describes intermolecular interaction accurately and also takes into account the detailed chemical structure of the molecule is necessary for reliable simulations of these systems.12 A few computational studies have been conducted that qualitatively describe the interactions of fluorocarbon systems.6,10,21-26 An early united-atom model,21-23 in which each carbon and its attached fluorines were treated as a single “atom”, developed to study the packing structure of monolayers of partially fluorinated amphiphiles supported on water, was able to qualitatively capture a few experimentally observed properties. However, when the model was used for perfluorinated alkanes,27 it underestimated the saturated liquid density and failed to predict the vapor-liquid equilibria because of a poor description of both intramolecular and intermolecular interactions. Later united models24,26 suggested that reliable unitedatom force fields can be derived from ab initio quantumchemical calculations. The united-atom force field of Cui et al.,24 which is based on density functional and ab initio calculations and optimized for fluid phase equilibrium of fluorinated alkanes, has been successful in describing pure-liquid properties such as saturation pressures, densities, and vapor-liquid equilibria. This force field was also found to predict phase equilibria accurately in fluorinated alkane/carbon dioxide mixtures.4 However, it was later found that this model is inadequate for predicting transport properties such as the viscosity and diffusivity of perfluorinated alkanes.12 Fluorocarbons are found to be very sensitive to molecular details and an all-atom description of the molecular structure, where the fluorine atoms are treated explicitly, may be required in order to capture, in principle, all of the required properties.
10.1021/jp710678w CCC: $40.75 © 2008 American Chemical Society Published on Web 03/07/2008
5030 J. Phys. Chem. C, Vol. 112, No. 13, 2008
Tsige and Grest
TABLE 1: Force Field Parameters for Perfluorinated Alkanes from the exp-6 Model. Taken from Reference 20 Uexp-6(rij) ) Aij exp(-Bijrij) + Cijr-6 ij nonbonded pair
A (kcal/mol)
B (Å-1)
C (kcal Å6/mol)
C-C C-F F-F
165 734.9 38 941.5 9693.8
3.8867 3.7155 3.5861
332.4 233.6 164.2
0 2 uB(rij) ) 1/2kbond ij (rij - rij)
bond
kbond ij
C-C CF3-F CF2-F
618.0 722.0 722.0
r0ij
(kcal/mol)
(Å)
1.560 1.339 1.351 ang 0 2 uA(θijk) ) 1/2kijk (θijk - θijk )
angle C-C-C C-C-F F-CF3-F F-CF2-F
ang kijk
θ0ij (degrees)
(kcal/mol)
160.6 180.0 240.0 240.0
115.6 109.0 109.1 110.3 7 tor uT(φijkl) ) 1/2 ∑n)1 kijkl (n)[1 - cos(nφijkl)]
torsion (kcal/mol)
tor kijkl (1)
tor kijkl (2)
tor kijkl (3)
tor kijkl (4)
tor kijkl (5)
tor kijkl (6)
tor kijkl (7)
C-C-C-C F-C-C-F
-1.463 0.000
-0.722 0.000
1.413 -0.480
-1.1415 0.000
-0.220 0.000
0.000 0.000
0.761 0.000
An advantage of using an all-atom representation of the molecules rather than a united-atom lies in the ability to examine the role of molecular flexibility and internal vibrations in affecting the molecular orientation in the system. Because the fluorine atoms are treated explicitly in all-atom models, the requirement in terms of computer time is significantly larger than that for united-atom models. This has been yet another reason for the limited number of atomistic simulation studies of fluorocarbons. The development of parallel molecular dynamics (MD) algorithms and the availability of increased computational power have now reached the point where detailed atomistic simulation study of large fluorocarbon systems is possible. Currently, there are two commonly used atomistic force fields optimized for simulation of perfluorinated alkanes: the Optimized Parameter for Liquid Simulation-All Atom (OPLSAA) force field of Watkins and Jorgensen6 (hereafter referred to as OPLS force field) and the force field of Borodin et al. (hereafter referred to as exp-6 force field).10,20 Both force fields are found to produce accurate descriptions of densities and heat of vaporization for perfluorinated alkanes. However, the OPLS force field is found to be less accurate in describing conformational energies of short chains,10 thus bringing into question the transferability of the force field to longer chains. The transferability of the exp-6 potential to longer chains has been demonstrated.10 At present the exp-6 force field is expected to perform better than the OPLS force field in describing structural and dynamic properties of a range of fluorinated alkane chains. The main objective of the present work is to elucidate the molecular basis for the structure and surface tension properties of the free surfaces of fluorocarbons. Atomistic MD simulations using both the OPLS and exp-6 force fields were performed to predict the surface properties of perfluorinated alkanes. This allows us to compare results from both force fields to each other and to available experimental results. Details of the models used and the simulation procedure are given in Section II. In Section III, simulation results from OPLS and exp-6 force fields are presented and discussed. Specifically, simulation results from the two force fields for density profiles, surface tensions, and orientations of chains at and away from
the interface are presented and compared with each other and with existing experimental results. Our findings are summarized in Section IV. II. Model and Simulation Details A. Force Field. In this work, we use two sets of all-atom force fields, OPLS and exp-6. Both of these force fields represent the total potential energy of the system as a sum of four types of potentials: nonbonded interactions, bond and angle stretching, and torsion potentials. The main difference between the two force fields is the way the repulsion term in the nonbonded interaction is modeled. OPLS treats the nonbonded interactions via a combination of a Lennard-Jones 6-12 potential functions and a Coulomb potential as
[( ) ( ) ]
UOPLS(rij) ) 4ij
σij rij
12
-
σij rij
6
+
qiqj rij
(1)
where rij is the distance between atoms i and j, ij and σij are the standard Lennard-Jones parameters, and qi and qj are the charges on atoms i and j. In the exp-6 force field the nonbond interactions are represented by a combination of a Buckingham potential (exponential-6) and a Coulomb potential20
[
Uexp-6(rij) ) Aij exp(-Bij rij) -
Cij rij6
]
+
qiqj rij
(2)
Parameters Aij, Bij, and Cij are given in Table 1 for convenience. In both force fields, nonbond interactions are evaluated for all intermolecular interactions and for intramolecular interactions between atoms separated by three or more bonds, that is, starting with 1-4 interactions. The 1-4 intramolecular interactions are reduced by a factor of 2 in the OPLS force field. In the OPLS force field a charge of -0.12 e, 0.24 e, and 0.36 e is assigned to F atoms, to fluromethylene C atoms, and to fluromethyl C atoms, respectively, while in the exp-6 force field charges of -0.1176 e, 0.2352 e, and 0.3528 e are assigned, respectively.
Surface Tension and Surface Orientation
J. Phys. Chem. C, Vol. 112, No. 13, 2008 5031
In both models, bond bending and bond stretching interactions are modeled as harmonic potentials 0 2 UB(rij) ) kbond ij (rij - rij)
(3)
0 2 UA(θijk) ) kang ijk (θijk - θijk)
(4)
and
where rij is the bond distance between atoms i and j, θijk is the angle formed by atoms i, j, and k, and r0ij and θ0ijk denote the corresponding equilibrium values with spring constants kbond ij and kang ijk , respectively. The torsional interactions for the OPLS model are defined by tors (φijkl) ) UOPLS
∑n kn cosn(φijkl)
(5)
TABLE 2: Parameters for the Different Systems Studieda
while for the exp-6 model they are defined by tors (φijkl) ) Uexp-6
1 2
∑n kn[1 - cos(nφijkl)]
Figure 1. Sample 2D view of the C10F22 system at the end of 7 ns run at 348 K. The left and right sides of the figure represent the surfaces of the system, and only a few of the 200 molecules are in the vapor phase.
(6)
where in both cases φijkl is the dihedral angle and kn are torsional parameters. There are several versions of the exp-6 force field parameters.10,20,28 For the present study, we used the latest parameters from ref 20 and for convenience they are given in Table 1. In this study, perfluorinated alkanes of short chain lengths C6F14-C10F22, intermediate chain length C20F42, and relatively long chain lengths C48F98 and C96F194 were studied. To generate initial configurations for the free surface simulations of these different chain lengths, we first generated bulk melts using TOWHEE,29 a Monte Carlo simulation code. Briefly, for a given perfluorinated alkane the TOWHEE code initially constructed a layered crystal with all chains in the all-trans conformation. Several different Monte Carlo moves were then used to sample the phase space at a temperature much higher than the melting temperature but lower than the boiling temperature of a given perfluorinated alkane. The equilibrated bulk configuration of the TOWHEE output for each perfluorinated alkane was further equilibrated using MD simulations with constant number of molecules, pressure, and temperature (NPT) at 1 atm with periodic boundary conditions in all directions for about 400 ps. After the NPT equilibration run, the simulation cell in the z direction was extended to generate the required free surfaces where the interface lies in the (x,y) plane. The increased system size in the z direction, as can be seen in the sample configuration shown in Figure 1, was large enough so that the long-range Coulombic interactions between perfluorinated alkane molecules at opposite surfaces was negligible. The system sizes, that is, number of chains in a simulation box for each perfluorinated alkane, were chosen to yield a liquid film of reasonable thickness, Lz, and lateral dimensions, Lx ) Ly, that were much larger than the cutoff distance so that the long-range Coulombic interaction between a chain and its mirror image in all periodic directions was minimal. Depending on the chain length of the perfluorinated alkane, the total number of molecules in the simulation cell was between 100 and 320 molecules. System sizes and dimensions of the simulated systems are reported in Table 2. For each configuration, MD simulations in the NVT ensemble using both OPLS and exp-6 force fields were performed at a range of temperatures above the experimental bulk melting
n N Lx (Å) Lz (Å)
6 270 45.29 97.29
7 200 42.56 89.43
8 320 50.94 99.08
9 320 52.48 103.89
10 200 47.52 97.29
20 100 47.42 110.29
48 100 66.21 121.58
96 100 78.00 241.62
a n is the number of carbon atoms per chain, N is the total number of chains in a system, Lx ) Ly are the surface dimensions, and Lz is simulation box dimension in the z direction.
temperature but below the boiling point of a given perfluorinated alkane system. All of the simulations were done using the LAMMPS simulation package.30 The equations of motion were integrated using the Verlet algorithm, and the temperature during the simulations was kept constant by using a Langevin thermostat with a damping constant of 0.01 fs-1. The effect of cutting off the dispersion terms on the results of the simulations was initially investigated by varying the cutoff radius rc between 10 and 16 Å. On the basis of the result of this investigation, discussed in the next section, we used rc ) 16 Å for the rest of the simulation. The Coulombic terms were calculated using the particle-particle/particle-mesh (PPPM) algorithm.31 An integration time step of 1 fs was used in all simulations. The minimum simulation time was 7 ns for short chains and 8 ns for intermediate and long chains; some systems were run as long as 12 ns. B. Surface Tension. The surface tension (γp) was calculated during the simulation by applying the Kirkwood-Buff formula32 to the geometry of our system where γp is expressed as an ensemble average33
γp )
[
]
Lz 〈 px〉 + 〈 py〉 〈 pz〉 2 2
(7)
where pz and (px + py)/2 are the normal and tangential pressure components, respectively. The factor of 1/2 outside the bracket takes into account the presence of two interfaces in the system. The surface tension, γp, calculated using eq 7 does not take into account the full nonbonded potential because interactions within only a cutoff sphere are considered. To estimate the total surface tension γ, we include a tail correction term γtail that represents contribution of the potential to the surface tension for distances longer than the cutoff. Tail correction was required for the Lennard-Jones interaction only because the long-range electrostatic interactions were handled by the PPPM method. The tail correction was determined from34,35
5032 J. Phys. Chem. C, Vol. 112, No. 13, 2008
γtail )
π 2
∫-∞∞ ∫-11 ∫r∞ r3 c
Tsige and Grest
dU(r) g(r)(1 - 3s2)(F(z)F(z dr sr) - (FG(z))2)dr ds dz (8)
where U(r) is the dispersion (r-6) portion of the pairwise potentials in eqs 1 and 2, g(r) is the radial distribution function, F(z) is the observed density profile, rc is the cutoff radius for the dispersion terms, and FG(z) is a Gibbs dividing surface
FG(z) ) Fc +
∆F sgn(z) 2
(9)
where Fc ) (FL + FV)/2 is the average density of the liquid and vapor phases and ∆F ) FL - FV is the difference between the average densities of the two phases. On the basis of previous simulation study of capillary waves in water,33 the density profile for the tail correction is approximated by an error function, as was shown by Lacasse et al.36 Figure 2 shows typical example of running time average of the surface tension, without tail correction, for two systems during simulations performed using the OPLS force field. The left side of the figure is for C9F20 at 323 K and right side is for C20F42 at 450 K. C9F20 at 323 K is a liquid and the simulation result indicates that at least 6 ns of simulation time is required in order to obtain a reliable value for its surface tension. This partly explains why previous simulations,25 on the order of 1-2 ns sampling, did not accurately determine the surface tension of short perfluorinated alkane chains. Great care, however, should be taken in interpreting the data for the C20F42 system. The sudden rise in surface tension at around 5 ns is an indication that the surface structure is experiencing a configurational change. Direct visualization of the system at 5 ns and beyond shows that the sudden rise is due to the system going to a crystalline state. The experimental melting temperature of C20F42 is 437 K, but we found that the OPLS force field overpredicted the melting temperature and the system crystallized at temperatures as high as 460 K. In general, our observation is that the deviation between the experimental and simulation melting temperatures using the OPLS force field increased as the chain length increased, which is an indication that the OPLS force field was not optimized for long chains. As a result, in order to determine the surface tension accurately the simulations were run for at least 7 ns after equilibration. It is important to mention that the simulation results using the exp-6 force field did not show any sign of crystallization within the range of temperatures we considered for each chain length. This range is above the experimental melting temperature of the system. However, surface tension calculations using the exp-6 force field required a longer equilibration time. C. Molecular Orientation. The effect of the interface on the orientation of molecules at and away from it can be studied using the orientation order parameter
1 S(z) ) 〈3 cos2 θ - 1〉 2
Figure 2. Running average of the surface tension as a function of time for C9F20 at 323 K (left side) and for C20F42 at 450 K (right side) using the OPLS force field. Data averaged over blocks of 1 ps.
different definitions for the molecular orientation vector as described below. To represent the orientation of the backbone segments, the vector connecting a pair of carbon atoms that are two units apart in a molecule was used. According to this definition, positive values of S are indications for the preference of the molecular segments to orient perpendicular to the interface while negative values of S refer to preference for parallel orientation. Knowing the orientations of the CF2 groups and the CF3 end-groups of a perfluorinated molecule should supplement our understanding of the orientation of the molecule near the interface. The other advantage of calculating and reporting these two orientations is that they can be observed experimentally using infrared-visible sum frequency spectroscopy techniques.37 To determine the orientation of the CF3 end groups, we used the vector connecting the C atoms with the center of the equilateral triangle formed by the three F atoms in the tetrahedron. In this case, positive values of S indicate that the -CF3 groups tend to point perpendicular to the surface, which in turn is an indication of the preference for the molecule end segments to orient perpendicular to the surface. Negative values of S imply preference of the molecular segments to orient parallel to the surface. Similarly, the orientation vector for the -CF2 groups was defined to be the vector connecting the C atoms with the middle of the line formed by the two F atoms. In this case, positive values of S are directly correlated to the tendency of the molecule segments to orient parallel to the surface with the carbon-backbone plane normal to the surface. The implications of negative values of S in this case, however, should be interpreted in conjunction with the result of the backbone segment calculation described in the previous paragraph. If the tendency of the orientation of the backbone segments is to be normal to the surface, then S determined from -CF2 groups’ orientation should be negative. If the orientation of the chain segments is parallel to the surface, then negative values of S imply that the carbon-backbone plane of the chain segments is parallel to the surface. III. Results and Discussions
(10)
where θ is the angle between a vector that represents the molecular orientation and the unit vector normal to the interface, that is, zˆ. Here 〈‚ ‚ ‚〉 represents an ensemble average over all vectors within a specified slab in the z direction. The value of the order parameter is zero for a set of randomly oriented vectors, 1 for all vectors aligned perpendicular to the interface, and -0.5 for all vectors aligned parallel to the interface. To characterize molecular orientation of a system, we use three
A. Density Profiles. The density profile is one of the fundamental parameters for characterizing the interface. Here we have determined the influence of the cutoff for the nonbonded interaction on the behavior of the density profile. Figure 3 shows typical density profiles, F(z), of C7F16 melt from simulation using the OPLS force field at 298 and 348 K for different cutoff values. The density at a given value of z, F(z), is calculated by partitioning the simulation cell into bins along the z direction and adding the masses of atoms in each bin per partitioned volume. The density profiles shown in the figure
Surface Tension and Surface Orientation
Figure 3. Effect of cutoff distance on density profiles. Shown is sample density profiles of C7F16 at 298 K (top) and at 348 K (bottom) as a function of z, the depth in the film, obtained from simulations using the OPLS force field with different cutoff (rc) values. The experimental density at 298 K from Caco et al.40 is also shown in the top figure.
Figure 4. Density of C6F14 (circles), C7F16 (up triangles), C8F18 (diamonds), C9F20 (down triangles), and C10F22 (squares), as a function of temperature obtained from simulations using the OPLS force field (top) and exp-6 force field (bottom). Experimental data from Caco et al.40 are plotted as open symbols.
demonstrate that there are two symmetric liquid-vapor interfaces. The density profiles in all cases are fairly flat in the center of the film and smoothly decay to zero at the free surfaces. Similar dependence on cutoff is also observed using the exp-6 force field. Let us focus first on the 298 K case, Figure 3a, where density profile for different cutoff values was investigated. The change in the density profile from 10 to 12 Å cutoff is considerable; the width of the density profile decreased while the liquid density increased. Increasing the cutoff further to 14 Å did not change the behavior of the density profile considerably but resulted in a slight increase in the liquid density, which was then comparable to the experimental density within the error of the simulation. The 10 Å cutoff density profile shows higher vapor density and thus larger interfacial thickness compared to the larger cutoffs. The difference between cutoff 14 and 15 Å (not shown) is very minor, and increasing the cutoff to 16 Å
J. Phys. Chem. C, Vol. 112, No. 13, 2008 5033
Figure 5. Surface tension of short perfluorinated alkane chains as a function of temperature obtained from simulations using the OPLS force field (top) and exp-6 force field (bottom). Symbols are as indicated in Figure 4.
has no observable change. However, to minimize finite size effects on density and surface tension calculations we used a cutoff of 16 Å for the rest of the simulations. The liquid density of C7F16 at 298 K for different cutoff values is reported in Table 4. The behavior of the density profile also depends on temperature. In Figure 3, for a given cutoff value, as the temperature is increased from 298 to 348 K the interfacial thickness increased, and as expected the liquid density decreased. The influence of the cutoff on density profile is particularly crucial at higher temperatures as shown in Figure 3b. For the higher temperature, 348 K, the use of large cutoff value resulted in a better density profile and clearly justifies the use of large cutoff value, 16 Å for the rest of the simulations. For the range of temperatures and chain lengths we investigated, the shape of the liquid-vapor interface can be described well by a monotonically decreasing function of the form of an error function as was found for most of the liquid-vapor interfaces.33,38,39 Figure 4 shows average liquid densities for short perfluorinated alkanes as a function of temperature obtained from simulations using both force fields and rc ) 16 Å. The average densities were calculated from the bulk region of the density profiles, that is the region in the center of the film over which the density is fairly a constant. Available experimental density values are also shown in the figure. All of the liquid densities determined from our simulations using the two force fields and available experimental densities for the different systems we studied are compiled in Table 3. The standard error in the mean of the simulated density at each state point is less than 1% for short chains and less than 2% for longer chains. The experimental densities we could find in the literature are limited to short chain lengths. For the relatively long chain, C98F198, simulation results using only the exp-6 force field are reported because simulations using the OPLS force field for these systems resulted in a crystalline state for the range of temperatures studied. In general, the simulation results are in excellent agreement with the experimental results, with the exp-6 force field consistently giving better agreement with the experimental data. Both force fields slightly underpredict the densities and the maximum error in comparison with experiment is 1.8% for the OPLS force field and 1.2% for the exp-6 force field. The difference in density values from the two force fields is on
5034 J. Phys. Chem. C, Vol. 112, No. 13, 2008
Tsige and Grest
TABLE 3: Simulation Results for Liquid Densities of Perfluorinated Alkanesa C6F14
C7F16 Fopls
Fopls
Fexp6
Fexp
Fexp6 T/K F opls Fexp6
288 298 313 323 348
1.68 1.65 1.60 1.57 1.47
1.69 1.70 1.73 1.74 1.76 1.77 1.78 1.79 1.80 1.80 1.81 1.83 1.83 450 1.67 1.68 1.70 1.72 1.73 1.73 1.75 1.77 1.77 1.78 1.79 1.80 1.80 475 1.61 1.63 1.65 1.67 1.69 1.70 1.71 1.73 1.73 1.74 1.75 1.77 1.76 500 1.58 1.62 1.64 1.67 1.69 1.71 1.72 1.74 1.74 525 1.47 1.53 1.56 1.59 1.62 1.63 1.67 550
1.56 1.47 1.38 1.29
1.60 1.53 1.49 1.41 1.32
exp
F (g/cm3
γp (mN/m)
1.59 1.65 1.68 1.70 1.70
6.22 8.45 10.23 10.65 10.84
10 12 14 15 16
γtail (mN/m)
γp + γtail (mN/m)
3.21 2.71 2.19 1.98 1.64
Fexp
Fopls
550 565 575 600 625
C96F194 Fexp6
1.51 1.50 1.48 1.56 1.43 1.50 1.37 1.47
1.52 1.46
experimental data taken from Caco et al.40
TABLE 4: Simulation Results for the Surface Tension and Liquid Density of C7F16 Melt as a Function of Cutoff at 298 Ka rc (Å)
Fexp6
C48F98
Fopls
using exp-6 force field.
Fopls
C20F42
Fexp6 T/K
exp6
Fexp
C10F22
Fexp6
using OPLS force field.
Fexp6
C9F20
Fopls
a opls
Fexp
C8F18
T/K
9.42 11.15 12.43 12.63 12.48
a The experimental values are 1.73 g/cm3 and 12.53 mN/m for the liquid density and surface tension, respectively, taken from Caco et al.40
average 0.6% (which is within the error of the simulation). However, the discrepancy in density between the two force fields for the C48F98 system is large, as high as 6%. B. Surface Tension. Simulation results for the surface tension γp, without tail correction, and γtail, the corresponding tail correction calculated using eq 8, for C7F16 system at 298 K and for different cutoff values are reported in Table 4. The surface tension γp increases with rc while γtail, as expected, decreases with increasing rc. For cutoff values rc ) 10 and 12 Å the sum, γ ) γp + γtail, differs appreciably from those γ values obtained with larger cutoff values. The discrepancy in γ between the larger cutoff values rc ) 14, 15, and 16 Å is small, less than 2%, and is within the error of the simulation. For protonated alkanes we have shown that for cutoff values on the order of 16 Å the liquid density and surface tension, γ, agree very well with results using long-range Ewald summation for dispersion forces.41 We also confirmed this in our present study with a few selected points (not shown here). Figure 5 shows the surface tension (γ ) γp + γtail) of short perfluorinated alkanes as a function of temperature for both OPLS and exp-6 force fields with rc ) 16 Å. Available experimental surface tension results are also included in the figure. The simulation results show that the surface tension decreases with increasing temperature and for a given temperature it increases with increasing chain length, both in accordance with experiment. All of the surface tension values we calculated for different perfluorinated alkanes and their corresponding known experimental surface tension values are reported in Table 5. The standard error in the mean of the simulated surface tension at each state point is less than 4%. In general, the simulation results from both force fields are in good agreement with experiment. Specifically, for C6F14 and C7F16
the simulation results from both force fields are in excellent agreement with experiment. For chain lengths longer than C7F16 the exp-6 force field substantially overpredicts the surface tension while the OPLS force field slightly overpredicts it and gives better agreement with experiment. However, for C20F42 and C48F98 the OPLS force field predicts higher surface tension than the corresponding state point from the exp-6 force field. Because of the large melting temperature and high vapor pressure, measuring the surface tension of C20F42 experimentally for temperatures more than a few degrees above the melting temperature has proven to be very difficult.42 As a result, the available surface tension experimental data for C20F42 is limited to only one temperature and is about 10 mN/m at 443 K;42 the simulated surface tension value at this temperature using the exp-6 force field is 10.96 mN/m. This indicates that for long chains the exp-6 force field may do a better job in predicting surface tension values than the OPLS force field, which predicted relatively higher surface tension values. For protonated alkanes, we showed that the OPLS force field gives excellent agreement with experimental data for surface tension at low temperatures, while the exp-6 force field agrees extremely well at lower molecular weights. However, unlike the perfluorinated case, the exp-6 force field was found to substantially overpredict the surface tension for protonated alkanes longer than hexane. C. Molecular Orientation. Figure 6a-c shows the three different orientation parameters discussed in Section IIC for C10F22 at 298 and 348 K and for C20F42 at 475 K, respectively, from simulations using the OPLS force field. Because the films are symmetric about the center of the film, which is set at z ) 0 for the order parameter calculation, the orientation parameters shown in the figure are averages of the orientation parameters from the two sides of the film. In all cases, the orientation parameters consistently show that the chain segments tend to point perpendicular to the surface with the -CF3 groups exposed to the surface while the orientation is random in the middle of the film, which is out of the range of the influence of the surface. The C10F22 chains show significant surface ordering at 298 K where the ordering for z > 30 Å represents the vapor region of the film. The surface ordering shows a significant decrease in ordering when the temperature is increased from 298 to 348 K. Furthermore, increasing the chain length to C20F42 increased the chain ordering at the interface significantly. Note that the orientation of the C20F42 system shown in the figure is at 475
TABLE 5: Simulation Results for Surface Tension (γ ) γp + γtail) of Perfluorinated Alkanes for rc ) 16 Åa C6F14
C7F16 γopls
γopls
γexp6
γexp6
T/K
γopls
288 298 313 323 348
12.45 11.25 10.07 8.87 6.91
12.68 12.47 13.48 13.81 13.60 14.66 15.25 14.48 15.55 15.85 15.22 17.71 18.01 450 11.70 11.36 12.48 12.86 12.53 13.79 14.26 13.54 14.50 14.70 14.14 16.55 16.70 475 10.43 11.12 11.50 12.23 13.01 11.93 13.19 13.59 12.77 15.47 15.52 500 9.69 10.27 10.70 11.47 12.02 12.26 12.56 14.09 14.28 525 7.32 8.18 8.70 9.56 9.95 10.26 12.14 550
9.12 7.23 6.08 4.81
10.32 9.01 7.20 5.91 4.55
550 565 575 600 625
9.61 9.02 7.36 7.12 6.69 5.79 5.61
exp
γopls
γexp6
γexp
γopls
γexp6
C48F98
γopls
using exp-6 force field.
γexp
C20F42 T/K
exp6
γexp
C10F22
γexp6
using OPLS force field.
γexp6
C9F20
γopls
a opls
γexp
C8F18
T/K
experimental data taken from Caco et al.40
γexp6
C96F194 γexp6
9.51 8.78 6.18
Surface Tension and Surface Orientation
J. Phys. Chem. C, Vol. 112, No. 13, 2008 5035 Sciences user facility. Sandia National Laboratories is a multiprogram laboratory operated by Sandia Corporation, a Lockheed-Martin Company, for the U.S. Department of Energy under Contract No. DE-AC04-94AL85000. References and Notes
Figure 6. Profile of the order parameter for (a) C10F22 at 298 K, (b) C10F22 at 348 K, and (c) C20F42 at 475 K.
K, much higher than the temperatures used for the C10F22 orientation calculations. In general, the amount of surface orientation of the chains depends strongly on temperature, as expected, decreasing with increasing temperature, and also on chain length, where at a given temperature the ordering increases with increasing chain length. However, the ordering at the interface for all of the systems we studied is not that strong because the value of the order parameter we calculated is small (S < 0.35). IV. Conclusions The surface properties of various perfluorinated alkanes were studied by molecular dynamics simulations using two all-atom force fields. The computed liquid densities from both models are in excellent agreement with existing experimental data for short chains. Longer chains at low temperatures are better described by the exp-6 force field because the OPLS force field overpredicted the melting temperature of the chains substantially. For longer chains at high temperatures, the discrepancy in the liquid density between the two force fields is very large compared to the discrepancy for shorter chains. We have found that a longer simulation time of at least 7 ns was required in order to obtain reliable surface tension values. Both models yield good agreement with experimental data. Despite the fact that both models predicted liquid densities with sufficient accuracy, the exp-6 force field overpredicted the surface tension considerably with chain length (with discrepancy of 1.6-9%) while the OPLS force field overpredicted it slightly (with discrepancy of 0.2-4%) and thus gave better agreement with experiment. For longer chains, however, the OPLS force field predicted higher surface tension than the corresponding state point from the exp-6 force field. At the surface, the chain segments are found to point perpendicular to the surface with the -CF3 end groups at the surface. The amount of surface orientation of the chains strongly depends on temperature and also on chain length, where at a given temperature the ordering increased with increasing chain length. Acknowledgment. We thank Ahmed Ismail and Flint Pierce for a critical reading of the manuscript. This work was performed, in part, at the Center for Integrated Nanotechnologies, a U.S. Department of Energy, Office of Basic Energy
(1) Burger, L. L.; Cady, G. H. J. Am. Chem. Soc. 1951, 73, 4243. (2) Rothlistberger, U.; Laasonon, K.; Klein, M. L. J. Chem. Phys. 1996, 104, 3692. (3) Sprik, M.; Rothlistberger, U.; Klein, M. L. Mol. Phys. 1998, 97, 355. (4) Cui, S. T.; Cochran, H. D.; Cummings, P. T. J. Phys. Chem. 1999, 103, 4485. (5) McCabe, C.; Bedrov, D.; Smith, G. D.; Cummings, P. T. Ind. Eng. Chem. Res. 2001, 40, 473. (6) Watkins, E. K.; Jorgensen, W. L. J. Phys. Chem. A 2001, 105, 4118. (7) Friedmann, R.; Naumann, S.; Brickmann, J. Phys. Chem. Chem. Phys. 2001, 3, 4195. (8) Hoerauf, H.; Kobuch, K.; Dresp, J.; Menz, D. H. Lab. InVest. 2001, 239, 373. (9) Roider, J.; Hoerauf, H.; Lorbuch, K.; Menz, D. H. Clin. InVest. 2002, 240, 965. (10) Borodin, O.; Smith, G. D.; Bedrov, D. J. Phys. Chem. B 2002, 106, 9912. (11) Sadford, G. Tetrahedron 2003, 59, 437. (12) McCabe, C.; Bedrov, D.; Borodin, O.; Smith, G. D.; Cummings, P. T. Ind. Eng. Chem. Res. 2003, 42, 6956. (13) Pashkevich, D. S.; Shelopin, G. G.; Mukhortov, D. A.; Petrov, V. B.; Alekseev, Y. I.; Asovich, V. S. Russ. J. Appl. Chem. 2004, 77, 1847. (14) Tsuzuki, S.; Uchimary, T.; Mikami, M. J. Chem. Phys. 2004, 121, 9917. (15) Dias, A.; Caco, A. I.; Coutinho, J.; Santos, L.; Pineiro, M. M.; Vega, L. F.; Gomes, M. F. C.; Marucho, I. M. Fluid Phase Equilibra 2004, 225, 39. (16) Broniatowski, M.; Minones, J.; Macho, J.; Dynarowicz-Latka, P. Pol. J. Chem. 2005, 79, 1047. (17) Takai, K.; Takagi, T.; Baba, T.; Kanamori, T. J. Fluorine Chem. 2007, 128, 10. (18) Ameduri, B.; Boutevin, B. Well-Architectured Fluoropolymers: Synthesis, Properties and Applications; Elsevier Science: New York, 2004. (19) Sakka, T.; Ogata, Y. H. J. Fluorine Chem. 2005, 126, 371. (20) Borodin, O.; Smith, G. D.; Bedrov, D. J. Phys. Chem. B 2006, 110, 6279. (21) Shin, S.; Collazo, N.; Rice, S. A. J. Chem. Phys. 1992, 96, 1352. (22) Collazo, N.; Shin, S.; Rice, S. A. J. Chem. Phys. 1992, 96, 4735. (23) Shin, S.; Collazo, N.; Rice, S. A. J. Chem. Phys. 1993, 98, 3469. (24) Cui, S. T.; Siepmann, J. I.; Cochran, H. D.; Cummings, P. T. Fluid Phase Equilibria 1998, 146, 51. (25) Hariharan, A.; Harris, J. G. J. Chem. Phys. 1994, 101, 4156. (26) Holt, D. B.; Farmer, B. L.; Macturk, K. S.; Eby, R. K. Polymer 1996, 37, 1847. (27) Siepmann, J. J.; Karaborni, S.; Smit, B.; Klein, M. L. AIChE 1994 Spring Meeting, paper 94c, 1994. (28) Byutner, O. G.; Smith, G. D. Macromolecules 2000, 33, 4264. (29) Martin, M. G.; Siepmann, J. J. J. Phys. Chem. B 1999, 103, 4508. (30) Plimpton, S. J. J. Comput. Phys. 1995, 117, 1. (31) Hockney, R.; Eastwood, J. Simulation using Particles; Adam Hilger: New York, 1988. (32) Kirkwood, J. G.; Buff, F. P. J. Chem. Phys. 1949, 17, 338. (33) Ismail, A. E.; Grest, G. S.; Stevens, M. J. J. Chem. Phys. 2006, 125, 14702. (34) Chapela, G. A.; Savilles, G.; Thompson, S. M.; Rowlinson, J. Chem. Soc., Faraday Trans. 1977, 73, 1133. (35) Blokhuis, E. M.; Bedeaux, D.; Holcomb, C. D.; Zollweg, J. A. Mol. Phys. 1995, 85, 665. (36) Lacasse, M. D. Grest, G. S.; Levine, A. J. Phys. ReV. Lett. 1998, 80, 309. (37) Yurdumakan, B.; Harp, G. P.; Tsige, M.; Dhinojwala, A. Langmuir 2005, 21, 10316. (38) Huang, J. S.; Webb, W. W. J. Chem. Phys. 1969, 50, 3677. (39) Beysens, D.; Robert, M. J. Chem. Phys. 1987, 87, 3056. (40) Caco, A. I.; Dias, A. A.; Marrucho, I. M.; Pineiro, M. M.; Santos, L. B.; Coutinho, J. P Fifteenth Symposium on Thermophysical Properties 2003, 15, 134. (41) Ismail, A. E.; Tsige, M.; Veld, P. J. in ’t; Grest, G. S. Mol. Phys. 2007, 105, 3155. (42) Sauer, B. B.; Dee, G. T. Macrmolecules 1994, 27, 6112.
