Lengthscale-Dependent Solvation and Density Fluctuations in n-Octane

Article pubs.acs.org/JPCB

Lengthscale-Dependent Solvation and Density Fluctuations in n‑Octane Eugene Wu and Shekhar Garde* Howard P. Isermann Department of Chemical and Biological Engineering and Center for Biotechnology and Interdisciplinary Studies, Rensselaer Polytechnic Institute, Troy, New York 12180, United States ABSTRACT: Much attention has been focused on the solvation and density fluctuations in water over the past decade. These studies have brought to light interesting physical features of solvation in condensed media, especially the dependence of solvation on the solute lengthscale, which may be general to many fluids. Here, we focus on the lengthscaledependent solvation and density fluctuations in n-octane, a simple organic liquid. Using extensive molecular simulations, we show a crossover in the solvation of solvophobic solutes with increasing size in n-octane, with the specifics of the crossover depending on the shape of the solute. Large lengthscale solvation, which is dominated by interface formation, emerges over subnanoscopic lengthscales. The crossover in n-octane occurs at smaller lengthscales than that in water. We connect the lengthscale of crossover to the range of attractive interactions in the fluid. The onset of the crossover is accompanied by the emergence of non-Gaussian tails in density fluctuations in solute shaped observation volumes. Simulations over a range of temperatures highlight a corresponding thermodynamic crossover in solvation. Qualitative similarities between lengthscale-dependent solvation in water, n-octane, and Lennard-Jones fluids highlight the generality of the underlying physics of solvation.

■

Lennard-Jones (LJ) fluid.10 Garde et al. also studied density fluctuations in polyethylene i-mer chains of various lengths from i = 8 to 64 and observed the onset of non-Gaussian density fluctuations in the tails of Pv(N) distributions and the corresponding dewetting of the solute surface in those polymeric systems.11 The same authors also studied 64-mer polymer chains and observed non-Gaussian behavior in the Pv(N) profiles.12 Ashbaugh and Pratt used scaled particle theory (SPT) supplemented with simulations to study the various aspects of the crossover in n-hexane. They point out similarities as well as interesting differences between solvation in water and in n-hexane. Specifically, the crossover occurs at a smaller solute size in n-hexane than in water, similar to that observed by Huang and Chandler and Huang et al. in a LJ fluid.10,13−15 Thus, the onset of dewetting, non-Gaussian behavior in density fluctuations, and changes in the sign of entropy all occur at smaller lengthscales in n-hexane and LJ fluid compared to that in water. Here we explore, in much detail, these solvation phenomena in liquid n-octane. We employ a combination of test particle insertions and the indirect umbrella sampling (INDUS) method to systematically study the solute-size-dependent crossover in solvation for spherical and cuboidal solutes. We

INTRODUCTION Lengthscale-dependent solvation of hydrophobic solutes in water displays several remarkable features.1−3 Small hydrophobic solutes are accommodated in the thermally fluctuating hydrogen-bonding network of water. The hydration of such solutes is unfavorable, and is dominated by a large negative entropy contribution. Conversely, the hydration of large solutes in water is governed by different physics, that of interface formation. This is also free energetically unfavorable but dominated by enthalpy rather than entropy.1,2 Further, there is a dewetting or drying transition with increasing size of the solute. This is especially relevant for an idealized hard-particle solute. As the hard-particle solute increases in size, water molecules in its vicinity lose significant interactions with neighboring waters, which leads to the dewetting of the solute surface.4,5 This lengthscale-dependent crossover is also reflected in the nature of density fluctuations in solute-shaped observation volumes in water. In small observation volumes, density fluctuations, as measured by the probability, Pv(N), of observing N water molecules in an observation volume of interest, v, are Gaussian.6 With increasing v, the Pv(N) becomes non-Gaussian, displaying fat tails in the low-N region.3 The lengthscale-dependent crossover described above has been studied primarily in water.1−4,7−9 However, the physics is general and applicable to any liquid with some solvent−solvent attractions. Indeed, Huang and Chandler employed Monte Carlo simulations and the Lum−Chandler−Weeks (LCW) theory2 to highlight the aspects of such a crossover, including the drying transition and non-Gaussian density fluctuations in a © XXXX American Chemical Society

Special Issue: Branka M. Ladanyi Festschrift Received: September 30, 2014 Revised: November 15, 2014

A

dx.doi.org/10.1021/jp509912v | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

examine the nature of density fluctuations, solvation thermodynamics, contact densities and dewetting, and asymptotic approach to macroscopic solvation in liquid n-octane and comment on its differences with hydrophobic hydration phenomena. Our study, combined with the past understanding of hydrophobic solvation sheds further light on the universality of lengthscale-dependent solvation in associating fluids.

