Solubility and Molecular Conformations of n-Alkane Chains in Water

J. Phys. Chem. B 2009, 113, 6405–6414

6405

Solubility and Molecular Conformations of n-Alkane Chains in Water Andrew L. Ferguson, Pablo G. Debenedetti, and Athanassios Z. Panagiotopoulos* Department of Chemical Engineering, Princeton UniVersity, Princeton, New Jersey 08544 ReceiVed: December 19, 2008

We employ molecular dynamics simulations to study the solubility and molecular conformations of n-alkane chains in water. We find nearly exponential decrease in solubility with carbon number up to n-eicosane (C20), and excellent agreement with experiment up to n-dodecane (C12). We detect no sharp break in the dependence of the solubility upon carbon number. A free energy landscape analysis of the chain conformations reveals remarkable similarities between the ideal gas and solvated phase landscapes, suggesting that solvated chain conformations are driven primarily by ideal gas statistics. We find no evidence for hydrophobic collapse of n-alkane chains shorter than n-eicosane (C20). The primary effect of the solvent is the appearance of a barrier on the order kBT, not present in the ideal gas, between the free energy basins corresponding to compact and extended chain conformations, and destabilization of the most extended conformations. Our findings are robust to nontrivial modification of the potential model, suggesting that the absence of strong solvent effects on the free energy landscapes is fundamental to relatively short (e20-mer) chains composed of small hydrophobic monomers, and does not depend on the precise nature of the chain interactions. I. Introduction The preferential adoption of collapsed conformations by sufficiently long hydrophobic molecules in water has long been interpreted1,2 as a classic manifestation of the hydrophobic effect3 and has been observed by simulation.2,4,5 Hydrophobic collapse is an important driving force in protein folding6,7 with recent experimental8,9 and simulation9 studies providing strong support for this view. Normal alkanes, which are strongly hydrophobic, have been the subject of numerous simulation studies investigating their structure and solubility in aqueous solution.2,4,5,10-18 The solubility of n-alkanes is of practical import for chemical process design,19 and its relationship to the underlying chain conformations is also of fundamental scientific interest. Numerous experimental studies have been conducted to determine the solubility of n-alkanes in water.20-29 As illustrated in Figure 1, the agreement between different sources is poor for chains heavier than n-undecane (C11). The approximately exponential decrease in the equilibrium aqueous concentration of n-alkanes with carbon number renders its experimental determination exceedingly difficult for long chains,23 explaining the discrepancies between the reported solubilities from different sources. The apparent break in some solubility data sets at n-undecane (C11) has been speculatively attributed to the presence of metastable colloidal suspensions of n-alkanes leading to an artificial enhancement in the measured solubility,23,24 or to a change in the equilibrium microstructure from an extended to a collapsed conformation.24 Group contribution models developed specifically to describe hydrocarbons in water have been successful in correlating n-alkane solubilities in the short-chain region where reliable experimental data exist,30,31 and the results of one such model31 are presented in Figure 1. If, however, a microscopic conformational change does occur, then the extrapolation of such models is unjustified.19 In more recent years, a limited number of molecular simulation studies have been conducted to obtain solubilities of * Corresponding author. E-mail: [email protected].

Figure 1. Solubilities of n-alkanes in water at 298 K and 1 bar. Our results are plotted as red circles; error bars corresponding to the formal uncertainties in the Henry’s constant results described in the caption to Figure 2 are smaller than symbol size. Experimental data: green squares, Sander (2000);20 blue squares, Mackay and Shiu (1981);28 purple up triangle, Tsonopoulos (1999); 24 gray up triangle, Franks (1966);26 pink down triangle, Baker (1959);25 blue down triangle, Khadikar et al. (2003); brown left-pointing triangle, Tolls et al. (2002);23 green left-pointing triangle, McAuliffe (1966, 1969);22,29 yellow rightpointing triangle, Sutton and Calder (1974).27 Group contribution model: solid line, Plyasunov and Shock (2000).31

n-alkanes in water.14,16-18 Studies of methane (C1), ethane (C2), and n-butane (C4) show good agreement with experiment,16,18 but simulations of heavier chains have been less successful.14,16 The poor performance for longer chains has been attributed to inadequacies in the potential functions implemented,14 but it is also due to the difficulties in performing free energy calculations for longer chains.14,16 Simulations, on the other hand, can provide information on the thermodynamic properties as well as the microscopic conformations of the participating molecules, and are not subject to artifacts such as contamination, permitting,

10.1021/jp811229q CCC: $40.75  2009 American Chemical Society Published on Web 04/10/2009

6406

J. Phys. Chem. B, Vol. 113, No. 18, 2009

Ferguson et al. In section V, we analyze free energy landscapes of ideal gas and solvated chains, parametrized by the radius of gyration, to determine the influence of solvent interactions on chain conformations. In section VI we study solute-induced cavitation. In section VII, we summarize our results and present our conclusions. Additional results are provided in the Supporting Information. II. Simulation Methodology

Figure 2. Calculated Henry’s constants for n-alkanes in water at 298 K and 1 bar. The uncertainty associated with each point was estimated by taking the standard deviation of the results of five independent Widom simulations for n-tetradecane (C14), and is smaller than the symbol size. The n-docosane (C22) error bar represents the largest formal uncertainty associated with any single result, computed from the propagation the uncertainties through the incremental Widom insertion calculation.

in principle, the determination of n-alkane solubilities under precisely specified conditions. In this work, we compute the aqueous solubilities of n-alkanes ranging from methane (C1) to n-docosane (C22) at 298 K and 1 bar. We employ replica exchange molecular dynamics (REMD) simulations32 and the incremental Widom insertion technique33 to obviate sampling difficulties. To the best of our knowledge, we have made the first reliable determination from either experiment or simulation of Henry’s constant andsunder appropriate assumptions and the specific conditions described belowsthe solubility of n-alkanes heavier than n-hexadecane (C16) in water. To leading order, our results show an exponential decrease of solubility with carbon number, which results in a linear relationship when plotted on semilogarithmic axes. We observe no evidence of a discontinuity in the slope of this trend, but do observe interesting second-order deviations from linearity attributed to a maximum in Henry’s constant around noctadecane (C18). We compare the ensemble of structures of ideal gas and solvated n-alkane chains up to n-eicosane (C20), where we use the term “ideal gas” to describe simulations of a single n-alkane chain that interacts only with itself. This comparison enables us to identify the effect of solute-solvent interactions on chain structure, relate trends in solubility to chain conformations, and investigate the solvent contributions to chain collapse. We find the observed collapse of sufficiently long chains to be largely an ideal gas effect, with hydrophobic collapse apparently limited to chains longer than the upper limit considered here. The primary effect of the solvent is the appearance of a free energy barrier separating the collapsed and extended chain conformations, and a destabilization of the most extended conformations. These findings are robust to nontrivial modifications in the force field, suggesting that collapse driven by ideal gas statistics is fundamental to polymer chains composed of small hydrophobic monomers. The structure of this paper is as follows. In section II we provide details of our simulation methodology. In section III, we discuss the results of the solubility calculations. In section IV, we compare the chain structure in the ideal gas and in water.

