14520
J. Phys. Chem. C 2010, 114, 14520–14527
Computational Investigation of the Influence of Surfactants on the Air-Water Interfacial Behavior of Polycylic Aromatic Hydrocarbons Collin D. Wick,*,† Bin Chen,‡ and Kalliat T. Valsaraj‡ Louisiana Tech UniVersity, Ruston, LA, and Louisiana State UniVersity, Baton Rouge, LA ReceiVed: April 30, 2010; ReVised Manuscript ReceiVed: June 28, 2010
A combination of Monte Carlo and molecular dynamics simulations was carried out to investigate the effect of 1-octanol surface coverage on the interfacial partitioning and behavior of polyaromatic hydrocarbons. Also, how the surface coverage of 1-octanol is related to its gas-phase density was examined. It was found that 1-octanol partitioned sparsely, preferring to lie flat at the water surface at low gas-phase concentrations. As the gas-phase concentration increased, the 1-octanol surface coverage increased rapidly, and it oriented more perpendicular to the water surface. The interfacial partitioning of polyaromatic hydrocarbons was enhanced significantly for larger ones, such as naphthalene and anthracene, when 1-octanol was present at the air-water interface, but benzene had a similar portioning for both cases. Naphthalene preferred to lie flat on the bare water surface, but with 1-octanol present, it oriented fairly perpendicular to the surface when close to the water surface and parallel to the surface when at the interface of 1-octanol and air. I. Introduction Water droplets or hydrometeors are ubiquitous in the atmosphere as fog, mist, clouds, and other forms of atmospheric moisture. These droplets can serve as mediums for the chemical transformation of atmospheric trace gases that have important environmental consequences, such as the creation of smog components, water acidification, and carcinogen transport and creation. Many trace gases readily solvate into these droplets, such as formaldehyde, but some molecules may be found at the droplet surface. Some surface-active gas species include ozone, oxygen, and hydroxyl radicals.1–4 For these species, the air-water interface may be a significant factor in their reactivity. This becomes more important as the droplets’ size decreases, as the surface-to-volume ratio becomes larger. For instance, an interfacial reaction mechanism has been proposed for the creation of chlorine from hydroxyl radicals while in contact with very small saltwater aerosols in the marine boundary layer.5 Hydrophobic species, such as many pesticides, herbicides, and polyaromatic hydrocarbons (PAHs), though, are known to strongly accumulate at the air-water interface, and surface effects are dominant.6–10 PAHs are common in the atmosphere, as there are a variety of sources, including organism byproducts and the combustion of an array of organic matter.11,12 Many PAHs are suspected carcinogens, and their oxidation products can be even more toxic than the parent compound. Experimental and computational investigations have found that, for a variety of PAHs, they have an interfacial free energy that is much greater than in the bulk,13–17 and their oxidation into potentially environmentally harmful species occurs readily at the air-water interface.18 These investigations have brought significant insight into the behavior of PAHs at the surface of distilled water. However, realistic hydrometeors are not made out of distilled water but contain many other components. * To whom correspondence should be addressed. E-mail: cwick@ latech.edu. † Louisiana Tech University. ‡ Louisiana State University.
Field measurements in the gulf coast region have found the presence of a variety of surfactants, including humic and fulvic acid-like species, present in fog and rain droplets.19 Surfactant species are known to influence the properties of droplets,14,20–23 and, specifically, the uptake and reactivity of neutral matter. For instance, the ozonolysis rate of some PAHs on water droplets and surfaces has been found to be enhanced by the presence of certain surfactants.18,24 On the other hand, some surfactant species, generally of lower molecular weight, have been found to reduce the uptake of different PAHs.25 There is a need to understand the factors that promote different adsorption behaviors that result with the presence of surfactants at the surface of droplets. Experimental measurements of the interfacial partitioning of trace gases at the air-water interface with surfactants would be of significant benefit for the investigation of the influence of surfactants. Unfortunately, to the knowledge of the authors, no such studies exist as of yet. Computational methods can be a very useful tool to calculate interfacial free energy differences between multitudes of states and have the added benefit of providing additional molecular level insight in an unambiguous manner. However, computational methods of the scale needed to calculate these properties generally rely on a molecular model that has to be parametrized. As a result, one needs to be careful that any computational study uses a molecular model in which its strengths and weaknesses are understood, and direct experimental verification is desired. Molecular dynamics (MD) simulations have been used to understand PAH partitioning at the air-water interface and have been shown to provide good agreement with experiment for their interfacial partitioning.15,16 The authors of this referenced work found that interfacial benzene molecules preferred to orient with their aromatic plane parallel to the water surface. With the addition of surfactants, interfacial partitioning and PAH orientation are likely to be different than at the neat air-water interface. The addition of surfactants adds an additional challenge to MD simulations in how one determines the surface coverage of the surfactant. Monte Carlo (MC) simulations in the open ensembles, Gibbs or grand canonical, can be used to determine the surface coverage of a surfactant under saturation
10.1021/jp1039578 2010 American Chemical Society Published on Web 08/02/2010
Computational Investigation of the Air-H2O Interface or other vapor pressure conditions.26 With these techniques, and quality molecular models, a new insight of how PAH molecules behave at the surface of droplets with surfactants present can be developed, which will improve the understanding of atmospheric processes. While the atmospheric droplets found in the gulf coast region have a very complex assortment of surfactant species,19 1-octanol (OCT) can be used as a model surfactant that has a moderate molecular weight. This paper reports work on a computational investigation using both MD and MC simulations methods to elucidate the effect of OCT on the interfacial partitioning, structure, and orientation of droplets with PAHs present. II. Molecular Models As stated earlier, the surfactant investigated in this work was OCT, in which the TraPPE-UA force field was used.27 The TraPPE-UA model utilizes pseudoatoms for all alkyl groups located at the carbon centers (i.e., no alkyl hydrogens are present), and all other atoms are treated explicitly, including the alcohol oxygen and hydrogen. For the PAH molecules, the TraPPE-EH model was used,28 which includes atomic sites for all atoms. The TraPPE force field, in general, is parametrized to reproduce the vapor-liquid coexistence curves for a variety of molecules. Because of this, most of the TraPPE molecular models exhibit good agreement with experimental densities and thermodynamic properties over a fairly wide range of conditions, often doing a good job of predicting multicomponent phase equilibria.29–31 In addition, the MD simulations did not use a rigid model, as is the case for the TraPPE-EH model,28 due to challenges in simulating large rigid models with MD. Instead, the MD simulations used bonded parameters extracted from the Amber force field,32 except that the MD simulations used the bond lengths prescribed in the TraPPE-EH models and fixed carbon-hydrogen bonds. The water model used for this investigation was the TIP4P model,33 which has been used, along with the TraPPE model, previously to predict partitioning between water and OCT.34 Furthermore, these two models have been used together to study the partitioning of 1-butanol in bulk water and at the air-water interface.26 It should be noted that other models have been used to understand PAH portioning at the air-water interface,15 but the TraPPE model along with the TIP4P model has a history of providing reasonable experimental agreement for thermodynamic properties.26,34 For the modeling of the PAH molecules in aqueous environments, the TraPPE-EH force field had to be modified. The force field was parametrized in a relatively nonpolar environment,28 which may not be transferable to the more polar environment of water and its interface. As a result, the charges were increased from 0.095 to 0.13 for the H atoms to reproduce the experimental free energy of solvation in water. All C atoms that were bonded to H atoms had charges of -0.13, and for all tertiary bonded C atoms (those only bonded to other C’s), their charge was neutral. This allowed a very transferable model to be used for the PAH molecules, in which each PAH molecule had only three possible interaction sites no matter what the PAH structure was. The specific PAH molecules investigated include benzene (BEN), naphthalene (NAP), and anthracene (ANT). III. Simulation Details A. Monte Carlo Simulations. Coupled-decoupled configurational-bias Monte Carlo (CBMC) simulations35–38 in the grand canonical ensemble were used to obtain the equilibrium surface coverage of OCT and also for the interfacial partitioning of the PAH molecules at 300 K. Similar to the Gibbs ensemble setup
J. Phys. Chem. C, Vol. 114, No. 34, 2010 14521 used previously for simulating the 1-butanol coated surfaces,26 these simulations were carried out using separate, but thermodynamically connected, boxes (to simulate the different phases involved) except that the vapor box is represented by an idealgas particle reservoir. As one part of the simulation conditions, the chemical potential of the reservoir was specified, but to represent an idea gas, the particles in this system did not interact. Using this setup, the chemical potential of different species can be related to the ideal gas pressure or density.39 Before running simulations in the grand canonical ensemble, the liquid water was equilibrated in a box in the NpT ensemble. The simulation box was then elongated in the z axis to create two vapor-liquid interfaces with this water liquid slab occupying one third of the box. Once the interfaces were created, the volume of this simulation box and the interfacial area were held constant, and further equilibration in the NVT ensemble was carried out. Five types of MC moves were employed to sample the phase space: translation, rotation, CBMC conformational, CBMC swap, and CBMC switch moves.34,37,40 The first three were used to ensure thermal equilibrium, whereas the last two were performed to equilibrate the chemical potential with the ideal gas reservoir. In particular, the swap move involves deleting one particle from one phase and regrowing it with CBMC in another phase to allow for the direct equalization of the chemical potential for this species between the two phases. In addition, the switch move was used to exchange two molecules that shared a similar structure.40 For instance, BEN can be exchanged with a large part of NAP, and NAP with ANT, requiring only the growth of four more carbon atoms and hydrogen atoms for the exchange. This gives a significantly higher acceptance rate than the swap move, which has to regrow the whole molecule, especially for larger solutes, such as ANT. The x and y dimensions of the system, which were parallel to the interfaces, were both near 64 Å. A total of 6912 water molecules were used, and they were all kept in the liquid box (i.e., no swap/switch moves were employed for water) to yield an average thickness along the z dimension of ∼40 Å for this liquid slab. The nonbonded interactions were truncated at 28 Å. No tail corrections were included for the Lennard-Jones (LJ) interactions, and the Ewald sum was employed to deal with the long-range electrostatic interactions.41 A range of OCT ideal gas densities were specified, including the following values: 0, 2 × 10-10, 5 × 10-10, 1 × 10-9, 2 × 10-9, 5 × 10-9, and 1 × 10-8 molecules per Å3. This allows the creation of surfaces of different OCT coverages (i.e., from bare water surfaces to the fully OCT coated ones). The partitioning of the PAH molecules was investigated for only two of the surface coverages: the bare water surface and a surface saturated by OCT, which was obtained at an OCT ideal gas-phase density of 5 × 10-9 molecules per Å3. To investigate the partitioning of the PAH solutes, their ideal gas densities were adjusted so that, on average, greater than 0.1 molecules were found in the liquid box at any one time. Given that these solutes tend to deposit at the interface (see the Results and Discussion section), to a good approximation, it can be assumed that all of the solute molecules found in the liquid box are located at the interface. Therefore, the average number of solute molecules found in the liquid box was used to compute its surface concentration and, correspondingly, the air-interface partition coefficient, KIA, which can be subsequently converted to the free energy transfer from the air to the interface as follows: ∆GIA ) -RT ln(KIA/δ0), where δ0 refers to the standard state surface thickness and typically a value of 6 × 10-10 m was used.14,17
14522
J. Phys. Chem. C, Vol. 114, No. 34, 2010
Additional grand canonical Monte Carlo simulations were carried out to evaluate the solvation free energy of the PAH molecules in bulk water and in OCT. These simulations used a pure liquid water or OCT box without an interface coupled to an ideal gas particle reservoir. Because the long-range interactions can be conveniently included for both LJ (through tail corrections) and electrostatic interactions (via Ewald sum), a relatively smaller cutoff of 14 Å and a smaller system size containing 864 water molecules or 240 OCT molecules were used. However, compared to the interfacial simulations described above, particle insertions in the interior of the bulk liquid water or OCT are much more difficult to sample than at the interfacial region. To alleviate this problem, many intermediate structures were introduced in an expanded ensemble approach.42 For example, BEN was regrown inside the liquid phases by adding one C-H group at a time and the free energy was integrated along this growth trajectory. Also, the charges were initially turned off for these intermediates until BEN was fully grown. The charges were then gradually turned on from 0 to 0.13 e, and many intermediate charge values were used for a smooth evaluation of the free energy changes during this charge evolution. To improve the statistics for calculating the free energy of particle insertion, umbrella sampling was used,43 and the umbrella biasing potential was adjusted so that all the intermediate species involved had roughly an equal probability to appear in the bulk liquid water or OCT. Umbrella sampling was also used to calculate free energy profiles for the PAH molecular transfer across the interfaces.43 For these simulations, a biasing potential corresponding to the distance between the PAH center of mass and liquid water center of mass along the z axis (denoted zPAH) was used. This potential had no functional form for the MC simulations and was adjusted so that all species had roughly equal probabilities to be at any position within a certain range of zPAH. B. Molecular Dynamics Simulations. The coordinates from the MC simulations were used as an input for the MD simulations, but the MD system sizes were reduced to approximately 25% of their original size, and one PAH molecule was used in the simulations. The system used in the MD simulations had half the x and y dimensions of the MC simulations, giving lengths of 32 Å, but the same z dimensions as the MC simulations of around 40 Å were used, except where noted. The actual number of water and OCT molecules varied for each system, but there were between 1650 and 1704 water molecules and 85 and 90 OCT molecules present. To calculate the free energy profile of PAH molecules across the water surface, both a constrained approach and umbrella sampling were used. The constrained approach maps out a free energy profile across the surface in a manner similar to what has been used previously to study ions at the air-water interface.44 This method fixes the distance between the z coordinate (which is normal to the interface) of the PAH center of mass and the liquid center of mass (zPAH) and calculates the force acting on the molecular center of mass. This was carried out for multiple simulations in 1 Å increments. Integrating this force over the z direction gives the potential of mean force.44 This method was found to give a good (i.e., faster converging with lower uncertainty) free energy profile in the low-density regions of the system, but umbrella sampling43 was found to be more efficient when the solute was in contact with the water or OCT molecules. The reason for this was due to the very large changes in free energies over a small range of positions across the interface, which required many more points in the integration, making the constrained approach prohibitively expensive. As a
Wick et al. result, umbrella sampling was used to calculate the free energy profile when the PAH molecule was in contact with the OCT or water liquids. The umbrella sampling simulations utilized a harmonic biasing potential with a force constant of 1.0 kcal/mol. Between 30 and 40 different simulations were carried out for each system, each centered at a different zPAH value ranging from 10.0 to 40.0 Å for the neat air-water interface, and between 10.0 and 50.0 Å for the surfactant-coated system, all in 1.0 Å increments. Both the umbrella sampling and the constrained approaches were carried out for 1 ns for each simulation window, with at least 200 ps of equilibration preceding each production run. Equilibration was deemed to be reached when subsequent 200 ps simulation runs gave the same results. For instance, the ANT simulations with the surfactant-coated surface had 600 ps of equilibration preceding the production for the umbrella sampling simulations. The two methods were matched in overlapping regions to give one continuous curve, in which the results are discussed later. The free energy of solvation from the MD simulations was simply taken from the differences between the gas and the liquid free energy. However, to calculate the interfacial free energy is not as straightforward. To make a valid comparison with experimental data requires comparing the partition coefficient calculated from the surfactant excess (not the free energy minimum) with the experimental results. The methodology to calculate the surface excess has been described by other workers previously,45 and we use it as described in the reference. The surface excess can be used to define a partition coefficient (KIA), in which an interfacial free energy can be defined the same as described for the MC simulations (∆G ) -RT ln[K/6.0 Å]).17 The MD simulations were carried out in the NVT ensemble with the Berendsen thermostat to keep the temperature at 298 K,46 and a time step of 1 fs was used for all MD simulations. The water molecular geometry and PAH carbon-hydrogen bond lengths were kept fixed with the SHAKE and RATTLE algorithms.47 A potential truncation of 12 Å was enforced for the LJ interactions, and the long-ranged electrostatics were handled with the particle mesh Ewald summation technique.48 It should be noted that all free energy profiles had an extra term in them to account for long-ranged LJ interactions beyond the cutoff (as tail corrections would not account for this in an asymmetric system). Using the average water density from the simulations, the interaction of the PAH atoms with water at this density beyond 12 Å was added to the solvation free energy, and half of this energy was added to the interfacial free energy. For anthracene, this correction was less than 2 kJ/mol for the solvation free energy and less than that for the others. All MD simulations were carried out with the Amber 10 simulation package49 with a minor modification to implement umbrella sampling and the constrained method. IV. Results and Discussion A. OCT Partitioning. From the MC simulations in the grand canonical ensemble, the average number of OCT molecules () that partitioned to each water surface (as there are two) and the average area of surfactant headgroups () as a function of ideal gas density is given in Table 1. It can be observed that, at the two lowest ideal gas densities, there is little interfacial partitioning, but between 0.5 × 10-9 and 1 × 10-9 Å-3, increases very rapidly, about an order of magnitude. At higher ideal gas densities, the increase in is moderate and levels off to around 170. The area per surfactant group can be compared with experimental values of 26.4 Å2 at
Computational Investigation of the Air-H2O Interface
J. Phys. Chem. C, Vol. 114, No. 34, 2010 14523
TABLE 1: Average Number of OCT Molecules () at Each Air-Water Surface and the Area per Octanol Group as a Function of Ideal Gas Number Densitya Fideal (Å-3)
(Å2)
0 2 × 10-10 5 × 10-10 1 × 10-9 2 × 10-9 5 × 10-9 1 × 10-8
0 2.6 ( 0.1 8.85 ( 0.2 111.8 ( 0.5 140.4 ( 0.3 158.2 ( 0.2 168.9 ( 0.1
0 1576 462 36 29.2 25.9 24.26
a The ideal gas number density can be multiplied by 2.16 × 102 to give a specific density (in g/cm3).
Figure 1. Density profile of water and OCT in the pure water system, a system with an OCT ideal gas density of 5 × 10-10 Å-3 (low coverage), and the OCT system with an ideal gas density of 5 × 10-9 Å-3 (high coverage). Solid lines represent water, dashed lines alkyl groups, and the dotted lines OCT oxygens. The low coverage alkyl and OCT oxygen densities are increased by a factor of 10 for clarity, and 0 represents the GDS of water.
293 K50 and 27.6 Å-3 at 288 K.51 The saturated vapor density of our 1-octanol model is 3.75 × 10-9 Å-3,27 which would give us results around 27 Å-3 at saturation, very close to experiment. The density profiles for water, OCT alkyl groups, and oxygen atoms are given in Figure 1 for three of the OCT ideal gas densities used. Zero in these plots represents the Gibbs dividing surface (GDS) of water, negative values represent toward the water bulk, and positive values toward the vapor. It should be noted that, at the low coverage, the OCT densities were increased by a factor of 10 for clarity. It can be observed that, with the addition of surfactant, the water interface becomes broader, which is expected, due to the reduction in surface tension when surfactants are added to a system.52 In addition, the OCT oxygens are all present in regions with water density
Figure 3. Free energy profiles for the BEN, NAP, and ANT across the air-water interface for MC (dashed lines) and MD (solid lines) simulations. Zero represents the GDS of water.
present, except at the high coverage, in which there is a small amount that is outside this region (as the green dotted line stretches slightly past 5 Å). There is a clear difference in the behavior of the OCT alkyl groups between the low and high coverage systems. In the low coverage system, the alkyl density distribution is centered much closer toward the GDS than in the high coverage system. This is likely due to the fact that, at this low coverage, there is an order of magnitude fewer OCT molecules present, and as a result, interactions with other OCT molecules will be minimal. The only interaction of the OCT alkyl carbons, then, will be with the water surface, resulting in the OCT molecules lying flat to maximize interactions with the water surface. When the OCT coverage is much higher, OCT molecules can readily interact with other OCT molecules, causing them to orient in a fashion that maximizes interactions with other OCT alkyl groups. Also, this orientation contributes to a higher surface concentration, minimizing each OCT molecule’s individual surface coverage. Snapshots of the low and high concentration OCT systems are given in Figure 2. These snapshots show that, at low coverages, the OCT molecules can lie flat, but at high coverages, they stand up, but not necessarily in a highly ordered fashion. B. Free Energy Profiles. The free energy profiles from the MC and MD simulations for the bare air-water interface are given in Figure 3. The free energy profiles start to become negative when they come within 15 Å of the GDS of water. When the molecules approach around 2.5 Å from the GDS, they all reach a free energy minimum at pretty much the same position. It is interesting to note that this is despite the fact that the PAH molecules have significantly different geometries, with ANT being much larger than BEN. How-
Figure 2. Snapshots of the low and high coverage systems described in Figure 1.
14524
J. Phys. Chem. C, Vol. 114, No. 34, 2010
Wick et al.
TABLE 2: Free Energies of Solvation and Interfacial Free Energies at the Bare and Surfactant-Coated Surfaces Calculated from MD and MC Simulations and from Experimenta MD
MC
experiment
molecule
∆solvG
bare
surf.
∆solvG
bare
surf.
∆solvG61
bare17
benzene naphthalene anthracene
-4.3 -10.6 -15.5
-13.2 -20.8 -30.7
-15.8 -29.6 -54.7
-4.2 -10.5 -16.1
-13.7 -21.1 -28.5
-16.2 -28.2 -40.2
-3.7 -9.5 to -10.1 -14.7 to -16.8
-16.3 -26.5 N/A
a
All values are in kJ/mol.
ever, if one takes into account that they are all planar and if they all had preferred orientations with their ring planes being parallel to the surface, as has been argued previously,7 then this result would be expected. The orientations will be described in greater detail later. The MC free energy profiles generally agree with those extracted from MD, with some differences for NAP and ANT near their free energy minimums. In particular, the ANT results have a lower free energy minimum from the MC simulations than the MD simulations. Some of the discrepancies for the larger ANT molecules may be the result of the fact that the MD simulations utilized ANT molecules that had a small degree of flexibility. Furthermore, the MC results were found to be more difficult to converge for the free energy profiles when the large PAHs were in contact with water. The free energy of solvation from the MD simulations can be determined by taking the difference between the free energy when the molecules are far away from the GDS with those when they reach a steady value in the bulk water. These were used to calculate the free energy of solvation from the MD simulations. As said in section III.A, the free energies for the MC simulations were extracted from expanded ensemble simulations in the grand canonical ensemble.42 Also, the interfacial free energies were calculated from the MD profiles and separate MC simulations as described in the simulation details. The process for calculating the MD interfacial free energies is given in section III.B. Both the solvation and the interfacial free energies calculated by MD and MC are given in Table 2. It should be noted that the solvation free energies for MD were only extracted from the free energy profiles for the bare water surface. The free energy of solvation agrees very well between simulation and experiment for all of the PAHs. This is despite the fact that the charge distribution is similar for all three PAHs and shows that the strategy employed for parametrizing the PAHs allows for a good degree of transferability. In addition, both the free energies of solvation and the interfacial free energies for the bare surface have a close agreement between the MD and MC simulation results, showing a good consistency between the two. For BEN and NAP, the interfacial free energy can be compared between simulation and experiment. The simulation results underpredict the interfacial free energy by a few kilojoules per mole. This was somewhat unexpected, given the good agreement for the free energy of solvation. The reason for this disagreement may be in the molecular model itself. A PAH molecule at an air-water interface would likely have a significant induced dipole as a result of the broken symmetry of the interface. This has been observed in perpetuity for large halide anions at the air-water interface.53–56 Because the molecular models used in this study are not polarizable, this effect is not included, which may account for the discrepancy between simulation and experiment. Nonetheless, this should not affect the qualitative behavior of the PAHs at the air-water interface, but developing polarizable models to address these is ongoing. The free energy profile across the surfactant-coated surface is given in Figure 4, and the corresponding interfacial free
Figure 4. Free energy profiles for the BEN, NAP, and ANT across the surfactant-coated interface for MC (dashed lines) and MD (solid lines) simulations. The dotted green line is from an additional MD simulation described in the text.
energies are given in Table 2. The presence of OCT increases the magnitude of the interfacial free energy strongly for NAP and ANT but has relatively little effect on BEN. Furthermore, the free energy minimum positions are shifted farther away from the GDS than at the neat air-water interface, from around 2.5 to 5.0 Å. As a result, it would be expected that the presence of OCT would have a minor effect on the accommodation of BEN, but a significant one for the accommodation of NAP and ANT. There is experimental evidence showing that the presence of OCT at the surface of water significantly increases the uptake of ANT in comparison with neat water.57 For BEN, to our knowledge, there are no experimental results available as to the effect of surfactants on its uptake. The MC and MD simulation results agree very well for both BEN and NAP, but there is a significant discrepancy for ANT. The dotted line in Figure 4 represents the original free energy profile for ANT at the surfactant-coated surface for the system size described in the simulation details (1650 water molecules and 85 OCT molecules). It can be observed that it agrees very well with the MC free energy profile. However, as the free energy profile approaches the middle of the water liquid, it plateaus near 0 kJ/mol. For NAP and for BEN, the free energy of solvation from the MD free energy profiles agreed very well with those calculated from the MC bulk simulations. Also, the free energy of solvation for all PAH molecules at the air-water interface agrees well with the bulk MC simulations. As a result, these all appeared to reach a bulklike state before being 10 Å from the GDS. The ANT system in the surfactant-coated system was the only situation where this was not the case. To check if this was a system size effect, the number of water molecules in the MD simulations was increased from 1008 to 2668, elongating the water liquid slab by approximately 40 Å to around 80 Å, with the new total box length being 140 Å, but with the same dimensions in the x and y directions (34 Å). With the elongated water slab, the free energy of solvation for ANT from the MD simulations (-14.0 kJ/mol) approaches the value calculated from the MC simulations in bulk water
Computational Investigation of the Air-H2O Interface TABLE 3: Comparison of the Free Energies for the Transfer of PAHs from Water to OCTa
a
molecule
simulation
experiment59
benzene naphthalene anthracene
-13.6 -21.4 -29.9
-12.2 -19.1 -25.7
J. Phys. Chem. C, Vol. 114, No. 34, 2010 14525 SCHEME 1: Vectors Used To Calculate the Angular Distribution for NAPa
All values are in kJ/mol.
(-16.1 kJ/mol). They do not have perfect agreement but are reasonably close. Of interest is that the free energy profiles for the two system sizes diverge as they approach 10 Å from the GDS toward the air, until around 5 Å from the GDS. The reason we found for this was that, as the large ANT molecule came in contact with the surface, it significantly perturbed the interfacial structure that appeared to percolate through the entire smaller system with the surfactant-coated surface. It appears to be a result of the combination of a large hydrophobic molecule solvating in a polar solvent with a reduced surface tension, but this is speculative. It should be noted that we found system size effects for interfacial solvation of different solutes into squalane.58 To provide confidence in our results for the effect of OCT on the interfacial partitioning of PAHs, we calculated the PAH Gibbs free energy of transfer from water to OCT to compare with experiment. The experimental values were taken from converting their KOW59 values (which were the octanol-water partition coefficient) to ∆G by ∆G ) -RT ln KOW. The simulation values were calculated by employing a thermodynamic cycle, by subtracting the gas-OCT free energy from the gas-water one. Although the effect of mutual solubilities of these two solvents (i.e., presence of water in octanol or octanol in water) was neglected by the simulation, for hydrophobic solutes, the presence of the other solvent has a negligible influence on the KOW.34,59,60 Table 3 gives our simulation results, along with results from experiment. Our free energies are overestimated in magnitude by around 1 kJ/mol per aromatic ring. This 1 kJ/mol per aromatic ring difference reflects mostly the deviation on the gas-OCT partitioning leg of the thermodynamic cycle, given the good agreement observed on the aqueous solvation free energy results. Interestingly, it was found to almost match with the amount of decrease on the free energy of transfer to the OCT phase when adding the 0.13 charges to the PAH molecules (i.e., the free energy for the expanded ensemble42 step from neutral PAH to charged). This means that the simulation results on the gas-OCT partitioning can be improved if the original lower-charged (or a neutral) TraPPE-EH PAH model was used, but this comes with a sacrifice on the gas/water partitioning results. A polarizable model would be desired for such an investigation. However, we are reasonably close to experiment to have confidence in the trends in interfacial partitioning we observe with the inclusion of OCT on the surface, but future experiments on these would be useful. C. Orientational Profiles. The orientational profiles were extracted from the MC free energy profiles to give the angular distribution between the xy plane (which is normal to the interface) and the vector forming the moments of inertia for NAP. The three moments for NAP are shown in Scheme 1 and are denoted as follows. NORM is the vector normal to the NAP plane, LONG is the vector from one end of NAP to the other (bisecting three carbon-carbon bonds), and SHORT is the vector normal to the other two. For the distributions, a value of cos θ ) 0 represents the specified vector as perpendicular to the interface, and cos θ ) 1.0 is parallel to the interface.
a
The dashed line is coming out of the plane.
Figure 5. Orientational profile for NAP at its free energy minimum at the neat air-water interface.
Figure 5 gives the orientational profile for NAP at the bare air-water interface while located at its free energy minimum (2.5 Å from the GDS). For the NORM direction, there is a very strong preference for it to be perpendicular to the interface, which represents the NAP having a preferential orientation with its plane parallel to the surface. It should be noted that a similar behavior was found for ANT and BEN, but we only showed NAP for clarity. In addition, LONG and SHORT prefer to be parallel to the interface, with the LONG direction showing the strongest preference. Clearly, NAP and the other PAHs prefer orientations in which contact with the air-water interface is maximized. This also has been described previously for PAHs at the bare air-water interface.7 Figure 6 gives the orientational distribution for the three vectors described at the surfactant-coated surface at the free energy minimum (5.0 Å from the GDS). The distributions shown here vary considerably from those at the bare air-water surface, showing an opposite behavior. The preferred orientation of the LONG and SHORT vectors is perpendicular to the air-water interface, and the NORM vector strongly prefers to be parallel to the interface. Interestingly, the LONG direction is more strongly oriented perpendicular to the surface, showing a likely geometry in which the SHORT axis of NAP is oriented perpendicular to the interface. As a result, NAP will not maximize contact with the water surface, as it did for the bare system. This is not unexpected, as for NAPH to maximize contact with the water surface, it would have to displace OCT molecules, which have stronger hydrogen bonds with water. In
14526
J. Phys. Chem. C, Vol. 114, No. 34, 2010
Figure 6. Orientational profile for NAP at its free energy minimum at the surfactant-coated surface.
Wick et al. For instance, the NORM distribution reaches a maximum of around 1.6 in Figure 7, whereas at the bare surface (Figure 5), the same orientational distribution reaches a maximum of 2.25, showing a much stronger preference for NORM to be perpendicular to the surface. Figure 8 shows a series of snapshots for NAP at the bare water and surfactant-coated surfaces. The water and surfactant species were taken from the snapshot where NAP was at the free energy minimum. The snapshots illustrate what was described for the orientational profiles. The NAP molecules lie down flat on the bare water surface, whereas NAP molecules prefer to be oriented with the LONG axis perpendicular to the surface in the surfactant-coated system. Also, it can be observed that, when NAP is present at the interface between the surfactant and air, its preferred orientation is to be relatively flat to the surface. V. Conclusions
Figure 7. Orientational profile for NAP at the surfactant-coated surface 10 Å toward the air from its free energy minimum.
addition, NAP interactions with the OCT molecules should already be fairly strong, giving a smaller energetic benefit for NAP to maximize interactions with the water surface. Again, the same behavior was found to be present for BEN and ANT in the OCT system but was not shown for clarity. Figure 7 gives the orientational distribution of NAP when it is 15 Å from the GDS of water, or 10 Å toward the air from its free energy minimum. It can be observed that the orientational behavior of NAP is markedly different than when it is at the free energy minimum. The NORM direction is now oriented perpendicular to the surface, and the LONG and SHORT directions are oriented parallel to the surface. This shares many similarities with what is observed at the bare water surface, except that the degree of orientational preference is not as strong.
A combination of MC and MD simulations was carried out to investigate the behavior of PAH molecules at bare and OCTcoated water surfaces. The presence of OCT enhanced the interfacial partitioning of NAP and ANT but only had a minor effect on the partitioning of BEN. OCT also had a noticeable effect on the interfacial orientation of NAP, in which NAP went from a preferred orientation that was parallel to the air-water interface at the bare water surface to orienting perpendicular to the air-water interface when OCT was present. Furthermore, when NAP was at the interface between OCT and air, NAP then preferred to orient parallel to the interface, but not to the degree that was found at the air-water interface. This work begins to bring insight into how the presence of surfactants will influence the interfacial partitioning of PAHs. Acknowledgment. Part of this research was funded by the Louisiana Board of Regents Research Competitiveness Subprogram, Contract No. 3LEQSF(2008-11)-RD-A-21, and part of it from the National Science Foundation (CHE/MCB0448918). The calculations were carried out using the resources from the Louisiana Optical Network Initiative (LONI), the Center for Computation and Technology, and the Office of Computing Services at LSU.
Figure 8. Snapshots of NAP at different regions at the bare surface and at the surfactant-coated surface. Configurations in which NAP is not in its free energy minimum show it as transparent.
Computational Investigation of the Air-H2O Interface References and Notes (1) Vacha, R.; Slavicek, P.; Mucha, M.; Finlayson-Pitts, B. J.; Jungwirth, P. J. Phys. Chem. A 2004, 108, 11573. (2) Roeselova, M.; Jungwirth, P.; Tobias, D. J.; Gerber, R. B. J. Phys. Chem. B 2003, 107, 12690. (3) Roeselova, M.; Vieceli, J.; Dang, L. X.; Garrett, B. C.; Tobias, D. J. J. Am. Chem. Soc. 2004, 126, 16308. (4) Wick, C. D.; Dang, L. X. J. Phys. Chem. B 2006, 110, 8917. (5) Knipping, E. M.; Lakin, M. J.; Foster, K. L.; Jungwirth, P.; Tobias, D. J.; Gerber, R. B.; Dabdub, D.; Finlayson-Pitts, B. J. Science 2000, 288, 301. (6) Chen, J.; Ehrenhauser, F. S.; Valsaraj, K. T.; Wornat, M. J. J. Phys. Chem. A 2006, 110, 9161. (7) Vacha, R.; Jungwirth, P.; Chenb, J.; Valsaraj, K. Phys. Chem. Chem. Phys. 2006, 8, 4461. (8) Vacha, R.; Cwiklik, L.; Rezac, J.; Hobza, P.; Jungwirth, P.; Valsaraj, K.; Bahr, S.; Kempter, V. J. Phys. Chem. A 2008, 112, 4942. (9) Valsaraj, K. T.; Thoma, G. J.; Reible, D. D.; Thibodeaux, L. J. Atmos. EnViron., Part A 1993, 27, 203. (10) Djikaev, Y. S.; Tabazadeh, A. J. Geophys. Res., [Atmos.] 2003, 108, 4689. (11) Saxena, P.; Hildemann, L. M. J. Atmos. Chem. 1996, 24, 57. (12) Lima, A. L. C.; Eglinton, T. I.; Reddy, C. M. EnViron. Sci. Technol. 2003, 37, 53. (13) Raja, S.; Valsaraj, K. T.; Ravikrishna, R. Chem. Abstr. 2003, 225, U821. (14) Raja, S.; Valsaraj, K. T. EnViron. Sci. Technol. 2004, 38, 763. (15) Vacha, R.; Jungwirth, P.; Chen, J.; Valsaraj, K. Phys. Chem. Chem. Phys. 2006, 8, 4461. (16) Vacha, R.; Cwiklik, L.; Rezac, J.; Honza, P.; Jungwirth, P.; Valsaraj, K.; Bahr, S.; Kempter, V. J. Phys. Chem. A 2008, 112, 4942. (17) Raja, S.; Yaccone, F. S.; Ravikrishna, R.; Valsaraj, K. T. J. Chem. Eng. Data 2002, 47, 1213. (18) Raja, S.; Valsaraj, K. T. J. Air Waste Manage. Assoc. 2005, 55, 1345. (19) Raja, S.; Raghunathan, R.; Yu, X. Y.; Lee, T. Y.; Chen, J.; Kommalapati, R. R.; Murugesan, K.; Shen, X.; Qingzhong, Y.; Valsaraj, K. T.; Collett, J. L. Atmos. EnViron. 2008, 42, 2048. (20) Gill, P. S.; Graedel, T. E.; Weschler, C. J. ReV. Geophys. 1983, 21, 903. (21) Cappiello, A.; De Simoni, E.; Fiorucci, C.; Mangani, F.; Palma, P.; Trufelli, H.; Decesari, S.; Facchini, M. C.; Fuzzi, S. EnViron. Sci. Technol. 2003, 37, 1229. (22) Rajamani, S.; Ghosh, T.; Garde, S. J. Chem. Phys. 2004, 120, 4457. (23) Bertram, A. K.; Ivanov, A. V.; Hunter, M.; Molina, L. T.; Molina, M. J. J. Phys. Chem. A 2001, 105, 9415. (24) Mmereki, B. T.; Donaldson, D. J.; Gilman, J. B.; Eliason, T. L.; Vaida, V. Atmos. EnViron. 2004, 38, 6091. (25) Donaldson, D. J.; Mmereki, B. T.; Chaudhuri, S. R.; Handley, S.; Oh, M. Faraday Discuss. 2005, 130, 227. (26) Chen, B.; Siepmann, J. I.; Klein, M. L. J. Am. Chem. Soc. 2002, 124, 12232. (27) Chen, B.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2002, 105, 3093.
J. Phys. Chem. C, Vol. 114, No. 34, 2010 14527 (28) Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2007, 111, 10790. (29) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2004, 108, 17596. (30) Wick, C. D.; Siepmann, J. I.; Theodorou, D. N. J. Am. Chem. Soc. 2005, 127, 12338. (31) Kamath, G.; Potoff, J. J. J. Phys. Chem. C 2007, 111, 1451. (32) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (33) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (34) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 2006, 110, 3555. (35) Siepmann, J. I. Mol. Phys. 1990, 70, 1145. (36) Siepmann, J. I.; Frenkel, D. Mol. Phys. 1992, 75, 59. (37) Wick, C. D.; Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 8008. (38) De Pablo, J. J.; Laso, M.; Suter, U. W. J. Chem. Phys. 1992, 96, 6157. (39) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: San Diego, CA, 2002. (40) Martin, M. G.; Siepmann, J. I. J. Am. Chem. Soc. 1997, 119, 8921. (41) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1987. (42) Lyubartsev, A. P.; Martsinovski, A. A.; Shevkunov, S. V.; Vorontsov-Velyaminov, P. N. J. Chem. Phys. 1992, 96, 1776. (43) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187. (44) Dang, L. X. J. Phys. Chem. B 2002, 106, 10388. (45) Wilson, M. A.; Pohorille, A. J. Phys. Chem. B 1997, 101, 3130. (46) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (47) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (48) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (49) Case, D. A.; et al. University of California: San Francisco, CA, 2008. (50) Hommelen, J. R. J. Colloid Sci. 1959, 14, 385. (51) Vochten, R.; Petre, G. J. Colloid Interface Sci. 1973, 42, 320. (52) Atkins, P.; de Paula, J. Atkins’ Physical Chemistry, 5th ed.; University Press: Oxford, U.K., 1994. (53) Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2002, 106, 6361. (54) Jungwirth, P.; Tobias, D. J. Chem. ReV. 2006, 106, 1259. (55) Tobias, D. J.; Hemminger, J. C. Science 2008, 319, 1197. (56) Chang, T. M.; Dang, L. X. Chem. ReV. 2006, 106, 1305. (57) Mmereki, B. T.; Chaudhuri, S. R.; Donaldson, D. J. J. Phys. Chem. A 2003, 107, 2264. (58) Wick, C. D.; Siepmann, J. I.; Schure, M. R. Anal. Chem. 2002, 74, 3518. (59) Sangster, J. Octanol-Water Partition Coefficients: Fundamentals and Physical Chemistry; John Wiley & Sons: West Sussex, U.K., 1997. (60) Chen, B.; Siepmann, J. I. J. Am. Chem. Soc. 2000, 122, 6464. (61) webbook.nist.gov. 2009.
JP1039578