0.5, 0.75, 1.0, and 2.0 nm. To adequately sample the range of N in each volume, we placed INDUS windows from Ñ * = −20 to Ñ * = 60, 70, 80, 110, 140, and 245 in increments of 5 for hv = 0.3, 0.325, 0.35, 0.375, 0.4, 0.5, 0.75, 1.0, and 2.0 nm, respectively. We used negative Ñ * windows to provide a strong potential to ensure complete emptying of the probe volume. A κ value of 0.5 kJ/mol was used in all simulations. The only exception is for hv = 2.0 nm, where the large probe volume can cause the n-octane slab to split apart. To maintain the slab configuration, κ was reduced to 0.25 kJ/mol at 298 K and to 0.125 kJ/mol at higher temperatures for all low Ñ * windows. For spherical probe volumes, we used radii, rv, of 0.3, 0.325, 0.35, 0.375, 0.4, 0.5, 0.75, and 1.0 nm to sample a range of solute lengthscales. We placed INDUS windows from Ñ * = −5 to 5, −6 to 6, −7 to 7, −8 to 8, and −10 to 10 in increments of 1 for rv = 0.3, 0.325, 0.35, 0.375, and 0.4, respectively. For rv = 0.5 nm, windows ranged from Ñ * = −16 to Ñ * = 16 in increments of 2. The window range was Ñ * = −20 to Ñ * = 60 and 130 in increments of 5 for rv = 0.75 nm and rv = 1.0 nm, respectively. A κ value of 1.0 kJ/mol was used for rv = 0.3, 0.325, 0.35, and 0.375 nm, and 0.75 kJ/mol was used for all other radii. Most of the above INDUS calculations were repeated at five other temperatures (310, 322, 334, 346, and 358 K) to quantify the thermodynamics of solvation, especially for large volumes. All INDUS simulations were run for 3 ns each with frames stored every 2 ps for the high Ñ * windows. For low Ñ * windows, frames were stored every 0.5 ps for analysis. The INDUS production run (500−3000 ps) was used for analysis, with the first 500 ps used for equilibration. Standard deviations for quantities were calculated using blocks of 500 ps.

■

METHODS System and Simulation Details. To calculate density fluctuations in large volumes, the original implementation of the INDUS method requires the presence of a buffering vapor zone, which allows the liquid phase to expand into that zone to keep the system pressure constant.3,16 Thus, our simulation system includes a 4.5 nm thick slab of n-octane (with a 4.5 nm × 4.5 nm xy cross section) in a 10.5 nm long box along the z direction. The n-octane molecules were represented explicitly using the AMBER94 force field.17 All simulations were performed in the canonical ensemble (N, V, T) using the molecular dynamics package GROMACS.18 The standard version was used for equilibrium simulations, while a version modified by Patel et al. was used for INDUS simulations.16 The leapfrog algorithm was used with a 2 fs timestep to integrate the equations of motion. The temperature was maintained using a Nosé−Hoover thermostat.19 In the AMBER94 force field, both carbons and hydrogens carry small partial charges.17 The particle mesh Ewald (PME) algorithm was used to calculate electrostatic interactions20 with a grid spacing of 0.12 nm and a cutoff of 1.0 nm. For LennardJones interactions, a 1.0 nm cutoff was applied and Lorentz− Berthelot mixing rules were used for cross Lennard-Jones interactions.21 The P-LINCS algorithm was used to constrain bonds containing hydrogen atoms.22 Six equilibrium simulations were run for a total of 5 ns each (500 ps equilibration and 4.5 ns production run) at different temperatures (298, 310, 322, 334, 346, and 358 K) with configurations being stored every 1 ps for analysis. The excess chemical potential of small hard-sphere solutes (with radius, rv, smaller than 0.5 nm) was calculated from equilibrium simulations using the Widom test particle insertion method.23 From the production run (500− 4500 ps), 1003 insertions were attempted every 2 ps for a total of 2 × 109 insertions at each temperature for each cavity size. Details of the INDUS Method. The technical details of the INDUS methods have been discussed by Patel et al. recently.3,16 The method calculates the probability distribution, Pv(N), in observation volume v, using a modified umbrella sampling algorithm. Applying an umbrella potential on particle numbers, N, leads to impulsive forces (because N is a discrete variable, i.e., integer), which is difficult to accommodate in a typical MD simulation algorithm.16 The INDUS method instead applies a harmonic potential, UI = 0.5κ(Ñ − Ñ *)2, on a continuous variable Ñ , that is obtained by coarse-graining N using a Gaussian truncation and recalculates unbiased distributions of N through appropriate transformations and WHAM.24 The other terms in the INDUS harmonic potential are κ, a spring-like constant, and Ñ *, the reference coarsegrained N for each window. We employed the INDUS method to calculate Pv(N) in spherical as well as cuboidal volumes placed at the center of the n-octane slab at 298 K. The INDUS potential is applied to only the carbon atoms of the n-octane molecule. The cuboidal observation volumes have a constant cross-sectional area of 2 nm × 2 nm and thicknesses, hv, of 0.3, 0.325, 0.35, 0.375, 0.4,

■