Our use of REMD was precipitated by our experience of structural hysteresis in canonical molecular dynamics (MD) simulations of n-hexadecane (C16) and longer chains: simulations initialized with the chain in an extended (compact) conformation remained extended (compact) for the duration of a 1 ns trajectory. A 1 ns REMD simulation of n-hexadecane (C16), in contrast, sampled compact conformations never visited in an MD simulation of the same duration and exhibited far more frequent transitions between extended and compact conformations (Figure S1 in the Supporting Information). REMD facilitates escape from such kinetic traps and extraction of proper thermodynamic averages. Polymeric chains longer than a few monomers in length present an additional sampling difficulty with regard to the solubility calculation via the conventional Widom method in a dense solvent, since the vast majority of test insertions of the chain result in unphysical atomic overlaps which yield large interaction energies contributing very little to the canonical average. More sophisticated sampling techniques are necessary for chains longer than n-butane (C4).14,16 The incremental Widom insertion method described below permits, in principle, free energy calculations for arbitrarily long chains. II.A. Replica Exchange Molecular Dynamics. Simulations in the aqueous phase were conducted using an in-house REMD code and implementing the TraPPE potential for the n-alkane chains34 and the SPC/E model of water.35 A modified potential for the chains was also considered in which the bending and torsional contributions of the TraPPE model were turned off, and will henceforward be referred to as the “soft” potential. Throughout this paper, chains governed by the “soft” potential will henceforth be referred to by their carbon number only (e.g., C10), whereas those simulated with the full TraPPE potential will be referred to by name followed by the carbon number in parentheses (e.g., n-decane (C10)). Integration of the equations of motion was performed with the velocity Verlet36 algorithm in Cartesian coordinates with a time step of 1 fs. A Nose´-Hoover thermostat chain37,38 and Andersen-Hoover barostat39 were implemented,40 and the RATTLE algorithm41 was used to maintain bond and angle constraints. The Ewald summation technique42,43 was used to treat the electrostatic interactions with a real space cutoff identical to the Lennard-Jones cutoff of 9.875 Å; the reciprocal space cutoff was chosen to give similar errors in the reciprocal and real space force calculations.44,45 Lorentz-Berthelot combining rules were employed to compute interaction parameters between the chains and water molecules.46 Starting configurations were prepared by inserting the chains into a cubic box with periodic boundary conditions filled with water molecules at 298 K and a density of 1 g/cm3 where the size of the simulation cell was sufficiently large to preclude direct Lennard-Jones or real space Ewald interactions between periodic images of the chain. Water molecules closest to the chain were deleted one by one until a number of solvent molecules equal to the number of carbon atoms in the chain had been eliminated.10 High energy overlaps were removed

n-Alkane Chains in Water

J. Phys. Chem. B, Vol. 113, No. 18, 2009 6407

using steepest descent energy minimization. A 100 ps NPT run equilibrated the system to 298 K and 1 bar followed by a further NVT equilibration run at 298 K of at least 200 ps. Production runs were conducted under NVT conditions where the REMD temperature replicas were initialized from the final state of the equilibration run and given 10 ps to adapt to the new temperature before attempting temperature swaps between neighboring replicas every 0.1 ps. The replica exchange temperature ladder was iteratively refined over the range 298-350 K to result in approximately equal temperature acceptances of ∼20% between neighboring replicas;32,47 data were harvested from replicas occupying the 298 K rung. Simulations used to construct free energy landscapes (section II.D) were conducted for 1 ns, and extended by 0.5 ns in the case of the chains simulated with the full TraPPE potential to improve the quality of the landscapes; sampling was performed every 10 fs. In the case of the solubility calculations, 100 Widom test insertions were conducted every 10 fs over the course of a 1 ns trajectory. II.B. Ideal Gas Monte Carlo. To obviate the temperature control problems associated with molecular dynamics simulations of isolated molecules, simulations of n-alkane chains in the ideal gas phase at 298 K were conducted using an in-house Monte Carlo (MC) code. The only trial move performed was chain regrowth using continuous space conformational biasing to improve sampling.45 The sampling frequency was set to the reciprocal of the regrowth acceptance rate which, for the longer chains considered, required runs of up to 1 million MC steps to gather sufficient conformational statistics. Five independent 100 000 step simulations were run for each solubility calculation, with 100 incremental Widom insertions performed at a frequency equal to the reciprocal of the chain regrowth acceptance rate, resulting in around 5 million insertions per simulation. II.C. Incremental Widom Insertion Method. Aqueous solubility of the n-alkane chains was calculated using a generalization of the incremental Widom insertion technique of Kumar et al.,33,48 where we treat the CH4/CH3/CH2 united atoms as monomers. The heteropolymeric nature of n-alkane chains necessitates a slight generalization of the method. The calculation of the excess chemical potential of methane (C1) requires consideration of a CH4 monomer, whereas the calculation for longer chains requires a CH3 monomer in the starting position. As such, Widom insertions for i ) 1 must be conducted with both CH4 and CH3 monomers, and by extension with both CH3 and CH2 monomers for i ) 2,..., N. Monomeric Widom test insertions were conducted in the REMD simulations of an “uncapped” n-alkane chainsa chain in which the terminal monomer has been removedscomprising (i - 1) monomers at infinite dilution in aqueous solvent, to compute the quantity33,48

exp(-βµ˜ i) ) 〈exp(-β∆Ui)〉FULL

(1)

where µ˜ i is the nonideal contribution to the chemical potential of monomer i in the fully interacting system, ∆Ui is the interaction energy of the Widom test particle with the rest of the system, and β ) 1/(kBT), where kB is Boltzmann’s constant and T is the absolute temperature. The angled brackets on the right-hand side represent the double average over insertion locations and canonical microstates of the fully interacting system, for which sufficient statistics may be gathered by repeated Widom test particle insertions until convergence is attained. Similar insertions were conducted in the ideal gas phase MC simulations to compute the intrachain contribution to the chemical potential of monomer i in the ideal gas state at the same temperature. In these calculations ∆Ui is the interaction

energy of the Widom test particle with the rest of the ideal gas system which, by definition, reduces to the intramolecular interaction of the test particle with the rest of the isolated chain. Conformational biasing45 of the insertion locations toward favorable bending and torsional angles improved sampling. From values of µ˜ i for each monomer in the n-alkane chain, the excess ex , is computed chemical potential of the complete chain, µchain as33 N

ex µchain ) -kBT

〈exp(-β∆U )〉

FULL ∑ ln 〈exp(-β∆Uii)〉IDEAL

(2)

i)1