RESULTS AND DISCUSSION Density Fluctuations in Bulk n-Octane. We are interested in the nature of density fluctuations in volumes of different shapes and sizes in bulk n-octane at 298 K. To this end, Figure 1 shows the probability, Pv(N), of finding N carbon centers inside an observation volume, v, obtained using the INDUS method.16 We calculated Pv(N)’s in spherical volumes with radii ranging from molecular (0.3 nm) to nanoscopic (1.0 nm). We also calculated Pv(N)’s in cuboidal volumes with a cross section of 2 nm × 2 nm with thicknesses ranging from molecular (0.3 nm) to nanoscopic (2.0 nm). In small spheres as well as in thin cuboids, Pv(N) curves appear Gaussian (i.e., parabolic on a log scale). As the observation volumes increase in size, Pv(N) becomes nonGaussian, especially in the low-N tail of the distribution. To better compare the various curves on the same footing and to understand the Gaussian to non-Gaussian transition, we plot the data in Figure 1 by transforming the variable N to (N − ⟨N⟩)/⟨δN2⟩1/2 (where ⟨N⟩ is the mean and ⟨δN2⟩ is the variance of N) and ln Pv(N) to ln Pv(N) + ln (⟨δN2⟩1/2), which we denote as ln Pv for simplicity (see Figure 2). Figure 2 shows that, in all spherical and cuboidal observation volumes, the Pv(N) distribution is roughly Gaussian near the mean of the distribution. However, for spheres with radii 0.5 nm and larger and cuboids with thicknesses of 0.4 nm and larger, density fluctuations begin to display a clear nonGaussian nature, especially in the low-N tail. Such fat tails in the Pv(N) profile have been well studied in water and indicate the onset of a transition from small to large lengthscale solvation physics.2,3,16 Small solvophobic solutes are accommodated into cavities in the solvent that open up by random thermal B



Article

Figure 3. Probability distribution for finer grained spherical (A) and cuboidal (B) observation volumes. The Gaussian to non-Gaussian transition occurs between radii of 0.4 and 0.5 nm for spheres and between thicknesses of 0.325 and 0.35 nm for cuboids. The error bars are smaller than the size of the symbols.

Figure 1. Natural logarithm of the probability distribution ln Pv(N) of finding N carbon centers in a probe volume v using the INDUS method.16 Panel A shows results for spherical observation volumes with radii, rv, ranging from 0.3 to 1.0 nm. Panel B shows similar data for cuboidal observation volumes with a constant 2 nm × 2 nm cross section, and width hv varying from 0.3 to 1.0 nm. The error bars are smaller than the size of the symbols.

To better understand the lengthscale at which the Gaussian to non-Gaussian transition in the Pv(N) distribution begins, we performed the INDUS calculations over fine-grained observation volume sizes. Figure 3 shows that the Gaussian to nonGaussian transition occurs between radii of 0.4 and 0.5 nm for spherical volumes and between thicknesses of 0.325 and 0.35 nm for cuboidal volumes. Thus, the transition occurs at a smaller lengthscale for the cuboids than for the spheres. The cavity lengthscale at the onset of this transition also corresponds to the lengthscale at which interface formation becomes relevant. At what radius does the surface of a spherical cavity start displaying the properties of an extended surface? Or at what thickness does the surface of a cuboid with a 2 nm × 2 nm cross section begin to resemble an extended surface? For a cuboidal cavity, the two 2 nm × 2 nm surfaces are already sufficiently large. However, the carbon atoms of n-octane molecules can interact across the thickness of the cuboid, as is the case for a thin cuboid. As the thickness of the cuboid increases, such cross-interactions are reduced. The results in Figure 3 show that, once the thickness exceeds about one carbon−carbon σ (where σ is the Lennard-Jones parameter), approximately equal to 0.34 nm in our force field, those interactions are sufficiently small, and the density fluctuations show a distinctly non-Gaussian character. As we show below, local dewetting of n-octane from the cavity surface also occurs at approximately the same cavity thickness. In contrast to the cuboidal volume, the curvature of a spherical cavity leads to better interactions between the carbons of n-octane molecules in the vicinity of the cavity. As a result, a spherical cavity needs to grow to about 1.5 times the carbon− carbon σ to display the onset of small to large lengthscale transition. These observations are consistent with the results of Huang and Chandler for density fluctuations in a LJ fluid.10 Lengthscale-Dependent Solvation in n-Octane. We expect the dependence of the free energy of solvation of solutes with increasing size to be consistent with the nature of density fluctuations reported above. We calculated the excess chemical

Figure 2. Probability distribution, ln Pv, as a function of the scaled coordinate (N − ⟨N⟩)/⟨δN2⟩1/2 for (A) spherical and (B) cuboidal observation volumes. The corresponding Gaussian reference curve is also shown (black dashed line). The legend is in nm for both spheres (rv) and cuboids (hv). The error bars are smaller than the size of the symbols.

fluctuations.25 In contrast, the solvation of large solutes requires the formation of a vapor−liquid-like interface. The onset of interface formation being the dominant contribution is apparent in the fat low-N tails of the distribution.3 Such nonGaussian distributions are also observed in a LJ fluid by Huang and Chandler10 and are consistent with predictions from LCW theory.2 This non-Gaussian behavior is also seen for spherical solutes in n-alkanes and in a 64-mer polymer melt with increasing chain size.11,12 C



Article

potential, μex, of solvation of spherical and cuboidal solutes using μex = −kBT ln Pv(0), where Pv(0) is the probability of finding zero carbon centers in the observation volume of interest, and kBT is the product of the Boltzmann constant and temperature.6 We obtained Pv(0) using the INDUS method as described above for both spherical and cuboidal observation volumes. For small spheres (rv < 0.3 nm), we used direct test particle insertions23 because that method is computationally efficient. We use the parameter ξ to represent both the radii (rv) for spherical volumes and thicknesses (hv) for cuboidal volumes.

Figure 5. (A) The excess chemical potential per unit area, μex/Av, as a function of ξ (=rv for spheres; =hv for cuboids). The horizontal dashed line shows the surface tension of n-octane at 298 K obtained from experimental data in units of kJ/mol-nm2.26 The other dashed line is a guide to the eye. (B) μex as a function of total surface area, Av. Linear fits for Av > π nm2 for spheres (green) and >12 nm2 for cuboids (blue) were used to obtain the apparent surface tension values (units of kJ/ mol-nm2) in the large lengthscale limit. The error bars are smaller than the size of the symbols.

Figure 4. Excess chemical potential of cavity formation in kBT units, μex/kBT, as a function of solute size ξ (=rv for spheres; =hv for cuboids). The results from test particle insertions for spheres (red) and INDUS for spheres (green) and cuboids (blue) are shown.

Figure 4 shows the free energy of solvation of spherical and cuboidal solutes in kBT units as a function of the solute size ξ (either the radius for spheres or thickness for cuboids). The excess chemical potential of both spheres and cuboids increases with their size monotonically, although the nature of that increase appears qualitatively different for spheres and cuboids. For spheres, μex increases rapidly with size, whereas for cuboids (with a constant cross-sectional area) the increases becomes linear for larger solutes. Such a variation is not surprising given that the volumes and surface areas of these two objects vary differently with ξ. For spherical observation volumes, the trend in μex as a function of solute size is qualitatively consistent with the trend seen in water,2 in a LJ fluid,10,15 and in organic liquids.11,12 To better understand this variation and compare the results on the same footing, we plot the excess chemical potential normalized by solvent accessible surface area, μex/Av, for the two solutes, as shown in Figure 5A. When plotted in this manner, the dependence with solute size for both shapes is qualitatively similar, and is consistent with that observed for water and LJ fluids.1,2,10 For small solutes, μex/Av varies approximately linearly with ξ, and reaches a plateau in the large lengthscale limit. A crossover from small to large lengthscale solvation occurs at ∼0.5 nm for spherical solutes, which coincides with the Gaussian to non-Gaussian transition shown in Figure 3A. The crossover from small to large lengthscale solvation suggests that, for sufficiently large solutes, the thermodynamics of solvation is governed by interface formation. Figure 5B shows the variation of μex with the solvent accessible surface area of spherical and cuboidal solutes. Indeed, μex varies linearly with Av for larger solutes. For areas greater than π nm2 for spheres and greater than 12 nm2 for cuboids (which roughly corresponds to solutes with radii or thicknesses larger than 0.5 nm), we used a linear fit to the data, μex = γAv, to obtain the apparent surface tension of the solute−n-octane interface. For spheres and cuboids, the apparent surface tension is 12.1 and

10.3 kJ/mol-nm2, respectively. We expect the surface tension values to increase gradually for larger solutes (as expected from the curvature dependence of surface tension). Additionally, we expect that the apparent surface tension for cuboids is smaller than that for spheres of similar areas. Cuboidal solutes have eight corners, and as such, n-octane molecules near the corners can interact with each other through the corner region. Such interactions reduce the energetic penalty of creating the surface in the vicinity of the corner. The influence of the corner region will diminish for larger cuboids, or for cuboids where all dimensions are increased simultaneously (i.e., the cross section is not kept constant, as is done here). The values of surface tension obtained here by fitting μex linearly with area for spherical and cuboidal solutes are comparable to the surface tension of n-octane at 298 K (12.8 kJ/mol-nm2) obtained using fits to experimental surface tension data present in ref 26. A better comparison would be with the surface tension of the model of n-octane simulated here. To our knowledge, this value has not been reported for the explicit noctane force field used here. Using components of the pressure tensor from NVT simulations that contain an n-octane slab with two free interfaces,27,28 we estimate the surface tension to be about 7.7 kJ/mol-nm2, which is comparable to but smaller than that obtained using the INDUS data for spheres and cuboids. The apparent surface tension obtained from data from nanoscopic solutes depends on the extent of interface fluctuations in their vicinity as well as on the range of solvent−solvent attractions relative to solute size. These aspects were discussed in more detail previously in refs 3, 10, 15, and 29. The nature of the crossover in solvation is different between n-octane and water. The crossover with increasing solute size occurs more gradually in water1 compared to that in n-octane. We expect the lengthscale of the crossover to be related to the D



Article