where we define the excess chemical potential as the difference between the chemical potential of a solute in water and that in an ideal gas mixture at the same composition, density, and ex is related to temperature. In the limit of infinite dilution, µchain 49 the Henry’s constant for the chain, Hchain, by the thermodynamic equation ex Hchain ) lim kBTFsolvent exp(βµchain ) xchainf0

(3)

where xchain is the mole fraction of chain in solution and Fsolvent is the number density of the pure solvent. Hchain is the proportionality constant between the fugacity of the solvated aq , and its mole fraction.50 Under conditions of chain, ˆf chain vapor-liquid equilibrium, the fugacity of the n-alkane in vap , aqueous solution is identical to that in the vapor phase, ˆf chain and at atmospheric pressure the latter maysunder the excellent approximation that the hydrocarbon vapor behaves ideally under these conditionssbe replaced by the vapor pressure, Pvap chain. These manipulations lead to the following expression, linking the equilibrium solubility to the n-alkane vapor pressure and the calculated Henry’s constant: aq vap vap xchain ) ˆf chain ⁄ Hchain ) ˆf chain ⁄ Hchain ≈ Pchain ⁄ Hchain

(4)

Our method assumes that the hydrocarbon chains exist at infinite dilution in aqueous solution, and we make no provision for the possible presence of multimeric aggregates. The excellent agreement of our results with experiment, where reliable data exist (Figure 1), provides a posteriori validation of these assumptions. II.D. Free Energy Landscapes. Histograms in a (multidimensional) order parameter constructed from canonical ensemble simulation data may be used to calculate, within an additive constant, the Helmholtz free energy, using the statistical mechanical identity

A(ψ _ ) ) -kBT ln P(ψ _ ) + constant

(5)