range of solvent−solvent interactions. In n-octane, that range is governed by the van der Waals interactions, which is short compared to that for water, where electrostatics dominate. Rajamani et al. have shown that the crossover lengthscale is proportional to the Egelstaff−Widom lengthscale, the product of surface tension and compressibility, in fluids.7,30 This lengthscale is larger in fluids with long-range electrostatic interactions, compared to other fluids, as pointed out by Egelstaff and Widom.30 The Egelstaff−Widom lengthscale in noctane is roughly 70% of that in water (using compressibility and surface tension from ref 31 for n-octane), consistent with the smaller crossover lengthscale in octane. Thus, in essence, once a 0.5 nm cavity is formed in n-octane, the interactions that span across the cavity are minimal, and the free energy of cavity formation is governed by the apparent surface tension. This is consistent with the Gaussian to non-Gaussian transition observed in the 0.4−0.5 nm radius regime for spherical volumes in n-octane. Huang and Chandler observed a similar transition for spherical volumes in a LJ fluid.10 Dewetting of n-Octane and Its Relation to the Surface Tension. The scaled particle theory (SPT) perspective provides a connection between the local density of solvent near hard-particle solutes and the gradient of the free energy of solvation with respect to particle size (i.e., the force acting on the particle surface due to solvent). This force is integrated to obtain the total free energy of solvation.5,32 The variation of the local density of water near hydrophobic solutes has been studied extensively in the past, going as far back as 1973, in the seminal paper by Frank Stillinger.4 Near very small particles, the local water density is perturbed negligibly. With increasing solute size, water packs more efficiently, and the local density increases, reaching a maximum near methane-like solutes. Further increase in the solute size leads to breaking of the local hydrogen bonding network, and the hard-sphere solute dewets monotonically. SPT relates the size dependence of dewetting to the surface tension of water.4,5 Here we explore this connection for spherical and cuboidal solutes in n-octane. We used configurations generated by the INDUS method in conjunction with a perturbation treatment to obtain the density profile of carbons near hard spherical and cuboidal solutes in noctane. Thus, ⟨ρ(r)⟩HS =

⟨exp[β(UI − UHS)]ρ(r)⟩I ⟨exp[β(UI − UHS)]⟩I

Figure 6. Carbon density profiles surrounding (A) a hard-sphere and (B) a hard-cuboid in n-octane for solutes of increasing radii or thicknesses: 0.3 nm (red), 0.4 nm (green), 0.5 nm (blue), 0.75 nm (orange), and 1.0 nm (cyan). Density profiles for cuboids are calculated in the (z) direction perpendicular to the cross section and within volume elements having the same 2 nm × 2 nm (xy) cross section. The large symbols highlight the contact density (ρc) for each hard-solute size. These points are connected with a black dashed line as a guide to the eye.

Figure 6 shows the density profiles of n-octane carbons as a function of distance from the center of the observation volume. For spheres, these are radial density profiles of n-octane carbon atoms surrounding the solute, whereas, for cuboids, they are one-particle density profiles in the direction normal to the 2 nm × 2 nm surface. The density profiles are plotted around five spherical solutes and five cuboids sampling radii or thicknesses from 0.3 to 1 nm. The contact density (of interest to SPT) is highlighted with large points. Near spheres of radii 0.3 nm or near a cuboid with a thickness of 0.3 nm, n-octane molecules pack well and display layering. As the solute size increases, the density drops monotonically, and the layering becomes less pronounced. To understand the solute size dependence of density, especially the contact density, for smaller spherical solutes, we used test particle insertions in pure n-octane liquid. We did not perform similar insertions or INDUS simulations of thinner cuboidal solutes with the 2 nm × 2 nm cross section. Contact density profiles for spheres (red and green) and cuboids (blue) are shown in Figure 7. The variation of contact density, ρc, for small spheres in organic solvents contains interesting features that have been understood well previously.11,13,14,33 For example, the kink at 0.07 nm corresponds to half the C−C bond length. Two carbons diametrically across can accommodate a cavity of r = 0.07 nm between them. As the cavity size increases, there is a local minimum in density. Although these variations are interesting, they are not relevant to the solvation of realistic solutes, as a point hard-particle solute (rHS = 0) requires a cavity of about 0.2 nm radius (half of the hard-sphere size of carbon) to accommodate it. A more detailed treatment for this regime can be found in refs 13, 14 and 33. Beyond this radius (about 0.2 nm), the variation in the carbon contact density near hard-sphere solutes is qualitatively similar to that described above for water, displaying a maximum followed by a monotonic dewetting of larger solutes.4 The

(1)

where subscripts I and HS denote the INDUS (reference) and hard-particle (perturbed) ensembles, respectively, β is 1/kBT, ρ(r) is the carbon density profile in the vicinity of either spherical or cuboidal observation solutes, and U is the potential energy for a given configuration in an INDUS simulation. The UHS in the above equation is infinity, and the weight factor is zero for any configuration that contains a nonzero number of solvent atoms in the observation volume of interest. This is a consequence of perturbing from a soft (INDUS) to a hardparticle system. To obtain good statistics, therefore, only those configurations in which the observation volume of interest is completely empty (N = 0) are useful. To sample such configurations, we used INDUS simulations with low Ñ * values. Specifically, we applied this perturbation scheme to INDUS windows with the two lowest Ñ * values for spherical volumes and four lowest Ñ * windows for cuboidal volumes along with the appropriate parameters for UI in eq 1 to calculate the density. E



Article

Figure 8. Free energy of solvation per unit area, μex/Av, for (A) spheres and (B) cuboids at temperatures of 298, 310, 322, 334, 346, and 358 K (blue to red). In panel A, results from test particle insertions and the INDUS method are plotted together for each isotherm.

Figure 7. Contact density, ρc, as a function of r for hard spheres and h for hard cuboids. Results from test particle insertions for spheres (red) and the INDUS method (spheres in green and cuboids in blue) are shown. The green line is a 1/r fit for larger spheres. The dashed blue line serves as a guide to the eye for cuboids. The solid blue horizontal line shows the point at which ρc does not vary with h for cuboids. The bulk simulation density (a dashed black line ρbulk) is shown for reference. Values for the apparent surface tension are in units of kJ/ mol-nm2.

cuboidal observation volumes (Figure 8B) at large thicknesses. For small spherical observation volumes, μex/Av is small, and its temperature dependence is not clearly noticeable on this scale. We discuss it in detail later in this section. As above, we calculated the apparent surface tension, γ, for spheres and cuboids at each temperature by fitting μex versus Av data in the large lengthscale limit (rv ≥ 0.5 nm and hv ≥ 0.5 nm). Figure 9 shows a comparison of the calculated surface tensions with the reported experimental surface tension of noctane26 as a function of temperature.

maximum in the contact density, however, occurs for smaller solutes in n-octane, consistent with the crossover lengthscale being smaller than that in water.13,14 Because the smallest cuboidal solute studied here is 0.3 nm thick, we only observe its dewetting with increasing thickness. For thinner cuboids, we expect that features qualitatively similar to spherical solutes would be observed. For sufficiently large solutes, the solvation of which can be described using a surface tension model, SPT relates the variation of contact density with solute size to the surface tension of the solvent.25,32 For spherical solutes with a hardsphere radius, r, 2γ 1 ρc,sphere (r ) = kBT r (2)

Figure 9. Apparent surface tension obtained from simulations as a function of temperature for spherical and cuboidal observation volumes. Experimental surface tensions (Expt) are shown as well.26 The temperature dependence of surface tension is fit to γ = γ0(1 − T/ Tc)α, where γ0 is the surface tension coefficient, Tc is the critical temperature, and α = 1.26.34

and for l × l × h cuboidal solutes (where l is the side of the cross section and h is the thickness of a hard-cuboid), 4γ 1 ρc,cuboid (h) = kBT l (3) Thus, for large spherical solutes, the contact density varies inversely with the hard-sphere radius. In contrast, for the cuboids considered here (when they are sufficiently thick), the contact density is expected to be independent of the hardcuboid thickness. For spherical observation volumes, we use a 1/r fit to INDUS data and eq 2 to obtain a γ of 12.8 kJ/molnm2. For cuboidal observation volumes with l = 2.0 nm, we find the point at which ρc is independent of ξ = h and, with eq 3, obtained γ to be 11.4 kJ/mol-nm2. The γ values for both spheres and cuboids are comparable to those obtained from fitting μex versus Av, as seen in Figure 5B. Thermodynamics of the Crossover in n-Octane. How does the solute lengthscale-dependent crossover behavior depend on temperature? As above, we used a combination of test particle insertions (for small spheres) and INDUS (for larger spheres and cuboids) to calculate the free energy of solvation at six different temperatures ranging from 298 to 358 K. Figure 8A shows that at large solute lengthscale, μex/Av plateaus at each temperature with the plateau value decreasing with increasing temperature. The same trend is observed for

Consistent with results reported above, the cuboidal observation volumes have a lower γ at each temperature compared to spherical observation volumes. The surface tension for both shapes decreases monotonically with temperature. Following ref 34, we estimate the critical temperature, Tc, of n-octane by fitting the surface tension as a function of

Figure 10. Enthalpic and entropic contributions to the apparent surface tension in the large lengthscale limit, obtained from appropriate temperature derivatives of γ = γ0(1 − T/Tc)α. F



■

temperature using γ = γ0(1 − T/Tc)α, where γ0 is the surface tension coefficient and α = 1.26. We estimate Tc to be 590 and 580 K using data for spherical and cuboidal observation volumes, respectively, which is within 5% of the experimental value of 569 K for n-octane.26 The temperature dependence of γ allows us to calculate the enthalpic and entropic contributions to it. We used the same fit as above (γ = γ0(1 − T/Tc)α) to the temperature dependent γ in Figure 9 for spherical volumes to obtain these contributions, as shown in Figure 10, such that γ=

ΔG ΔH ΔS = −T Av Av Av

Article

CONCLUSIONS We employed molecular dynamics simulations coupled with a specialized umbrella sampling technique (INDUS) and test particle insertions to study the solute lengthscale-dependent crossover in solvation in n-octane. Specifically, we calculated solvent density fluctuations and free energies of cavity formation, and we studied aspects of solute surface dewetting as a function of solute size. Consistent with other studies, our calculations reveal the qualitative similarities between the nature of solvation in organic liquids,11−14 LJ fluid models,10,15 and water.2,3,35 Solvation of small cavities in n-octane is entropically unfavorable, although the magnitude of entropy is much smaller than that in water. Correspondingly, density fluctuations in small volumes are Gaussian in nature. With increasing solute size, density fluctuations become nonGaussian, especially in the low-N tail, and the thermodynamics of solvation becomes increasingly governed by the physics of interface formation. The local density of carbons of n-octane near the solute surface decreases with increasing solute size. The surface tension estimated using the observed dewetting (within the SPT perspective) is consistent with that obtained independently from the area dependence of free energy. The key features of lengthscale-dependent solvationthe Gaussian to non-Gaussian transition, dewetting of the solute surface, and the change in the nature of thermodynamicspresent in both water2−4 and n-octane point to the general underlying physics of solvation. The difference between solvation in n-octane and water is the location of the crossover, which occurs at a smaller lengthscale in n-octane than in water. We argue that this difference is related to the range of interactions in the two liquids, van der Waals in n-octane (short) and electrostatic in water (long). Quantifying the density fluctuations and the associated quantities in n-octane and other organic liquids near interfaces (e.g., graphene sheets, nanotubes, nanoparticles) may open up avenues for understanding the solvent−solute coupling in nanocomposites36,37 and their relation to macroscopic properties. Density fluctuations at soft liquid−liquid interfaces are also of interest from the perspective of self-assembly in interfacial environments.38−41