where A(ψ _ ) is the Helmholtz free energy associated with configurations in which the system has an order parameter ψ _ ( ∆ψ _ ; P(ψ _ ) is the probability, as determined from the histograms, of observing the system to exist in a state with order parameter ψ _ ( ∆ψ _. Section V presents one-dimensional Helmholtz free energy landscapes in the chain radius of gyration, Rg, constructed from simulation trajectories of n-alkane and “soft” chains at infinite dilution in aqueous solvent at 298 K and a density corresponding to 1 bar, and in the ideal gas phase at 298 K. The resolution of the landscapes was specified by the histogram bin size, ∆Rg, which was better than 0.0125 Å. The same data were used to construct two-dimensional landscapes in linear combinations of the torsional angles along the n-alkane chains with bin resolution better than 0.08 rad. These results are presented in section IV. II.E. Cavitation Calculations. Free volumes induced by the presence of the n-alkane and “soft” chains in water were

6408

J. Phys. Chem. B, Vol. 113, No. 18, 2009

measured by considering the O atoms of SPC/E water and the CH4/CH3/CH2 united atoms of the chains as hard spheres with radii equal to half the associated Lennard-Jones σ parameter. A three-dimensional cubic mesh was constructed in the simulation box to split the volume into cubic cells of side length 0.2 Å. Further refinement of the mesh was found not to affect results. The solvation cavity of the hydrocarbon was defined as consisting of those cells for which the test insertion of a 2.167 Å spherical probe at the center of the cell resulted in no overlaps with water O atoms. These cells were then subjected to test insertions of a 0.5 Å probe, and those cells for which no overlap occurred with the united atoms of the hydrocarbon were considered part of the cavitation bubble induced by the solvated chain. A single probe approach was found inadequate to provide the necessary resolution to detect cavitation associated with the hydrocarbon while being sufficiently large to exclude the detection of natural cavities in bulk water. The radius of the large probe was specified such that the probability of a successful insertion in bulk water51 was less than 10-7. The qualitative nature of the free volume results was robust to the radius of the small probe over the range 0-0.75 Å. For reference, the insertion probability of the small probe of radius 0.5 Å in bulk water is ∼10%51 and the most probable spherical cavity radius in bulk water is ∼0.25 Å.52 The computational expense of this procedure, combined with the large storage requirements of solvent configuration files, precluded the analysis of every frame of the REMD trajectories. As a result, a subset of 2500 frames evenly spaced throughout the trajectory for each chain length considered were subjected to offline analysis and the results were appropriately weighted to yield a thermally averaged cavity volume. For each frame considered, the bubble volume was computed and binned according to the corresponding radius of gyration of the chain into bins of width 0.25 Å. Averaging over the bins was then performed by weighting the mean bubble volume in each window according to Rg histograms computed over the course of the entire simulation trajectory. This approach returned better results than performing a simple average over the free volumes, since for most chain lengths the sample size was too small for the Rg distributions of the 2500 frame subset to correspond to that computed over the entire trajectory. III. Solubility of n-Alkanes The result of the our solubility calculations described in section II.C is Henry’s constant, Hchain, at 298 K and 1 bar for n-alkane chains ranging from methane (C1) to n-docosane (C22). The vapor pressures required for the calculation of solubilities using eq 4 were calculated as follows. For n-alkanes lighter vap was than n-pentane (C5), which are gases at 298 K,53,54 Pchain taken to be equal to the total pressure, 1 bar. For n-pentane (C5) through n-heptadecane (C17), which are liquids under these vap was computed from standard literature conditions,53,54 Pchain expressions.53 The stable states of the n-alkanes heavier than n-heptadecane (C17) under the same conditions are solids,53,54 but for the present solubility calculations, these species were considered supercooled liquids so as not to introduce a break in the solubility trend due to a phase change. Consequently, vap for n-octadecane (C18), nthe calculated values of Pchain nonadecane (C19), and n-eicosane (C20) were computed using Antoine equation53 extrapolations of 3, 7, and 11 K, respectively, into the supercooled region. Due to a lack of Antoine equation data, the vapor pressures of n-heneicosane (C21) and n-docosane (C22) were calculated using polynomial fits developed by Chickos and Hanshaw,55 and required supercooled extrapolations

Ferguson et al. of 15 and 20 K, respectively. Since the fugacity of the stable solid phase is lower than that of the supercooled liquid, our results represent an upper limit of the true equilibrium solubility. The calculated n-alkane equilibrium solubilities up to ndocosane (C22) are presented alongside data from experiment and a group contribution model in Figure 1. Our solubility calculations are in excellent agreement with experiment for n-alkanes shorter than n-decane (C10), where experimental data from different sources are consistent. Extending up to npentadecane (C15), our results show the best agreement with Tolls et al.23 and do not support the break at n-undecane (C11) observed in some experimental data sets, suggesting that such findings may be due to experimental artifacts. For longer chains, we predict solubilities up to 5 orders of magnitude lower than the existing experimental data, a discrepancy further increased by our treatment of n-octadecane (C18) and heavier chains as supercooled liquids, rendering these results an upper bound on the true solubility. Agreement of our results with the group contribution model of Plyasunov and Shock31 is excellent below n-pentadecane (C15), but for longer chains, our results depart from the exponentially decreasing solubility predicted by this model. The excellent agreement of our results with experiment and the group contribution model for short chains lends credence to our values for the longer n-alkanes in the region where the experimental data are sparse and scattered, and the validity of the group contribution model extrapolation is unknown. The interesting departure from exponentially decreasing solubility for n-pentadecane (C15) and longer chains is better understood by recognizing that eq 4 contains contributions from the vapor pressure and Henry’s constant. Since the decrease in n-alkane vapor pressures is log-linear in carbon number beyond n-butane (C4), the deviation toward higher equilibrium mole fractions is a consequence of our finding that Henry’s constant increases log-linearly with carbon number up to n-tetradecane (C14), but, as shown in Figure 2, breaks from this trend for longer chains and exhibits a maximum at approximately n-octadecane (C18). The vapor pressure and Henry’s constant contributions to the solubility are presented on semilogarithmic axes in Figure S2 in the Supporting Information. One compilation of experimental data20 displays a Henry’s constant maximum around n-dodecane (C12), with this trend reproduced in a qualitative manner by simulations of simplified hydrophobic polymer chains.56 Our data show better agreement with experimental data from another source,28 but since the data do not extend beyond n-dodecane (C12), they cannot be used to corroborate our maximum. The inability of the group contribution model31 to capture the long-chain behavior of Henry’s constant suggests that extrapolation of this model beyond n-tetradecane (C14) is probably unwarranted. Our solubility calculations also furnish values of µex chainswhich are closely related to Henry’s constant via eq 3sand which may be interpreted as the free energy of transfer from the ideal gas to the aqueous phase. Experimental data for, and Pratt-Chandler theory applied to, the transfer of n-alkanes ranging from n-butane (C4) to n-decane (C10) from n-hexane solvent to water,57 and experimental data for the partitioning of n-carboxylic acids up to C21COOH between n-heptane and water,58 both show a linear increase with carbon number of the free energy of transfer from neat hydrocarbon solvent to aqueous phase with a slope of ∼0.85 kcal/mol per carbon atom. As illustrated in Figure 3, we find a corresponding linear trend in the free energy of transfer of n-alkanes from the ideal gas to aqueous phase up to n-tetradecane (C14) with a slope of ∼0.15 kcal/mol per carbon atom, but for longer chains the trend breaks from linearity and

n-Alkane Chains in Water

J. Phys. Chem. B, Vol. 113, No. 18, 2009 6409

Figure 3. Excess chemical potential of n-alkane chains in water at 298 K and 1 bar with respect to ideal gas phase at the same density ex aq IG and temperature, µchain ) µchain - µchain , divided by kBT. By our definition, this quantity is identical to the free energy of transfer from the ideal gas to the aqueous phase. The error bars represent the formal uncertainties computed from the propagation the uncertainties through the incremental Widom insertion calculation.

exhibits a maximum at n-octadecane (C18). Hydrophobicity theory predicts that the transfer free energy of chain molecules from the ideal gas to aqueous phase scales linearly with chain volume (carbon number) up to some critical chain length, beyond which the solvent interaction causes the chain to preferentially adopt collapsed conformations, causing the transfer free energy to scale with surface area, or sublinearly in carbon number.1 Crucially, as demonstrated below, we find no evidence that the observed break from linearity at n-tetradecane (C14) is mediated by the hydrophobic effect; rather, we find that it is driven by ideal gas statistics. IV. Chain Structure We chose to characterize the structure of the n-alkane chains by the radius of gyration, Rg, and the state of the torsional angles. Figure 4a displays Rg as a function of chain length for the n-alkane chains, from which it is apparent that, within measurement uncertainty, the radii of gyration of n-alkane chains up to n-eicosane (C20) do not change upon transfer from the ideal gas to the solvated phase. This result was found also to hold true for each of the principal moments of the gyration tensor. An analogous plot for the “soft” chains is presented in Figure 4b. For C10 chains and shorter, the reduction in Rg upon solvation arises from the unphysical nature of the “soft” potential in which the only interaction between united atoms separated by less than four bonds is the specification of fixed bond lengths. Consequently, we observe artifacts of the model at short chain lengths whereby the chains collapse into highly compact conformations containingoverlapsbetweennon-neighborunitedatomsspresumably to maintain the aqueous hydrogen bonding network1swith little energetic penalty. For chains longer than C10 these artifacts disappear and, similar to our findings for n-alkane chains, the radii of gyration and principal moments of the gyration tensor do not measurably differ between the ideal gas and solvated phases up to C20. Figure 5 shows the relative occupancies of gauche and trans torsional states for n-octadecane (C18) in the ideal gas and solvated phase averaged over all torsional angles in the chain. Consistent with work by Sun et al.,59 we find marked similarity in the torsional angle behavior between the ideal gas and solvated phases. For all chain lengths considered up to neicosane (C20), the fractional occupancies were identical within measurement uncertainty.

Figure 4. Radii of gyration of n-alkane (a) and “soft” (b) chains at 298 K and 1 bar in ideal gas phase (open circles) and in water (filled squares). Ideal gas phase error bars are smaller than the symbol size and were computed from the standard deviations of the averages over five independent MC runs. Solvated phase error bars represent the standard deviation computed via block averaging of the full 1 ns MD simulation trajectory; block lengths were determined by the decay of the chain averaged torsional angle autocorrelation function, and were typically 20-40 ps.

Figure 5. Fraction of samples taken during solvated REMD simulations at 298 K and 1 bar and ideal gas MC simulations at 298 K in which each torsional angle in n-octadecane (C18) was observed in a “low” gauche [0°,120°) (black bar), trans [120°,240°) (gray bar), or “high” gauche [240°,360°) (white bar) state, averaged over all torsional angles in the chain. Error bars were computed in the manner described in the caption to Figure 4.

The n-alkane torsional angles were further analyzed by applying principal components analysis60 to the set of torsional

6410

J. Phys. Chem. B, Vol. 113, No. 18, 2009

angle data extracted from the solvated and ideal gas simulations. The purpose of the analysis was to find the two orthogonal linear combinations of the torsional angles for which the data exhibit the most variance, under the assumption that such combinations represent the important dynamics of the chain. For most chain lengths, the coefficients for the linear combination of torsional angles of the first and second principal components suggested symmetric and antisymmetric combinations respectively of the torsional angles about the midpoint of the chain. Exploiting the inherent symmetry of the chains in an isotropic solvent, the first principal component (PC1) was then manually symmeterized by taking the mean of the coefficients pertaining to torsions i and n - i + 1, where i is the index of the torsional angle in the chain and n is the total number of torsional angles in the chain. The second principal component (PC2) was similarly antisymmeterized by taking the mean of the absolute values of the i and n - i + 1 coefficients and setting the i coefficient to the negative mean and the n - i + 1 coefficient to the positive mean; for chains with an even number of carbon atoms the contribution of the central torsion in PC2 was set to zero. The modified PC1 and PC2, measured in radians, were then used to parametrize two-dimensional free energy landscapes as described in section II.D. For chains up to n-decane (C10) the ideal gas and solvated phase landscapes exhibit multiple local free energy minima corresponding to chain conformations in which the torsional angles occupy gauche (∼60° and ∼300°) and trans (180°) states corresponding to the three minima in the torsional potential of the TraPPE model.34 The similarities in both the location and depth of the free energy minima in the ideal gas and solvated phase landscapes for these chain lengths is striking, providing further evidence that both the conformations exploredsat least in torsional angle phase spacesand the relative populations of such conformations are highly conserved between the ideal gas and solvated phases. The free energy landscapes for n-hexane (C6) are presented in Figure 6. The two axes of symmetry arise from the longitudinal and transverse symmetry of n-alkane chains in an isotropic solvent. The axes have been shifted so that the origin corresponds to the all-trans global free energy minimum. To a good approximation, the local free energy minima correspond to increasing numbers of gauche defects (one, two, and three) with increasing distance from the all-trans well. For chains between n-decane (C10) and n-eicosane (C20) the number of stable torsional states of the chain increases as 3n, where n is the number of torsional angles in the chain. While not all of these states are significantly populated, the density of free energy minima on the landscape vastly increases, leading to merging of these basins into a single global funnel centered on the all-trans conformation. The absence of distinct features, coupled with progressively more incomplete sampling of the conformational phase space for solvated chains of increasing length, makes comparison of ideal gas and solvated phase landscapes difficult for longer chains. Qualitatively, however, the funnel-shaped landscape is conserved upon solvation from the ideal gas. Due to the absence of a torsional potential function, a similar analysis applied to the “soft” chains resulted in relatively featureless funneled landscapes centered on the alltrans conformation for all chain lengths up to n-eicosane (C20) in both the ideal gas and solvated phase. These findings suggest that the ideal gas chain statistics of n-alkanes and sufficiently long “soft” chains dominate the direct solvent interaction and dictate the ensemble of conformations explored by these molecules in water. The idea that the structure

Ferguson et al.

Figure 6. Helmholtz free energy landscapes of n-hexane (C6) in the ideal gas at 298 K (a) and solvated phase at 298 K and 1 bar (b) as a function of the first two principal components of the chain torsional angles which have been modified into symmetric (PC1) and antisymmetric (PC2) combinations of the torsional angles. A is the Helmholtz free energy, and β ) 1/(kBT).