(4)

As expected from the physics of interface formation, the free energy of surface formation (γ) is dominated by unfavorable enthalpy and favored by entropy, with the magnitude of enthalpy being approximately twice that of entropy. The variation of chemical potential of small spheres with temperature is not noticeable on the scale shown in Figure 8A. To focus on thermodynamics at smaller lengthscales, we follow the same treatment used by Huang and Chandler35 and plot the free energy of solvation (obtained either using test particle insertions or the INDUS method) normalized by cavity volume, μex/v, over the range of temperatures studied, as shown in Figure 11.

■

Figure 11. Free energy of solvation per unit volume, μex/v, calculated from test particle insertions (red) or the INDUS method (green) for spherical volumes. The radius of each cavity is labeled to the right of each curve. Lines serves as guides to the eye.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

Unlike that for large cavities, μex increases with temperature for small cavities, suggesting an unfavorable entropy of solvation in n-octane, qualitatively similar to that for small hydrophobic solutes in water. Specifically, for cavities with a radius of 0.1 nm, μex increases monotonically with temperature over the range studied here. For cavities with a radius of 0.2 nm, μex displays a parabolic dependence on temperature, whereas, for cavities with radii of 0.3 nm and larger, we observe large lengthscale-like solvation thermodynamics. If the lengthscale for thermodynamic crossover is defined as the size at which solvation entropy changes sign from negative to positive, that lengthscale is between 0.2 and 0.3 nm in n-octane. Our observation is qualitatively consistent with that of Ashbaugh and Pratt, where they observed a negative surface entropy for spheres with radii smaller than 0.25 nm in n-hexane using scale particle theory.13,14 Interestingly, this size is even smaller than the cavity size at which the onset of non-Gaussian behavior is observed, as shown above.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Prof. Rajesh Khare and Dr. Ketan S. Khare for helpful discussions during the beginning stage of this project and comments on an earlier version of this manuscript. We also thank Drs. Hari Acharya, Srivathsan Vembanur, and Vasudevan Venkateshwaran for helpful discussions. Finally, we thank Prof. Amish J. Patel for his implementation of INDUS in GROMACS and for his comments on an earlier version of this manuscript. This work was partially supported by the National Science Foundation (Grant Nos. CBET-0967937 and CBET-1159990).

■

REFERENCES

(1) Chandler, D. Interfaces and the Driving Force of Hydrophobic Assembly. Nature 2005, 437, 640−647.

G



Article

(2) Lum, K.; Chandler, D.; Weeks, J. D. Hydrophobicity at Small and Large Length Scales. J. Phys. Chem. B 1999, 103, 4570−4577. (3) Patel, A. J.; Varilly, P.; Chandler, D. Fluctuations of Water Near Extended Hydrophobic and Hydrophilic Surfaces. J. Phys. Chem. B 2010, 114, 1632−1637. (4) Stillinger, F. H. Structure in Aqueous Solutions of Nonpolar Solutes from the Standpoint of Scaled-Particle Theory. J. Solution Chem. 1973, 2, 141−158. (5) Ashbaugh, H. S.; Pratt, L. R. Colloquium: Scaled Particle Theory and the Length Scales of Hydrophobicity. Rev. Mod. Phys. 2006, 78, 159−178. (6) Hummer, G.; Garde, S.; García, A. E.; Pohorille, A.; Pratt, L. R. An Information Theory Model of Hydrophobic Interactions. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 8951−8955. (7) Rajamani, S.; Truskett, T. M.; Garde, S. Hydrophobic Hydration from Small to Large Lengthscales: Understanding and Manipulating the Crossover. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 9475−9480. (8) Varilly, P.; Patel, A. J.; Chandler, D. An Improved CoarseGrained Model of Solvation and the Hydrophobic Effect. J. Chem. Phys. 2011, 134, 074109. (9) Patel, A. J.; Varilly, P.; Jamadagni, S. N.; Acharya, H.; Garde, S.; Chandler, D. Extended Surfaces Module Hydrophobic Interactions of Neighboring Solutes. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 17678− 17683. (10) Huang, D. M.; Chandler, D. Cavity Formation and the Drying Transition in the Lennard-Jones Fluid. Phys. Rev. E 2000, 61, 1501− 1506. (11) Garde, S.; Hummer, G.; Khare, R.; Pratt, L. R. Effect of Chain Length on Microscopic Density Fluctuations and Solvation in Polymeric Fluids. Polym. Mater.: Sci. Eng. (1983-2001) 2001, 85, 449−451. (12) Garde, S.; Khare, R.; Hummer, G. Microscopic Density Fluctuations and Solvation in Polymeric Fluids. J. Chem. Phys. 2000, 112, 1574−1578. (13) Ashbaugh, H. S.; Pratt, L. R. Contrasting Nonaqueous Against Aqueous Solvation on the Basis of Scaled-Particle Theory. J. Phys. Chem. B 2007, 111, 9330−9336. (14) Ashbaugh, H. S. Entropy Crossover from Molecular to Macroscopic Cavity Hydration. Chem. Phys. Lett. 2009, 447, 109−111. (15) Huang, D. M.; Geissler, P. L.; Chandler, D. Scaling of Hydrophobic Solvation Free Energies. J. Phys. Chem. B 2001, 105, 6704−6709. (16) Patel, A. J.; Varilly, P.; Chandler, D.; Garde, S. Quantifying Density Fluctuations in Volumes of All Shapes and Sizes Using Indirect Umbrella Sampling. J. Stat. Phys. 2011, 145, 265−275. (17) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Kenneth, M.; Merz, J.; Gerguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (18) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; et al. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845−854. (19) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (20) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (21) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (22) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 116−122. (23) Widom, B. Some Topics in the Theory of Fluids. J. Chem. Phys. 1963, 39, 2808−2812. (24) Roux, B. The Calculation of the Potential of Mean Force Using Computer Simulations. Comput. Phys. Commun. 1995, 91, 275−282.