of n-alkane chains in water is an extended one has previously been suggested by molecular simulation of n-dodecane (C12) by Wallqvist and Covell12 and AM1-SM2 free energy calculations up to n-pentadecane (C15) conducted by Tolls et al.23 Monte Carlo simulations of n-octadecane (C18) by Sun et al.59 led these authors to further suggest that conformational differences between the ideal gas and solvated phase are “statistically significant but not dramatic.” In section V we determine more precisely the nature of the solvent interaction. V. Radius of Gyration Free Energy Landscapes While the mean radii of gyration of n-alkane chains shorter than n-eicosane (C20) are unchanged upon solvation from the ideal gas, free energy landscapes parametrized by Rg were constructed to determine the distribution of the conformational ensemble in this parameter. This parametrization provides a good characterization of the degree of compactness of the chains and provides a thermodynamically useful one-dimensional projection of the free energy hypersurface. Transition path sampling applied to simulations of a hard sphere polymer in coarse-grained solvent has, however, demonstrated that the radius of gyration alone is insufficient to adequately capture the dynamics of chain collapse and is not by itself a kinetically relevant reaction coordinate.2 Landscapes for ideal gas and solvated n-alkanes ranging from n-butane (C4) to n-eicosane (C20) are presented in Figure 7a, from which it is clear that the major features of the solvated landscape are well described by ideal gas chain statistics. The two-well landscape for both the ideal gas and solvated n-butane (C4) chains arises from geometric considerations: the low Rg basin corresponds to chain conformations in which the

n-Alkane Chains in Water

Figure 7. Helmholtz free energy landscapes for n-alkane and “soft” chains in the ideal gas at 298 K (thin black line) and solvated phase at 298 K and 1 bar (thick pink line). A is the Helmholtz free energy, β ) 1/(kBT), and Rg is the radius of gyration. In each case, the various curves were shifted in the y-direction by specifying the arbitrary constant in eq 5. (a) Landscapes for (bottom to top) n-butane (C4), n-hexane (C6), n-octane (C8), n-decane (C10), n-dodecane (C12), n-tetradecane (C14), n-hexadecane (C16), n-octadecane (C18), and n-eicosane (C20). (b) Landscapes for “soft” chains (bottom to top) C4, C6, C8, C10, C12, C14, C16, C18, and C20. (c) Ideal gas landscapes for (bottom to top) n-docosane (C22), n-tetracosane (C24), nhexacosane (C26), n-octacosane (C28), n-triacontane (C30), n-pentatriacontane (C35), n-tetracontane (C40), and n-pentacontane (C50).

J. Phys. Chem. B, Vol. 113, No. 18, 2009 6411 single torsional angle is in one of the two gauche states, and the high Rg basin corresponds to the trans conformation. The free energy barrier between the wells corresponds to transitions to and from the all-trans conformation. For longer chains, the all-trans state is separated by this barrier from torsional angle conformations with a single gauche defect in the first or last torsional angle (except in the case of n-hexane (C6) in which configurations with a single gauche defect in the central torsion have the same value of Rg as those with a head or tail defect). With an increasing number of torsional angles, the number of paths by which the chain may transition to and from the alltrans state increases, leading to reduction in the apparent free energy barrier height with chain length in this one-dimensional free energy projection. This all-trans barrier is discernible up to n-hexadecane (C16), for which a weak signature can be seen at Rg ≈ 5.7 Å, but is obscured by noise for longer chains. Both the ideal gas and solvated landscapes exhibit increasing stabilization of low Rg configurations with increasing chain length, resulting in the development of a distinct free energy basin by n-eicosane (C20) in the ideal gas, and by n-dodecane (C12) in the solvated phase. For n-dodecane (C12) chains and longer, the solvent interaction causes the appearance of a free energy barrier of height ∼kBT, measured from the compact basin. For n-hexadecane (C16) chains and longer, the effect of the solvent is additionally manifested in the destabilization of the most extended chain conformations compared to the ideal gas. Analogous free energy landscapes for chains governed by the “soft” potential are presented in Figure 7b. As described above, we observe artifacts of the unphysical potential model at shorter chain lengths where the chains collapse into highly compact conformations containing atomic overlaps, but for C10-C20 we find trends similar to those observed for the n-alkane chains, although the all-trans barrier is not observed due to the absence of a torsional potential. Specifically, we find remarkable similarity between the ideal gas and solvated phase landscapes with a solvent-induced free energy barrier of height ∼kBT, measured from the compact basin, observed for C14 chains and longer, and destabilization of the most extended chain conformations. Increasing stabilization of compact relative to extended conformations with increasing chain length precipitates a transition to a low Rg global free energy minimum at C18 in both the ideal gas and solvated phases. Switching off the bending and torsional potentials reduces the chain stiffness and makes accessible those chain configurations containing steric overlaps between proximate united atoms at no energetic cost. As a result, both the energetic cost and conformational entropy loss of collapse are diminished,6 causing a low Rg global free energy minimum to be adopted at relatively shorter chain lengths for the “soft” chains. Additional ideal gas MC simulations of n-alkanes up to n-pentacontane (C50) (Figure 7c) show that the global free energy minimum shifts to the low Rg basin for chains longer than n-hexacosane (C26). In the case of the “soft” chains, this transition occurred simultaneously in the ideal gas and solvated phases at C18, suggesting a similar simultaneous transition at n-hexacosane (C26) for the chains modeled with the full TraPPE potential. For n-tetracontane (C40) and n-pentacontane (C50), the ideal gas free energy landscapes become essentially single welled due to an overwhelming preference for compact chain conformations. The ideal gas free energy landscape determined for ndodecane (C12) is very similar to that determined by Wallqvist and Covell,12 but while we found remarkable similarity between the ideal gas and solvated phase landscapes, these authors found

6412

J. Phys. Chem. B, Vol. 113, No. 18, 2009

solvation to slightly destabilize extended and significantly stabilize compact chain conformations, leading them to conclude that “the hydrophobic effect is present and of great importance in the compaction of [n-dodecane] hydrocarbon chains.” Given that sampling difficulties arising from structural hysteresis necessitated the use of REMD over 1 ns trajectories to extract converged reliable thermodynamic averages, we speculate that the short 25 ps canonical MD simulations at each point in the perturbation pathway, orsas observed by the authors themselves12sthe artificial modification of the chain dynamics by this approach, may have prevented the sampling of true equilibrium free energies. Interestingly, although these authors demonstrated an approximately linear relationship between the head-to-tail distance and the radius of gyration of the chain,12 we find that plotting the free energy as a function of the former results in noisier landscapes for which the solvent-induced free energy barrier is not nearly as apparent as those plotted as a function of the latter. This illustrates the importance of selecting good order parameters to reveal the essential features of the data.61 Two-well free energy landscapes parametrized by Rg have also been observed for hydrophobic polymers in water by ten Wolde and Chandler2 and Athawale et al.4 While the latter study focused on the importance of the hydration contributionswhich the authors distinguish from the overall solvent interactionsin promoting polymer collapse, reproductions of the free energy landscapes for a 25-mer composed of methane-like monomers (Figure S3a in the Supporting Information) demonstrate consistency with the present work in that that the two-basin topology of the ideal gas landscape persists in the solvated phase, with the solvent interaction resulting in the appearance of a free energy barrier separating the compact and extended basins. For this potential model, solvation also results in destabilization of the compact basin relative to the extended one by ∼3kBT, quite contrary to what one would expect from chains subject to the hydrophobic effect. Athawale et al.4 also considered a 25-mer hydrophobic polymer constituted of ethane-like monomers (Figure S3b in the Supporting Information). Consistent with our results for n-pentacontane (C50) in the ideal gas phase (Figure 7c), the landscape is single-welled. The similarity of the landscapes in the ideal gas and solvated phase is striking, with the solvent effect manifest as an elevation of the free energy at intermediate Rg values. Simulations of n-octadecane (C18) by Sun et al.59 also demonstrated marked similarity in the ideal gas and solvated free energy landscapes, but did not find evidence of a two-well structure, possibly due to the relatively coarse resolution of the histograms. These authors also demonstrated very little disruption of the water hydrogen bonding network around n-octadecane (C18),59 which one would expect to find for a hydrophobically collapsed chain.1 The structural hysteresis between extended and compact chain conformations on nanosecond time scales in canonical MD simulations described at the beginning of section II leads us to conjecture that, in addition to the use of an n-alkane potential model with an excessive preference for collapsed conformations,59 the overwhelming preference for collapsed chain conformations observed in MD simulations of n-eicosane (C20) by Mountain and Thirumalai5 may also be attributable to kinetic trapping of a collapsed state. It should be noted that this work was primarily concerned with the denaturation mechanism of urea, in which the precise n-alkane potential model employed was not of primary import.

Ferguson et al.

Figure 8. Mean solute induced cavity volumes for n-alkanes (open circles) and “soft” chains (filled squares) in water at 298 K and 1 bar. Error bars represent the standard deviation computed via block averaging the cavity volume calculation into three equally sized blocks.

Free energy landscapes parametrized by the chain radius of gyration in the ideal gas and solvated phases show remarkable similarity, indicating that the behavior of n-alkane chains in water is governed primarily by ideal gas statistics. While the dominance of ideal gas behavior over the solvent contribution of solvated n-alkanes has been previously suggested,23,59 we have determined here that the effect of the solvent interaction is the appearance of a free energy barrier of order kBT separating the compact and extended free energy basins for sufficiently long chains, and a destabilization of the most extended chain conformations. The lack of a solvent-induced stabilization of the compact relative to the extended free energy basin indicates the absence of hydrophobic collapse for n-alkane chains up to n-eicosane (C20). That this finding is robust to a nontrivial modification of the potential function in which bending and torsional contributions are turned off suggests that the absence of strong solvation effects is fundamental to relatively short chains (e20-mer) composed of small hydrophobic monomers. We do expect the hydrophobic effect to become important for longer chains, but the point at which this occurs remains to be determined. Furthermore, in agreement with previous work on n-dodecane (C12),12 our free energy landscapes indicate that the most probable chain conformation for n-alkanes up to n-eicosane (C20) in length is an extended one. However, we do not find the alltrans conformation to be a good representative for either the ideal gas or solvated chains23,59 for chains longer than around n-decane (C10), where the free energy of the all-trans conformation is ∼2kBT or more less stable than the most probable chain conformation. Indeed, for n-hexadecane (C16) chains and longer, the all-trans free energy is higher than that of the compact free energy well. VI. Solute-Induced Cavitation A close association between cavitation and chain collapse was previously demonstrated by ten Wolde and Chandler2 for a simplified model of a hydrophobic polymer comprising 12 hard spheres in a coarse-grained solvent. These authors showed chain collapse to proceed by the growth of a vapor bubble into which the polymer collapses, and the nucleation of the bubble to be associated with a free energy barrier of several kBT. Cavitation volumes in our simulations of n-alkane and “soft” chains were computed as described in section II.E and found to increase approximately exponentially with carbon number as illustrated in Figure 8. Visualization of the chain conforma-

n-Alkane Chains in Water