(25) Pratt, L. R.; Pohorille, A. Theory of Hydrophobicity: Transient Cavities in Molecular Liquids. Proc. Natl. Acad. Sci. U.S.A. 1992, 89, 2995−2999. (26) Grigoryev, B. A.; Nemzer, B. V.; Kurumov, D. S.; Sengers, J. V. Surface Tension of Normal Pentane, Hexane, Heptane, and Octane. Int. J. Thermophys. 1992, 13, 453−464. (27) Alejandre, J.; Tildesley, D. J.; Chapela, G. A. Molecular Dynamics Simulation of the Orthobaric Densities and Surface Tension of Water. J. Chem. Phys. 1995, 102, 4574−4583. (28) Godawat, R.; Jamadagni, S. N.; Garde, S. Unfolding of Hydrophobic Polymers in Guanidinium Chloride Solutions. J. Phys. Chem. B 2010, 114, 2246−2254. (29) Mittal, J.; Hummer, G. Static and Dynamic Correlations in Water at Hydrophobic Interfaces. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 20130−20135. (30) Egelstaff, P. A.; Widom, B. Liquid Surface Tension near the Triple Point. J. Chem. Phys. 1970, 53, 2667−2669. (31) Sanchez, I. C. Liquids: Surface Tension, Compressibility, and Invariants. J. Chem. Phys. 1983, 79, 405−415. (32) Reiss, H.; Frisch, H. L.; Lebowitz, J. L. Statistical Mechanics of Rigid Spheres. J. Chem. Phys. 1959, 31, 369−380. (33) Jain, A.; Ashbaugh, H. S. Digging a Hole: Scaled-particle Theory and Cavity Solvation in Organic Solvents. J. Chem. Phys. 2008, 129, 174505-1−174505-8. (34) Adam, N. K. The Physics and Chemistry of Surfaces; 3rd ed., Oxford University Press: London, U.K. 1941. (35) Huang, D. M.; Chandler, D. Temperature and Length Scale Dependence of Hydrophobic Effects and Their Possible Implications for Protein Folding. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 8324−8327. (36) Khare, K. S.; Khare, R. Effect of Carbon Nanotube Disperion on Glass Transition in Cross-Linked Epoxy-Carbon Nanotube Nanocomposites: Role of Interfacial Interactions. J. Phys. Chem. B 2013, 117, 7444−7454. (37) Khare, K. S.; Khabaz, F.; Khare, R. Effect of Carbon Nanotube Functionalization on Mechanical and Thermal Properties of CrossLinked Epoxy-Carbon Nanotube Nanocomposites: Role of Strengthening the Interfacial Interactions. ACS Appl. Mater. Interfaces 2014, 6, 6098−6110. (38) Patel, H. A.; Nauman, E. B.; Garde, S. Molecular Structure and Hydrophobic Solvation Thermodynamics at an Octane-Water Interface. J. Chem. Phys. 2003, 119, 9199−9206. (39) Chowdhary, J.; Ladanyi, B. M. Water-Hydrocarbon Interfaces: Effect of Hydrocarbon Branching on Interfacial Structure. J. Phys. Chem. B 2006, 110, 15442−15453. (40) Chowdhary, J.; Ladanyi, B. M. Surface Fluctuations at the Liquid-liquid Interface. Phys. Rev. E 2008, 77, 031609-1−031609-14. (41) Chowdhary, J.; Ladanyi, B. M. Hydrogen Bond Dynamics at the Water/Hydrocarbon Interface. J. Phys. Chem. B 2009, 113, 4045− 4053.

H


Lengthscale-Dependent Solvation and Density Fluctuations in n-Octane

Recommend Documents