J. Phys. Chem. B, Vol. 113, No. 18, 2009 6413 conformations and the development of the solvent-induced free energy barrier. A possible mechanism is that the increasing population of compact chain conformations in water with increasingchainlengthsasdictatedbyidealgasstatisticsspromotes cavitation in the cores of these conformations, and the observed barrier separating the compact and extended free energy wells (Figure 7a,b) results from the cost of cavity formation. The aforementioned inadequacy of the radius of gyration as a kinetically relevant reaction coordinate, however, precludes a kinetic interpretation of this barrier. VII. Summary and Conclusions

Figure 9. Illustration of cavities (green volumes) induced by n-alkane chains in water at 298 K and 1 bar. The n-alkane united atoms have radii equal to half the corresponding Lennard-Jones σ parameter. The cubic cells of side length 0.2 Å defining the cavity were expressed as spheres of radius 0.1 Å and the cavity rendered using the “Surf” drawing method in VMD62 with a probe radius of 2.0 Å. The water solvent has been removed for clarity. (a) A compact n-octadecane (C18) chain conformation with Rg ) 3.8 Å supports a single very large cavity of volume 29 Å3 in its core; such large cavity volumes were rarely observed. (b) Two small cavities totaling 0.09 Å3 situated near kinks in an extended n-eicosane (C20) chain with Rg ) 5.7 Å.

tions with the cavity volumes superimposed suggests that the larger cavities exceeding 1 Å3 reside in the cores of the compact chains, while smaller volume cavities tend to exist around kinks in extended chain conformations. Figure 9 displays snapshots of cavities induced by n-octadecane (C18) and n-eicosane (C20) chains rendered using VMD.62 The most probable radius of a spherical cavity in bulk water, where the O atoms are treated as hard spheres of radius 1.35 Å, is around 0.25 Å,51 with a corresponding volume of 0.065 Å3. Treating this as a characteristic value, the cavitation volume associated with both n-alkane and “soft” chains breach this threshold for chains comprising between 12 and 14 carbon atoms, precisely the chain lengths at which the solvent-induced free energy barrier first appears and at which the transfer free energy of n-alkanes from the ideal gas to aqueous phase breaks from linearity in carbon number. This suggests a connection between the onset of cavitation in the cores of the compact chain

We have conducted extensive REMD and MC simulations of n-alkanes combined with an incremental Widom insertion technique to determine their solubilities in water at 298 K and 1 bar, as well as their structural properties as a function of chain length. Our calculated solubilities are in excellent agreement with experiment below n-decane (C10) and in good agreement with the experimental data of Tolls et al.23 and the group contribution model of Plyasunov and Shock31 up to n-pentadecane (C15). For longer chains, our results do not corroborate the limited experimental data, but neither do we find agreement with a simple log-linear extrapolation of the short-chain results as predicted by the group contribution model. We have determined that there is a maximum in Henry’s constant as a function of chain length around n-octadecane (C18) (Figure 2) which induces some curvature to the solubility trend beyond n-tetradecane (C14) leading to higher solubilities than would be expected from extrapolation of the short-chain results (Figure 1). Our findings do not support experimental data sets that find a sharp break in the trend at n-undecane (C11). Structural analysis of the torsional angles and radii of gyration of n-alkane and “soft” chains, in which the bending and torsional potentials were switched off, demonstrate remarkable similarities in the conformational ensemble in the ideal gas and solvated phases (Figures 4, 5, and 7). This was previously suggested by work on n-octadecane (C18) by Sun et al.59 Furthermore, this result is robust to the elimination of chain stiffness and even to the use of quite different potential models, suggesting that the absence of strong solvation effects on the free energy landscapes is fundamental to relatively short chains (∼20-mers or shorter) composed of small hydrophobic monomers, and does not depend on the precise nature of the chain interactions. While there is no question that the presence of water has a detectable impact on the free energy landscape of a solvated chain relative to the ideal gas state, its effect is relatively subtle. The primary influence of the solvent for both the n-alkanes and their “soft” analogues is the appearance of a free energy barrier of order kBT separating the extended and compact basins on free energy landscapes parametrized by the radius of gyration of the chains, and a destabilization of the most extended chain conformations. Nevertheless, the absence of solvent-induced stabilization of the compact relative to extended chain conformations provides no evidence for hydrophobic collapse of n-alkanes or “soft” chains up to 20 carbon atoms in length. We anticipate that longer chains will be subject to hydrophobic collapse, but the chain length at which this becomes apparent remains to be determined. We have demonstrated the presence of solute-induced cavities for solvated n-alkane and “soft” chains, and shown that the average size of the cavity increases approximately exponentially with chain length. It also appears that the degree of cavitation becomes significant at the same chain lengths at which the solvent-induced free energy barrier appears and, in the case of

6414

J. Phys. Chem. B, Vol. 113, No. 18, 2009

n-alkanes, the transfer free energy from the ideal gas to aqueous phase first breaks from linearity in carbon number. Our work suggests that the driving forces responsible for the collapse of hydrophobic polymers and short (e25 residue) peptides with a hydrophobic core63-65 may be governed more by the conformational statistics of the isolated chain than by hydrophobic collapse. In the case of peptides, it seems likely that specific intrachain interactions, such as salt bridges and hydrogen bonding, play an important role in driving the adoption of folded conformations. Acknowledgment. Financial support for this work was provided by the Department of Energy (Basic Energy Sciences Grant DE-FG02-01ER15121 to A.Z.P.), the Princeton Center for Complex Materials (NSF DMR Grant 0213706 to A.Z.P.) and the National Science Foundation (Collaborative Research in Chemistry Grant CHE0404699 to P.G.D.). We gratefully acknowledge the TIGRESS high performance computer center at Princeton University which is jointly supported by the Princeton Institute for Computational Science and Engineering and the Princeton University Office of Information Technology. Supporting Information Available: Figure S1, showing the enhancement in sampling of the conformational phase space of n-hexadecane (C16) using REMD instead of canonical MD. Figure S2, showing the vapor pressure and calculated Henry’s constants of n-alkane chains up to n-docosane (C22) on loglinear axes. Figure S3, reproducing free energy landscapes for hydrophobic polymers studied by Athawale et al.4 This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Chandler, D. Nature (London) 2005, 437, 640–647. (2) ten Wolde, P. R.; Chandler, D. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 6539–6543. (3) Widom, B.; Bhimalapuram, P.; Koga, K. Phys. Chem. Chem. Phys. 2003, 5, 3085–3093. (4) Athawale, M. V.; Goel, G.; Ghosh, T.; Truskett, T. M.; Garde, S. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 733–738. (5) Mountain, R. D.; Thirumalai, D. J. Am. Chem. Soc. 2003, 125, 1950–1957. (6) Dill, K. A. Biochemistry 1990, 29, 7133–7155. (7) Sharp, K. A. Curr. Opin. Struct. Biol. 1991, 1, 171–174. (8) Lapidus, L. J.; Yao, S.; McGarrity, K. S.; Hertzog, D. E.; Tubman, E.; Bakajin, O. Biophys. J. 2007, 93, 218–224. (9) Walther, K. A.; Grater, F.; Dougan, L.; Badilla, C. L.; Berne, B. J.; Fernandez, J. M. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 7916–7921. (10) Mountain, R. D.; Thirumalai, D. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 8436–8440. (11) Wallqvist, A.; Covell, D. G. Biophys. J. 1996, 71, 600–608. (12) Wallqvist, A.; Covell, D. G. J. Phys. Chem. 1995, 99, 13118–13125. (13) Ashbaugh, H. S.; Garde, S.; Hummer, G.; Kaler, E. W.; Paulaitis, M. E. Biophys. J. 1999, 77, 645–654. (14) Boulougouris, G. C.; Errington, J. R.; Economou, I. G.; Panagiotopoulos, A. Z.; Theodorou, D. N. J. Phys. Chem. B 2000, 104, 4958–4963. (15) Jorgensen, W. L. J. Chem. Phys. 1982, 77, 5757–5765. (16) Economou, I. G. Fluid Phase Equilib. 2001, 183-4, 259–269. (17) Yezdimer, E. M.; Chialvo, A. A.; Cummings, P. T. Fluid Phase Equilib. 2001, 183-184, 289–294. (18) Errington, J. R.; Boulougouris, G. C.; Economou, I. G.; Panagiotopoulos, A. Z.; Theodorou, D. N. J. Phys. Chem. B 1998, 102, 8865–8873. (19) Ma¸czyn´ski, A.; Go´ral, M.; Win´iewska-Gocłlowska, B.; Skrzecz, A.; Shaw, D. Monatsh. Chem./Chem. Mon. 2003, 134, 633–653. (20) Sander, R. Henry’s Law Constants. In NIST Chemistry WebBook, NIST Standard Reference Database No. 69; Linstrom, P. J., Mallard, W. G., Eds.; National Institute of Standards and Technology: Gaithersburg, MD, 2005; http://webbook.nist.gov.

Ferguson et al. (21) Khadikar, P. V.; Mandloi, D.; Bajaj, A. V.; Joshi, S. Bioorg. Med. Chem. Lett. 2003, 13, 419–422. (22) McAuliffe, C. J. Phys. Chem. 1966, 70, 1267–1275. (23) Tolls, J.; van Dijk, J.; Verbruggen, E. J. M.; Hermens, J. L. M.; Loeprecht, B.; Schuurmann, G. J. Phys. Chem. A 2002, 106, 2760–2765. (24) Tsonopoulos, C. Fluid Phase Equilib. 1999, 156, 21–33. (25) Baker, E. G. Science 1959, 129, 871–874. (26) Franks, F. Nature (London) 1966, 210, 87–88. (27) Sutton, C.; Calder, J. A. EnViron. Sci. Technol. 1974, 8, 654–657. (28) Mackay, D.; Shiu, W. Y. J. Phys. Chem. Ref. Data 1981, 10, 1175– 1199. (29) McAuliffe, C. Science 1969, 163, 478–479. (30) Voutsas, E. C.; Tassios, D. P. Ind. Eng. Chem. Res. 1997, 36, 4973– 4976. (31) Plyasunov, A. V.; Shock, E. L. Geochim. Cosmochim. Acta 2000, 64, 439–468. (32) Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 1999, 314, 141–151. (33) Kumar, S. K.; Szleifer, I.; Panagiotopoulos, A. Z. Phys. ReV. Lett. 1991, 66, 2935–2938. (34) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 2569– 2577. (35) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269–6271. (36) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. J. Chem. Phys. 1982, 76, 637–649. (37) Nose´, S. J. Chem. Phys. 1984, 81, 511–519. (38) Hoover, W. G. Phys. ReV. A 1985, 31, 1695–1697. (39) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384–2393. (40) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117. (41) Andersen, H. C. J. Comput. Chem. 1983, 52, 24–34. (42) Ewald, P. P. Ann. Phys. (Leipzig) 1921, 64, 253–287. (43) Nymand, T. M.; Linse, P. J. Chem. Phys. 2000, 112, 6152–6160. (44) Kolafa, J.; Perram, J. W. Mol. Simul. 1992, 9, 351–368. (45) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, 2002. (46) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University Press: New York, 1989. (47) Hukushima, K.; Nemoto, K. J. Phys. Soc. Jpn. 1996, 65, 1604– 1608. (48) Widom, B. J. Chem. Phys. 1963, 39, 2808–2812. (49) de Pablo, J. J.; Laso, M.; Suter, U. W. J. Chem. Phys. 1992, 96, 6157–6162. (50) Sandler, S. I. Chemical, Biochemical, and Engineering Thermodynamics; Wiley: New York, 2006. (51) Hummer, G.; Garde, S.; Garcia, A. E.; Paulaitis, M. E.; Pratt, L. R. J. Phys. Chem. B 1998, 102, 10469–10482. (52) Southall, N. T.; Dill, K. A.; Haymet, A. D. J. J. Phys. Chem. B 2002, 106, 521–533. (53) Perry, R. H.; Green, D. W. Perry’s Chemical Engineers’ Handbook; McGraw-Hill: New York, 1997. (54) Thermodynamics Research Center, NIST Boulder Laboratories, M. Frenkel, Director. Thermodynamics Source Database. In NIST Chemistry WebBook, NIST Standard Reference Database Number 69, eds Linstrom P. J., Mallard W. G., Eds.; National Institute of Standards and Technology: Gaithersburg MD, 2005; http://webbook.nist.gov. (55) Chickos, J. S.; Hanshaw, W. J. Chem. Eng. Data 2004, 49, 77–85. (56) Stone, M. T. Monte Carlo approaches to the protein folding problem; University of Texas: Austin, 2002. (57) Pratt, L. R.; Chandler, D. J. Chem. Phys. 1977, 67, 3683–3704. (58) Smith, R.; Tanford, C. Proc. Natl. Acad. Sci. U.S.A. 1973, 70, 289– 293. (59) Sun, L.; Siepmann, J. I.; Schure, M. R. J. Phys. Chem. B 2006, 110, 10519–10525. (60) Garcı´a, A. E.; Sanbonmatsu, K. Y. Proteins: Struct., Funct., Genet. 2001, 42, 345–354. (61) Bolhuis, P. G. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 12129– 12134. (62) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33–38. (63) Struthers, M. D.; Cheng, R. P.; Imperiali, B. Science 1996, 271, 342–345. (64) Rana, S.; Kundu, B.; Durani, S. Chem. Commun. 2005, 2, 207– 209. (65) Neidigh, J. W.; Fesinmeyer, R. M.; Andersen, N. H. Nat. Struct. Biol. 2002, 9, 425–430.

JP811229Q