Simulated Solvation of Organic Ions: Protonated Methylamines in

Article pubs.acs.org/JPCB

Simulated Solvation of Organic Ions: Protonated Methylamines in Water Nanodroplets. Convergence toward Bulk Properties and the Absolute Proton Solvation Enthalpy Céline Houriez,*,† Michael Meot-Ner (Mautner),‡,§ and Michel Masella∥ †

MINES ParisTech, Centre Thermodynamique des Procédés (CTP), 60 bd Saint-Michel, 75006 Paris, France Department of Chemistry, Virginia Commonwealth University, Richmond, Virginia 23284-2006, United States § Department of Chemistry, University of Canterbury, Christchurch, New Zealand 8001 ∥ Laboratoire de Chimie du Vivant, Service d’ingénierie moléculaire des protéines, Institut de biologie et de technologies de Saclay, CEA Saclay, F-91191 Gif sur Yvette Cedex, France ‡

S Supporting Information *

ABSTRACT: We applied an alternative, purely theoretical route to estimate thermodynamical properties of organic ions in bulk solution. The method performs a large ensemble of simulations of ions solvated in water nanodroplets of different sizes, using a polarizable molecular dynamics approach. We consider protonated ammonia and methylamines, and K+ for comparison, solvated in droplets of 50−1000 water molecules. The parameters of the model are assigned from high level quantum computations of small clusters. All the bulk phase results extrapolated from droplet simulations match, and confirm independently, the relative and absolute experiment-based ion solvation energies. Without using experiment-based parameters or assumptions, the results confirm independently the solvation enthalpy of the proton, as −270.3 ± 1.1 kcal mol−1. The calculated relative solvation enthalpies of these ions are constant from small water clusters, where only the ionic headgroups are solvated, up to bulk solution. This agrees with experimental thermochemistry, that the relative solvation energies of alkylammonium ions by only four H2O molecules reproduce the relative bulk solvation energies, although the small clusters lack major bulk solvation factors. The droplet results also show a slow convergence of ion solvation properties toward their bulk limit, and predict that the stepwise solvation enthalpies of ion/water droplets are very close to those of pure neutral water droplets already after 50 water molecules. Both the ionic and neutral clusters approach the bulk condensation energy very gradually up to 10 000 water molecules, consistent with the macroscopic liquid drop model for pure water droplets. Compared to standard computational methods based on infinite periodic systems, our protocol represents a new purely theoretical approach to investigate the solvation properties of ions. It is applicable to the solvation of organic ions, which are pivotal in environmental, industrial, and biophysical chemistry but have been little investigated theoretically up to the present.

I. INTRODUCTION The present paper concerns molecular modeling of the solvation of organic ions that are common in biology, industry, and natural environments. Organic ions, such as protonated amines, are central to acid/base chemistry as it applies, for example, in biology, such as the contribution of basic groups of amino acids to protein structure and solubility, and enzyme catalysis.1 In solution and biology, solvation affects the basicities of these ionic groups as much as the molecular structure itself.2,3 It is hard to unravel the interplay of microscopic forces that affect the solvation of these complex ions, but experimental gas phase ion chemistry can now identify the purely molecular structural effects on acidities and basicities and the interaction with the inner shell solvent molecules. The gas phase data can yield, through thermochemical cycles, the absolute solvation enthalpies ΔHg→aq of the ions. Comparison of the gas phase © 2014 American Chemical Society

and solution properties helps to quantify the major effects of solvation on organic ions and structural effects such as methylation on ion solvation properties. As an example, we may consider the protonated alkylamines, whose experimentbased data are summarized in Table 1.3 Similar to their pKb, the protonation enthalpies in solution of methylamines vary irregularly, over a small range of 4.4 kcal mol−1 with increasing methylation. In comparison, protonation enthalpies in the gas phase become regularly more negative with increasing methylation from NH4+ to (CH3)3NH+, varying over a much wider range of 22.8 kcal mol−1. The range of protonation enthalpies is compressed in solution because, compensatingly, the enthalpies of solvation of the protonated alkylamines Received: February 14, 2014 Revised: May 7, 2014 Published: May 9, 2014 6222

dx.doi.org/10.1021/jp501630q | J. Phys. Chem. B 2014, 118, 6222−6233

The Journal of Physical Chemistry B

Article

Table 1. Thermochemistry of Protonation and Solvation of Methylaminesa

a

+

base B

pKb

ΔHprot aq

ΔHprot g

ΔHBg→aq

ΔHBH g→aq

+

ΔH0,4

/w4 ΔHBH g→aq

ΔHc

ΔHIHB

ΔHhydro

NH3 CH3NH3 (CH3)NH (CH3)N

4.75 3.34 3.27 4.19

−12.6 −13.3 −12.1 −8.9

−204.0 −214.9 −222.2 −226.8

−8.0 −10.6 −12.7 −12.7

−88.6 −81.0 −74.6 −66.8

−62.5 −54.0 −49.4 −44.3

−68.1 −69.0 −67.2 −64.5

7.0 10.0 12.6 20.8

−25.2 −25.6 −20.2 −12.6

0.0 −4.9 −12.0 −17.9

All enthalpy data in kcal mol−1, from ref 3. Absolute ion solvation enthalpies referenced to ΔH0g→aq(H+) = −272 kcal mol−1.60 ΔH0,4: sum of BH+/w

stepwise solvation enthalpies of BH+ by four water molecules. ΔHg→aq 4: solvation enthalpy for transfer of the four-coordinated BH+/(H2O)4 cluster, from gas phase to solution.3 ΔHc (cavity), ΔHIHB (ionic hydrogen bonding), and ΔHhydro (hydrophobic) ion solvation enthalpy components from cluster-based analysis of solvation factors.3

become less negative with increasing methylation, with ΔHg→aq varying over a range of 21.6 kcal mol−1 among these ions. These data revealed the major effects of ion solvation that can compress and even reverse relative intrinsic molecular basicities, as became evident with the first studies of gas-phase basicities and proton affinities,2−15 and similarly for acidities.16,17 Gas-phase studies of solvation of alkylammonium ions in small clusters with four H2O molecules showed that this partial solvation parallels the bulk solvation enthalpies. The cluster binding energies in the BH+/(H2O)4 clusters become regularly and signifiantly less negative by 18.2 kcal mol−1 with increasing methylation from NH4+ to (CH3)3NH+, comparable to the range of bulk solvation enthalpies that decrease by 21.6 kcal mol−1 among these ions (see Table 1). These data show that the further solvation enthalpies for transferring the small clusters BH+/(H2O)4 to solution are nearly constant among the methylammonium ions. In other words, four water molecules reproduce the relative bulk solvation enthalpies of ions, although the small BH+/(H2O)4 clusters lack major bulk solvation terms that vary substantially among the ions (see Table 1). This result is unexpected, but it is well established, as it derives directly from experimental thermochemistry and thermochemical cycles. The similarity of partial and full solvation enthalpies also applies to further alkylammonium ions from NH4+ to (Et)3NH+, as well as to oxonium ions.3,15,18 There seems to be little physical interpretation of the above solvation effects, except some early analysis of solvation terms.7,18 Compared to monatomic ions, the solvation of alkylammonium ions has been much less investigated using theoretical approaches, and most of the reported studies concern NH4+ (see refs 19−22). However, we may quote some recent theoretical results concerning the effects of methylation on the alkylammonium ion solvation, in particular their hydration structure in neat water,23 as well as recent spectroscopic findings concerning small alkylammonium ions interacting in the gas phase with small water droplets (of about 20 water molecules).24 We shall address theoretically some basic questions concerning the solvation of protonated alkylamines, by investigating the interactions of these ions with water in clusters and droplets, that bridge from small clusters toward bulk solvation. For that purpose, we will use purely theoretical accurate polarizable molecular modeling, whose parameters were developed considering only high-level quantum mechanical data in the gas phase for small hydrated ion/water complexes, without any empirical experiment-based parameters. While this is a general model, our aim here is to gain physical understanding for the effects of methylation on ion solvation in small water droplets in terms of entropic and enthalpic effects. This will give also new insights into the transition from clusters to bulk solvation of organic ions, and allow comparison with

predictions by macroscopic-based theories, such as the liquid drop (LD) model,25−29 and with macroscopic cluster-based analysis of solvation factors18 Experimentally, investigating ion/water clusters of various sizes represents also an interesting alternative approach to study condensed-phase properties of charged species, from the hydrated electron up to hydrated trivalent ions such as Eu(III) and Ce(III) (see, among others, refs 30−34). We will show that simulation of ions solvated in water nanodroplets represents also an interesting purely theoretical approach to check the consistency of single-ion solvation energies, that have been up to now derived from experimental thermochemical data indirectly, using nonthermodynamic assumptions (for example, that the binding energies of hydrated anion and cation clusters are equal from 4−6 to infinite water molecules). The most common computational approach to simulate a solvation process at the microscopic scale is based on simulating an infinite periodic system. In the particular case of ions, the energetic data extracted from such simulations have to be postprocessed to remove the drawbacks introduced by periodicity (such as the absence of an explicit air/liquid water interface, for instance). However, the reliability of the correction schemes used is still discussed (see the recent discussions of ref 35, for instance). At the difference of the latter computational approaches, the data extracted from droplet simulations do not have to be postprocessed and they allow one to easily extract and explain the physical forces occurring during the solvation process in terms of microscopic interactions, such as polarization or dispersion effects. We will focus our study on the ions NH4+, CH3NH4+, (CH3)2NH2+, and (CH3)3NH+, and for simplicity, we also denote NH4+ as an alkylammonium ion. For comparison, we will also consider the monatomic cation K+. These ions will be investigated when interacting with different water droplets of 50, 100, 300, 600, and 1000 water molecules. The droplet results will allow us to get further insights concerning the convergence of ion solvation properties from clusters to bulk solvation. However, to check the accuracy of our model, we will also investigate the solvation of these ions in bulk water and interacting at the air/liquid water interface using a standard simulation protocol based on simulating an infinite periodic medium. We will compare our bulk simulation results on solvent structure near the ions, and ion solvation free energies, to earlier studies performed using periodic-based theoretical protocols, and to experimental data.

II. THEORETICAL METHODS In the following, N is the total number of atoms considered, M is the total number of atoms in an ion, Nμ is the total number of polarizable atoms, and Nw is the total number of water molecules. All molecular modeling computations were 6223



Article

case of a radial charge distribution ρ(r) ∝ exp(−c × r3) (c is a parameter and r the distance from an atomic center). Lastly, alkylammonium intramolecular degrees of freedom are handled using standard energy terms (they are all included in Urel). All intermolecular energy parameters are assigned to reproduce quantum data concerning monomer properties (like the ion molecular polarizabilities) and ion/water interaction properties, such as the binding energies and geometries of small ion/water clusters including up to four water molecules. The quantum data corresponding to alkylammonium/water systems are derived from computations performed at the MP2 level of theory using medium to extended basis sets, allowing for binding energy extrapolation at the complete basis set (CBS) limit (unpublished data, most of them are provided as Supporting Information). Note that they in turn reproduce experimental cluster binding energies from equilibrium measurements. For K+ /water systems, we considered the quantum data reported by Lee et al.38,39 The details of the ion/ water parameter assignment are provided as Supporting Information. Concerning small ion/water clusters, the agreement between our model and quantum computations for the binding energies is achieved within 0.3 kcal mol−1 on average, and within less than 0.05 Å for the ion/water oxygen distance. We checked also the sensitivity of our results to the model parameters, in particular those corresponding to methane/ water repulsive and dispersion interactions, as well as to the electrostatic charge set of methyl groups. We will briefly discuss the corresponding main results below. B. Molecular Dynamics, Free Energy Computations, Ion/Water Potential of Mean Force, and Ion Location Probability. The MD simulations are performed under ambient conditions according to a standard computational protocol (see the Supporting Information). Simulations of ion/ water droplet systems and at the air/liquid water interface are performed in the NVT ensemble, whereas bulk ion/water systems are simulated in the NPT ensemble. Bulk and air/ liquid water interface systems include about 1000 and 2000 water molecules, respectively. The simulation duration is 10 ns, and the trajectories are sampled each 2.5 ps, for all the simulated systems. Free energy quantities in the aqueous phase are computed from the standard thermodynamic integration (TI) scheme,40 based on 40 intermediate steps. Each step corresponds to a trajectory of 1.25 ns, whose last 1 ns segment is sampled each 100 fs to compute the TI statistical averages. Lastly, we computed the potential of mean force (PMF) corresponding to the interaction of an ion with the droplet center of mass (COM) and of an ion with the simulation cell center (SCC) according to a standard umbrella sampling technique (see the Supporting Information for details). C. Evaporation Phenomena. To minimize the impact of evaporation phenomena along our trajectories, the droplet systems are confined in a spherical cavity, whose radius Rcavity corresponds to the largest droplet COM/water oxygen distance in the droplet to which 12 Å is added. As a water molecule crosses the cavity boundary, it undergoes a reflection from a perfect elastic collision with the cavity wall. Hence, differently from the simulation protocol used by Caleman et al.41 to investigate the behavior of monatomic ions in a small water droplet, we do not consider any additional and physically questionable repulsive energy term centered at the cavity boundary to prevent evaporation. By defining the total number Ns of interacting water molecules in a droplet as the number of water molecules located within a sphere centered at the droplet

performed with our own code POLARIS(MD), that we used recently for modeling neat water.36 A. The Model. To handle water/water interactions, we consider our recent many-body model TCPE/2013 for water,36 which is shown to provide an accurate description of water clusters in the gas phase, of liquid water under various temperature and pressure conditions, and of the air/liquid water interface. It is based on an energy decomposition Uww close to what we consider for modeling ion/water interactions (see below). However, TCPE/2013 includes a short-range anisotropic many-body hydrogen bond (HB) energy term Uhb instead of a pairwise dispersion term Udisp. Note that the energy term Uww is the potential energy of a water cluster, with respect to individual unbound gas phase water molecules. The ion/water interaction energy Uiw is a sum of five components U iw = U rep + U qq ′ + U disp + U pol + U rel

(1)

corresponding to the repulsive, Coulombic, polarization, dispersion, and intramolecular relaxation energies, respectively. The total potential energy U of an ion/water system, i.e., the energy of the reaction nH2O + ion → ion/(H2O)n, is thus U = Uiw + Uww. The above first three terms of Uiw correspond to interactions between different atoms that are not chemically bonded or that are separated at least by two chemical bonds. These sums are denoted by the superscript * in the following, and rij is the distance between the atoms i and j. Urep is taken as an exponential-based repulsive energy term, whereas the Coulombic and dispersion terms are defined in a standard way: U rep =

*

∑ ∑

aij exp( −bijrij) (2)

i = 1, M j = 1, N

U qq ′ =

4πϵ0rij

(3)

6 * ⎛ rij* ⎞ ⎜ ⎟ ∑ ⎜ ⎟ r j = 1, N ⎝ ij ⎠

(4)

i = 1, M j = 1, N

U

disp

qiqj

*

∑ ∑

=−

∑ i = 1, M

The static charges {qi}i=1,N are located on atomic centers. aij, bij, and r*ij are adjustable parameters. Concerning the dispersion, we consider only nitrogen/oxygen and carbon/oxygen intermolecular interactions. The polarization term is based on an induced dipole moment approach. Only the non-hydrogen atoms are polarizable centers. Their isotropic polarizability is denoted by αi, and their induced dipole moment pi obeys *

pi = αi·(Eiq +

∑

Tij·pj) (5)

j = 1, Nμ

Here, Eqi is the electric field generated on the polarizable atom i by the surrounding static charges qj, and Tij is the dipolar tensor. The corresponding polarization energy term is U pol =

1 2

∑ i = 1, Nμ

pi 2 αi

−

∑ i = 1, Nμ

pi ·Eiq −

1 2

*

∑ ∑

piTij pj

i = 1, Nμ j = 1, Nμ

(6) pol

U accounts also for intermolecular short-range damping effects, according to the original ideas of B. T. Thole37 in the 6224



Article

COM and whose radius is Rcavity = 8 Å, the above simulation protocol yields Ns to differ on average by 0.3 molecules from the total number of water molecules Nw along all of our 10 ns trajectories for ion/water systems. For pure water droplets, Ns differs from Nw by no more than one molecule on average. When plotting the mean ion/water interaction energy U̅ iw computed along our umbrella sampling trajectories vs the reference distance rc, the corresponding profile is almost continuous. The U̅ iw estimates are thus almost insensitive to evaporation phenomena. The mean total water/water interaction energies U̅ ww are much more sensitive to evaporation effects. To minimize the uncertainty tied to these phenomena when averaging U̅ ww among the umbrella sampling trajectories, we fit first the values U̅ ww corresponding to a set of simulations for which the values Ns are the largest, by means of a fifth order polynomial function of rc. Then, we extrapolate the values U̅ ww corresponding to the other simulations from this function. Interestingly, the extrapolated values agree within 1 kcal mol−1 with those computed from the raw data of 50 ns simulations (see the Supporting Information). Lastly, the reference values U̅ ww of pure water droplets (computed from a large set of independent 10 ns simulations of pure water droplets of various sizes) are −389, −828, −2670, −5458, and −9225 kcal mol−1, for Nw varying from 50 to 1000, respectively. The uncertainty affecting these values ranges from 2 to 5 kcal mol−1 (see the Supporting Information).

Along the simulations in bulk water, we computed the energies ΔU̅ ww bulk corresponding to the water disorganization induced by the presence of the ion. They are estimated by calculating the difference between the mean water interaction energies computed from hydrated ion simulations and from a neat water simulation. They range from 34.9 (CH3NH3+) to 37.0 (K+) kcal mol−1. These strong destabilization energies originate from the cavity forming within liquid water and from the strong repulsive water/water interactions occurring in the ion first hydration shell (see the discussions provided in ref 36, for instance). The ion/water PMFs at the air/liquid water interface are plotted in Figure 3. The K+ PMF agrees with that reported by Dang and Chang,44 which was also computed from a polarizable approach. It is very flat within the bulk, and it starts to diverge at about 3 Å from the surface. All the alkylammonium ion/water PMFs are close to the K+ one, except that they all present a small barrier within the bulk, at about 5 Å from the surface (the barrier height is at most 0.3 kcal mol−1). Moreover, these PMFs start to diverge at a slightly smaller distance from the surface than for K+. The surface propensities of K+ and NH4+ are thus close, even if the NH4+ one seems to be slightly stronger than that for K+, in agreement with recent sum frequency generation spectroscopy results.45 Together with the experimental estimates reported by Kelly et al.,46 we summarize in Table 2 the differences in the solvation free energies among the ions here considered, computed from our model. With the exception of the reaction (CH3)2NH2+ → (CH3)3NH+, we note a good agreement between our values and the experimental ones, within at most 2 kcal mol−1. In particular, our model predicts very close solvation free energy for K+ and NH4+, as experimentally reported. Concerning the reaction (CH 3 ) 2 NH 2 + → (CH3)3NH+, our theoretical estimate differs slightly more noticeably from experiment, by 3 kcal mol−1. We also computed the K+ absolute ion solvation free energy using a TI approach, whose results were postprocessed according to the standard correction scheme mentioned in the Introduction. It agrees within 2 kcal mol−1 with the recent experimental estimate reported by Kelly and co-workers46 (see the Supporting Information). In all, the above results demonstrate that our approach can provide an accurate description of the solvation of these ions. B. Water Density and Droplet Radius. We computed the water density ρ(r) at the distance r from the droplet COM along all our droplet simulations. Their profiles are insensitive to the nature of the solvated ions, and they are very close to those computed for pure water droplets (for NH4+ droplets, see the plots shown in Figure 4). The water density ρd within the droplet is constant, and the thickness of the air/droplet interface (i.e., the region where the water density decreases from 95 to 5% of the ρd value) is about 4 Å, regardless of the droplet size, in agreement with all the theoretical estimates reported to date concerning the air/liquid water interface.36,47 All our density profiles in droplets can be accurately fitted with the function ρd + ρg ρd − ρg ρ (r ) = + tanh(ξ(r − R d)) (7) 2 2

III. RESULTS AND DISCUSSION In the following discussions, all ion/water distances are measured from the nitrogen nucleus or from the potassium one. A. Ion Solvation Process in Bulk Water and at the Air/ Liquid Water Interface. Ion/water radial distribution functions g(r) computed along simulations of single ions dissolved in bulk water are shown in Figure 1. The NH4+ and

Figure 1. Ion/water oxygen radial distribution functions in the bulk solution. Dashed gray line, K+; dark line, NH4+; blue line, CH3NH3+; green line, (CH3)2)NH2+; orange line, (CH3)3NH+. For alkylammonium ions, the functions are computed from the nitrogen nucleus. Ion coordination numbers (CNs) mentioned in the text are computed by integrating the g(r) functions up to their first minimum.

K+ coordination numbers (CNs) are 4.5 and 7.0, respectively, in good agreement with all the experimental and theoretical results reported to date (see refs 19, 42, and 43 and the references mentioned therein, for instance). When methylation increases, there are fewer protic hydrogen atoms to which water molecules can bond directly. This explains why alkylammonium CN values decrease with increasing methylation (see also Figure 2). The CNs are 3.5, 2.4, and 1.5 from CH3NH3+ to (CH3)3NH+, respectively.

Here, ρg is the water density in the gas phase, ξ is an adjustable parameter, and Rd is interpreted in the present study as the droplet radius. The Rd values are here 6.7, 8.7, 12.6, 16.0, and 19.0 Å for Nw ranging from 50 to 1000. Regardless of the ion 6225



Article

Figure 2. Hydration structures of ions dissolved in Nw = 300 water droplets. These structures correspond to snapshots extracted from MD trajectories where the ions are restrained close to the droplet surface (the rc parameter of the potential Ures is 10 Å). These instantaneous structures are not the mean solvation structures discussed in the text defined from the radial distribution functions shown in Figure 1. In orange, water molecules lying within 5 Å from any atoms of the ions. In red and white, water molecules belonging to the first hydration sphere of K+ or of the NHn moiety, i.e., located at less than 3 Å.

Figure 3. Ion/water PMFs at the air/liquid water interface. The liquid water surface is located at about 29 Å. Black, NH4+; blue, CH3NH3+; (CH3)2NH2+, green; (CH3)3NH+, orange; K+, dashed gray line.

Figure 4. Water density in droplets as a function of the distance to the droplet COM. Gray, Nw = 50; blue, Nw = 100; green, Nw = 300; orange, Nw = 600; red, Nw = 1000. Violet: fit of the water density with the function defined in eq 7. The vertical dashed lines are located at the droplet radii Rd, which may be interpreted as the location of the droplet surface.

Table 2. Differences in Ion Solvation Free Energy between ion A and ion Ba

a

ion A

ion B

model

experiment

K+ NH4+ CH3NH3+ (CH3)2NH2+

NH4+ CH3NH3+ (CH3)2NH2+ CH3)3NH+

−1.5 6.7 5.8 4.5

0.8 8.8 7.8 7.5

Hence, regardless of the droplet size, the water density within a droplet is constant and its value is close to the bulk one. That justifies the assumption made in different implementations of the LD model, which consists of estimating a droplet effective radius from the bulk water density. C. Probability of Ion Locations and Ion/Water PMF at Air/Droplet Interfaces. For a quasi-spherical droplet, the ion/ water PMF relates to the partition function Z

Experiment: data from Kelly et al.46 All data in kcal mol−1.

and the droplet size, the water density ρd matches the bulk water density ρbulk predicted by the water model TCPE/2013 under ambient conditions, within less than 1% (ρbulk = 0.0331 molecules per Å3 36).

Z∝ 6226

∫ r 2 exp(−PMF(r)/kBT) dr

(8)



Article

where r is the distance between the ion and the droplet COM. Here, the volume r2 dr allows one to account for the ion accessible volume when the ion is located in a spherical shell of thickness dr at a distance r from the droplet COM. From the above equation and by considering the component TSgeom(r) = kBT ln(r2), we get readily the relation providing the probability density P(r) of finding an ion at a distance r from the droplet COM P(r ) = exp( −[PMF(r ) − TSgeom (r )]/kBT )/Z

volume at the droplet core resulting from the component TSgeom. Note that the CH3NH3+ and (CH3)2NH2+ PMFs are close to the (CH3)3NH+ one. For Nw = 50, all the PMFs are characterized by a minimum located at 2 Å from the droplet surface (the droplet radius Rd is 6.75 Å for all of the ions). The PMF minimum magnitude ranges from 0.5 (K+ and NH4+) to 2.0 kcal mol−1 (CH3NH3+). For larger and larger droplets, all the PMFs present a flatter and flatter profile within the droplet core, whereas the PMF minimum close to the droplet surface is less and less accented for the alkylammonium ions. For the small ions K+ and NH4+, that PMF minimum even disappears as soon as Nw = 100. By checking the result sensitivity to the model parameters (see the Supporting Information), we showed that the alkylammonium PMFs, and therefore the corresponding probability densities P(r), are sensitive to the magnitude of the parameter used to account for methyl/water dispersion. Slightly altering this parameter may lead to reinforcing or removing the above-mentioned PMF minimum close to the droplet surface. Interestingly, the difference ΔR̅ between the droplet radius Rd and the most probable ion position within a droplet, and its corresponding variance ΔR̅ 2, are both much less sensitive to the magnitude of methyl/water dispersion interactions. For instance, the ΔR̅ values differ at most by 1 Å, and usually by less than 0.5 Å, when computed from all the methyl/water dispersion parameters considered. The quantities ΔR̅ and ΔR̅ 2 are thus reliable indices to quantify the behavior of an ion within a droplet, such as its surface propensity. For Nw ≥ 100, they both increase as Nw increases (see Figure 7), showing

(9)

TSgeom(r) may be interprated as an entropic component, which is responsible for an ion excluded volume centered at the droplet COM (it forces P(r = 0) to be zero). It may even dominate the effects resulting from microscopic interactions, especially when the ion/water PMF is flat within the droplet. In that case, it strongly favors the ion propensity for the droplet surface. Note that a bulk solution in an ideal infinitely large and deep reservoir would not have this factor because the volume of the bulk is infinitely larger than that of the boundary surface, and every cross section layer at any depth has equal volume. That leads in that case to the following probability density P(r) of finding an ion at a distance r from the bulk core P(r ) ∝ exp( −PMF(r )/kBT )

(10)

Geometry entropic terms analogous to TSgeom must be considered for other solutions of limited volume, as solutions in surface films. The PMFs and the corresponding probability densities P(r) for K+, NH4+, and (CH3)3NH+ are plotted as a function of the ion/droplet COM distance r and for the different droplets in Figures 5 and 6. These plots clearly show the ion excluded

Figure 5. Ion/water PMFs in droplets of size Nw, for K+ and NH4+ (left) and (CH3)3NH+ (right). The distance corresponds to the ion/ droplet COM one. K+ data are shown in dashed lines. Black, blue, green, orange, and red lines correspond to Nw = 50, 100, 300, 600, and 1000, respectively.

Figure 7. Mean alkylammonium distance ΔR̅ as a function of the droplet size Nw. The magnitudes of the error bars correspond to the square root of the variance ΔR̅ 2 of the latter distance. The latter quantities are computed according to ΔR̅ = Rd − ∫ rP(r) dr and ΔR̅ 2 = Rd2 − ∫ r2P(r) dr, with r being the ion/droplet COM distance and Rd the droplet radius.

Figure 6. Ion location probability P(r) in droplets of size Nw, for K+ and NH4+ (left) and (CH3)3NH+ (right). Same notations as for Figure 5. The vertical dashed lines are located at the droplet radii Rd.

the progressive transition of the ion behavior from a “cluster” regime (i.e., ion positions are restrained to a small volume at the droplet surface) to a bulk regime (the ions explore freely a large portion of the droplet volume, farther and farther from the droplet surface). The sensitivity of ion/water PMFs at air/ droplet interfaces to the intensity of dispersion effects originates from the abrupt variation of the ion/water dispersion energy at the droplet surface (see the discussions below as well as those provided in ref 48, for instance). Concerning the methyl groups, by systematically checking the sets of structures extracted from our MD trajectories, they are repelled from the droplet core; i.e., they usually lie in between the nitrogen atom and the droplet surface. Moreover, as ΔR̅ values larger than 4 Å are not observed before Nw = 300, 6227



Article

outside the droplet ranges from nc − 2 to nc (here, nc is the ion CN in bulk water). This explains why the mean ion/water energies U̅ iw(r), ̅ and their corresponding energy components, do not converge toward zero outside of the droplet. This also explains why the short-range ion/water repulsive energies U̅ iw rep(r)̅ vary so smoothly when the ions cross the droplet surface. From our simulations, the mean ion/water Coulombic energies U̅ iw qq′(r)̅ vary smoothly, even when the ions cross the droplet surface. For Nw = 300, the largest difference in the −1 for energies U̅ iw qq′(r)̅ between r ̅ = 0 and 20 Å is 10 kcal mol −1 large alkylammonium ions, and only 2 kcal mol for the small ions K+ and NH4+. We get similar results for all of the droplets. Hence, for the latter small ions, it seems that the water droplet reorganizes its structure to maximize the ion/water electrostatic interactions for all the ion locations. This is supported by the plots of the mean water/water energies U̅ ww(r)̅ shown in Figure 8. For K+ and NH4+, U̅ ww(r)̅ becomes linearly less negative as r ̅ increases, up to r ̅ = 16 Å. A different behavior of U̅ ww(r)̅ is observed for the large ions (CH3)2NH2+ and (CH3)3NH+. Their U̅ ww(r)̅ present a weak maximum within the droplet, at about 5 Å from its core for Nw = 300, which may originate from the water HB network destabilization induced by the ion methyl groups. However, the uncertainties affecting water/ water interaction energies in droplets prevent one from further discussion of these data. As for the Coulombic energies, the ion/water polarization iw U̅ iw pol(r)̅ and dispersion U̅ disp(r)̅ energies are constant within the droplet core for all of the ions. However, they start to change a few Å before the droplet surface, with the polarization energies becoming then more and more negative as r ̅ increases, while the contrary is observed for the dispersion. The ion/water dispersion interactions are thus centripetal, favoring ion solvation inside the droplets, whereas ion/water polarization interactions are centrifugal, favoring ion location at the droplet surface. Further, we note that methylation leads to strengthening the ion/water dispersion interactions and to weakening the ion/water polarization effects within the droplet. This was expected, as the ion/water dispersion energy magnitude is an increasing function of the number of methyl groups, whereas these groups prevent a large number of water molecules from interacting at short range with the charged NH headgroup of alkylammonium ions. This also explains why methylation strongly weakens the ion/water Coulombic interactions; see Figure 8. For the small ions K+ and NH4+, their total ion/water energies U̅ iw(r)̅ are almost constant within 1 kcal mol−1 from the droplet core up to a few Å above the droplet surface, whereas their corresponding water/water energies U̅ ww(r)̅ become linearly less negative as the ion/droplet distance increases. These ions are thus driven into the droplet to optimize water/water interactions. This result clearly disagrees with the earlier one reported by Caleman et al. in their study of monatomic ion solvation in a medium size water droplet.41 These authors concluded that monatomic ions are driven in droplets mainly by ion/water electrostatic effects. However, they considered an additional external force acting on the water molecules to prevent evaporation, which, in turn, also prevents droplet structural distortion, and which may thus largely affect their simulation results. For the large ions (CH3)2NH2+ and (CH3)3NH+, the sum iw iw U̅ rep(r)̅ + U̅ iw pol(r)̅ + U̅ disp(r)̅ is almost constant for 0 Å ≤ r ̅ ≤ 16 −1 Å, within 1 kcal mol (see the Supporting Information). Their

methyl groups are thus not fully solvated below that droplet size. This agrees with recent experimental findings concerning small alkylammonium/water systems (Nw = 19−21) derived from infrared photodissociation and blackbody infrared radiative dissociation at 133 K.24 D. Microscopic Interactions within a Droplet. Along each umbrella sampling trajectory, which corresponds to a specific mean ion/droplet COM distance r,̅ we computed the ion/water mean energies U̅ iw(r)̅ and their components U̅ iw rep(r), ̅ iw iw U̅ iw (r ), U (r ), and U (r ), as well as the water/water mean ̅ ̅ qq′ ̅ pol ̅ disp ̅ energies U̅ ww(r). ̅ For ions solvated in a Nw = 300 droplet, these components are plotted as a function of r ̅ in Figure 8 (see the

Figure 8. Ion/water interaction energies U̅ iw, their components U̅ rep, U̅ qq′, U̅ pol, and U̅ disp, and water/water interaction energies U̅ ww (extrapolated according to the protocol presented in section IIC) for NH4+ (black line), CH3NH3+ (blue), (CH3)2NH2+ (green), (CH3)3NH+ (orange), and K+ (dashed gray line) as a function of the ion/droplet COM distance. The droplet size is Nw = 300, and its radius Rd is 12.5 Å.

Supporting Information for the water/water energy components). The root-mean-square deviation (RMSD) of U̅ iw(r)̅ is about 6 kcal mol−1 for all the ion/droplet systems and for all the simulations (and the RMSD of U̅ ww(r)̅ scales as Nw1/2, as expected; see the Supporting Information). Averaging the above-mentioned energies according to the probabilities P(r)̅ provides the mean interaction energies for a given ion/droplet system, like the ion/droplet U̅ iw and the water/water U̅ ww ones. Note that the mean ion/water energy iw iw components U̅ iw rep, U̅ pol, and U̅ disp are already converged with respect to the droplet size as soon as Nw = 300, within 0.5 kcal mol−1. As already shown in earlier studies of ions at the air/liquid water interface,41,44 water molecules jump from the droplet onto the ions at the droplet surface. Along our simulations, the mean number of water molecules still interacting with an ion 6228



Article

ion/water Coulombic energies U̅ iw qq′(r)̅ start to diverge a few Å before the droplet surface, and their water/water energies U̅ ww(r)̅ present a minimum at about 1 Å above the surface. This may result from the desolvation of the alkylammonium hydrophobic methyl groups, which destabilize the water HB network within the droplets. Hence, contrary to the small ions NH4+ and K+, large alkylammonium ions are thus driven into the droplet by more favorable ion/water Coulombic interactions. That conclusion holds for CH3NH3+, even if its U̅ ww(r)̅ profile is much flatter compared to those of the latter two ions. As previously mentioned, the profiles of the probability densities P(r) are sensitive to the model parameters, in particular, to that handling the water/methyl dispersion interaction intensity. However, the overall flat profiles of U̅ iw(r)̅ and U̅ ww(r)̅ lead to mean energies U̅ iw and U̅ ww almost constant within less than 0.5 kcal mol−1, for all the alkylammonium ions and for all the model parameter sets considered (see the Supporting Information). This confirms the reliability of our conclusions concerning the energetics of ion solvation in droplets. Lastly, we may also note that the magnitudes of the ion/water dispersion energy U̅ iw disp for large alkylammonium ions match those of the sum ΔHIHB + ΔHhydro originally computed by one of us from macroscopic clusterbased analysis of solvation factors (see Table 1).18 This demonstrates that these macroscopic cluster-based terms reflect microscopic molecular interactions that are represented in the present molecular simulation parameters. Lastly, we computed the entropic components TS(r)̅ of our PMFs, i.e., the differences between the mean total potential energy U̅ (r)̅ = U̅ iw(r)̅ + U̅ ww(r)̅ and the PMF(r)̅ along our simulations (see the Supporting Information). The main features of our TS(r)̅ profiles are very close to those reported by Caleman et al. in their study dealing with the solvation of monatomic ions in a medium size water droplet.41 In particular, as already shown by the latter authors, TS clearly favors ion solvation inside droplets, contrary to TSgeom. E. Ion/Water Droplet Interaction Energy: Convergence toward Bulk Limit. The mean total ion/water droplet energies U̅ iw and the mean droplet water/water destabilization energies ΔU̅ ww (induced by the ion presence and computed similarly to ΔUww bulk, see section IIIA) are plotted as a function of the droplet size Nw in Figure 9 for all of the ions. Note that the sum U̅ iw + ΔU̅ ww gives the energy corresponding to the reaction (H2O)Nw + ion → ion/(H2O)Nw. The magnitude of the components ΔU̅ ww reflects the strong destabilization of the water droplet HB network by the introduction of the ion. For alkylammonium ions, these components seem to decrease as Nw increases, until they reach the values 30, 26, 24, and 23 kcal mol−1 for the four ions and Nw = 1000 (for K+, the ΔU̅ ww value is very close to the NH4+ one, 32 kcal mol−1). For the small ions K+ and NH4+, these values agree with the bulk ones ΔU̅ ww bulk, within the error bar originating from evaporation phenomena, about 5 kcal mol−1. Note that the mean water dipole moment in droplets is smaller than that in the bulk (still by about 2% for Nw = 1000, see the Supporting Information), which may also explain the difference in the energies ΔU̅ ww between droplets and the bulk. For large alkylammonium ions, the values ΔU̅ ww in droplets differ more noticeably from the bulk one, up to 12 kcal mol−1 for (CH3)3NH+, suggesting that the water HB network may accommodate the presence of hydrophobic groups more easily

Figure 9. Mean ion/water interaction energies U̅ iw (left) and mean water disorganization energies ΔU̅ ww (right) as a function of the droplet size Nw. Same notations as for Figure 5. For U̅ iw, the dashed lines are the results of the fit with power-law function U̅ iw(∞) + ϵNw−1/3. For ΔU̅ ww, the dashed line corresponds to the mean bulk −1 among the ions. value ΔU̅ ww bulk, which differs at most by 1 kcal mol The error bars are only shown for NH4+ and (CH3)3NH+. Their magnitude corresponds to the uncertainty tied to evaporation phenomena.

in droplets than in the bulk. Note that, for the small droplet systems, the components ΔU̅ ww are already smaller for large ions than for small ones. This may originate from the propensity of the ion methyl groups to be at the droplet surface for small systems (see the discussions of section IIIC). They thus disrupt less the water structure of small droplets than when they are inside large droplets or in the bulk. Our ΔU̅ ww and ΔU̅ ww bulk values are about twice as large as the corresponding cavity energy ΔHc computed from a continuum solvation approach (see the data summarized in Table 1). ΔHc is the enthalpic cost for creating an empty cavity within neat water. It does not account for the strong repulsive water/water interactions occurring in the ion first hydration shell (see, among others, the discussions provided in ref 36). That explains the difference in magnitude between ΔU̅ ww and ΔHc values. The mean ion/water droplet energies U̅ iw can be accurately fitted with a power-law function of Nw U̅ iw(Nw ) = U̅ iw(∞) + ϵN wγ

(11)

Here, U̅ (∞), ϵ, and γ are three parameters. By imposing no constraints on them, the fitting procedure leads to exponent γ values varying from −0.24 to −0.36, showing the U̅ iw terms to be inversely related to the droplet radius Rd. In other words, U̅ iw(Nw) becomes more negative and closer to the bulk limit with increasing droplet radius, in agreement with LD approaches.27,29 By setting γ to the ideal value proposed in these approaches (γ = −1/3; such a value was also assumed to extrapolate experimental Eu(III)/water droplet data to their bulk limit31), the fitting procedure leads to the U̅ iw(∞) values −122.2, −113.6, −106.4, and −100.3 kcal mol−1, respectively, for NH4+ to (CH3)3NH+ and −119.6 kcal mol−1 for K+. These values are thus the extrapolated ion/water interaction energies in bulk water, and they can be combined with the energy effects of the ions on the bulk water structure, as given by ΔU̅ ww bulk, to calculate the overall single-ion solvation enthalpies ΔHg→aq iw

ww ΔHg → aq = U̅ iw(∞) + ΔU̅bulk − kBT

(12)

This leads to the ΔHg→aq values −87.7, −79.3, −71.4, and −65.9 kcal mol−1 for NH4+ to (CH3)3NH+ and −83.2 kcal mol−1 for K+. For alkylammonium ions, the agreement with the experimental thermochemistry-based values ΔHexp g→aq, i.e., −88.8, −81.0, −74.6, and −66.8 kcal mol−1,3,18 is particularly good, 6229



Article

the ions, although they originate from independent physical forces. Further, similar results between clusters and bulk solvation apply also larger alkyloxonium and alkylammonium ions. These unexpected similarities between relative solvation energies in small clusters and in bulk solution need further theoretical analysis.3,15,18 While the model confirms the observations about relative solvation energies, the agreement of our extrapolated bulk single-ion solvation energies appears to rule out a previous suggestion in the same source,3,15,18 that the accepted single-ion solvation energies may be too small by tens of kcal mol−1. This was based assuming that the stepwise binding energies level off to the bulk water condensation energy already after 4−6 water molecules. Rather, the present results allow another suggestion there, that the apparent discrepancy arises from cumulative small differences between binding energies of the large cluster and bulk water condensation energy up to large clusters.3,15,18 We shall investigate these effects in further studies. F. Stepwise Water Cluster Enthalpies ΔHNw−1,Nw. From the ion/water interaction energy U̅ iw(Nw) and the total water/ water interaction energy U̅ ww(Nw), the stepwise water cluster enthalpies ΔHNw−1,Nw can be expressed as a sum of three terms

with the differences between the values being at most 4 kcal mol−1. The experiment-based absolute single-ion solvation enthalpies are referenced to the proton value ΔHg→aq(H+), which is derived indirectly and may have an uncertainty at least of ±3 kcal mol−1. For example, refs 3 and 18 use an older value of −271.6 kcal mol−1, while the most current value is lower, −274.8 kcal mol−1.49 By considering our absolute solvation energies for the alkylammonium ions and the data of Table 1, our simulation results yield, indirectly, the proton solvation enthalpy to be −270.3 ± 1.1 kcal mol−1, confirming both the accepted single-ion solvation enthalpies and the reference value for solvation of the proton, within a few kcal mol−1. This is significant because experiment-based single-ion solvation enthalpies can be obtained only indirectly using nonthermodynamics assumptions. The computational single-ion solvation energies are purely theoretical based on quantum mechanical calculations for small clusters and intermolecular parameters, and independent of the indirect experiment-based values, and therefore provide a fully independent confirmation of the accepted single-ion solvation energies. We may also note that the differences in the components U̅ iw and ΔU̅ ww, and thus the differences in their sum ΔHg→droplet, among the ions, appear to be constant (within 5 kcal mol−1, i.e., the uncertainty tied to evaporation) as soon as the smallest droplet here considered, thus at least for Nw = 50. We also investigated small alkylammonium ion/(H2O)4 clusters using our modeling approach. The cluster absolute minima were tentatively identified from systematic quenching of 100 ns gas phase simulations at 150 K. From NH4+ to (CH3)3NH+, the binding energies of four H2O molecules to the ions are −65.0, −59.1, −55.7, and −48.3 kcal mol−1, respectively. The differences among the ions are comparable to the differences among the experimental cluster binding enthalpies ΔH0,4 of −62.5, −54.0, −49.4, and −44.3 kcal mol−1 for the four ions, respectively (Table 1). Further, the differences among these experimental or computed partial solvation energies of the ions by four H2O molecules are similar to the differences among their experiment-based bulk solvation energies ΔHg→aq. In other words, the further solvation enthalpies beyond four H2O molecules, as in transferring the BH+(H2O)4 clusters from gas phase to bulk solution, are nearly constant among the present ions, according to both the computed and experimental data. Moreover, our computed data show that the differences in the solvation energies of the ions remain constant for every cluster size, from Nw = 4 to bulk, which has not been shown before because of the lack of information concerning large clusters. Further, the same solvation enthalpies also apply to larger alkylammonium ions, as well as to the oxonium ions.3,15,18 Note that this result is derived from four-solvated cluster enthalpies and from the differences among the bulk solvation enthalpies of the ions, all of which are well established experimental data. As discussed previously, this finding, that as few as four H2O molecules reproduce the relative bulk solvation enthalpies, is surprising because the small BH+(H 2O)4 clusters lack interactions between the alkyl groups and water, and also lack major bulk solvation factors such as dielectric charging, cavity forming, and hydrophobic solvation, each of which terms varies by up to 30 kcal mol−1 among these ions.3,15,18 The present computational results reproduce this relation between solvation in clusters and in bulk solution, giving independent support to the experimental findings. The results then suggest that the various solvation factors beyond the first BH+(H2O)4 molecules vary in a fortuitously compensating manner among

ΔH Nw − 1, Nw = u ̅Niww − 1, Nw + u ̅Nww + kBT w − 1, Nw

(13)

For large values of Nw, the components u̅ can be expressed as u ̅Niww − 1, Nw = U̅ iw(Nw ) − U̅ iw(Nw − 1) = ∂U̅ iw(Nw )/∂Nw (14a)

u ̅Nww = U̅ ww(Nw ) − U̅ ww(Nw − 1) = ∂U̅ ww(Nw )/∂Nw w − 1, Nw (14b)

The above components u̅ are computed by fitting first the quantities U̅ iw(Nw) and U̅ ww(Nw)/Nw with power-law functions, such as the one discussed in section IIIE (for the derivatives of U̅ iw(Nw), with the power-law function exponent set to −1/3). The fitted functions are then used to compute the above derivatives and the corresponding approximated water stepwise enthalpies ΔH̃ Nw−1,Nw. These enthalpies for the five ions are plotted as a function of Nw in Figure 10, together with the stepwise enthalpies ΔH̃ Nw−1,Nw corresponding to pure water droplets (computed from the derivatives of U̅ ww(Nw)/Nw)) and

Figure 10. Water stepwise enthalpy ΔHNw−1,Nw in ion/water droplets. Plain or dashed lines: data derived from eq 14 (same color definition as in the above figures). Horizontal dotted line: water vaporization enthalpy under ambient conditions. Dashed blue line: water stepwise enthalpy corresponding to pure water droplets. Empty diamonds: quantum data concerning small NH4+/(H2O)1−4 clusters (see the Supporting Information). 6230



Article

with the stepwise enthalpies ΔHNw−1,Nw corresponding to small NH4+/(water)Nw clusters (Nw = 1−4), computed from quantum data (for these small clusters, our model results differ from the quantum ones by less than 0.3 kcal mol−1). Note also that the values ΔH̃ Nw−1,Nw are summarized in the Supporting Information for Nw = 2−50. The ΔH̃ Nw−1,Nw values are insensitive to the ion nature as soon as Nw = 50, and they are then almost indistinguishable from those corresponding to pure water droplets. For small hydrated ion clusters, their ΔH̃ Nw−1,Nw values decrease to reach a minimum between Nw of 8 and 15, where they range then from 8.5 to 9.5 kcal mol−1, in agreement with the experimental stepwise enthalpies reported for small hydrated alkylammonium ions.3,15 In all, our stepwise enthalpies ΔH̃ Nw−1,Nw agree remarkably well with those computed from most of the LD models for a monovalent ion, as discussed by Peslherbe et al.27 and Donald and Williams.29 In particular, the convergence of our ΔH̃ Nw−1,Nw toward their bulk limit is slow. For instance, our ΔH̃ Nw−1,Nw reach 98% of their expected value, i.e., the water vaporization enthalpy ΔH0vap, not before Nw ≈ 104 molecules for all the ions (ΔH0vap = 10.5 kcal mol−1 under ambient conditions). The enthalpies ΔH̃ Nw−1,Nw computed from energy component derivatives are largely overestimated compared to their quantum mechanical and experimental values ΔHNw−1,Nw for Nw = 1−2, while both sets of values are in better agreement for Nw = 3−4. However, for intermediate Nw values (≤20), the values ΔH̃ Nw−1,Nw become larger positive from NH4+ to (CH 3 ) 3 NH + , while the opposite is observed for the experimental values ΔHNw−1,Nw for Nw = 1−4. The relations used for computing our ΔH̃ Nw−1,Nw can be considered as reliable not before Nw = 10. Experimentally, the decrease of the stepwise water addition enthalpies to smaller negative values of −8 to −10 kcal mol−1, compared with −10.5 kcal mol−1 for the bulk condensation enthalpy, has been observed in many series of small hydrated ion clusters of 6−8 water molecules.50 However, in the particular case of the hydrated hydronium H3O+/(H2O)Nw series, an experimental study found ΔHNw−1,Nw of a nearly constant value, 10−11 kcal mol−1, for Nw = 14−28,51 and these bulk-like values are not likely to change up to bulk solution. A Monte Carlo simulation also found nearly constant ΔHNw−1,Nw values for this series at these sizes.52 The stepwise binding enthalpies of these large clusters are still uncertain, and no stepwise energies are available experimentally for any systems past Nw = 28. Presently, the agreement between our molecular models and the LD models, that derives from macroscopic water density and surface tension, supports both approaches. However, the thermochemistry of large clusters is critical for bridging from small clusters to bulk solution, and the properties of these large clusters require further studies.

results can therefore test independently experiment-based data, especially single-ion solvation energies that are obtained from experimental data indirectly and with nonthermodynamic assumptions. All our bulk phase results, using parameters based on quantum chemistry, and droplet data extrapolated to the bulk limit, match the available experiment-based data for absolute single-ion solvation enthalpies. There exist several different approaches for modeling ion solvation and analyzing solvation terms. Some are based on macroscopic properties such as the LD models,27,29 and a cluster-based analysis,18 while simulations build macroscopic and thermochemical properties from microscopic molecular models. It is of interest to summarize here briefly the relations between the present molecular model, the macroscopic models, and experiment: (1) The magnitudes of the ion/water droplet dispersion energy U̅ iw disp match the sum ΔHIHB + ΔHhydro computed from the cluster-based analysis that involves macroscopic continuum dielectric charging and cavity surface tension terms. This supports the microscopic basis of these terms. (2) Our water destabilization energies ΔU̅ ww, in droplets and in bulk solution, are larger in magnitude than the surface tension energy in creating the solvent cavity for the ions. This results because the water/water interaction in our model includes both this term and the strong water/ water repulsion in the first hydration shells of the ions (see, among others, the discussions provided in ref 36). (3) The sum of our ion/water and water/water energies extrapolated to bulk solvation reproduces the absolute experiment-based single-ion solvation enthalpies, mostly within better than 5 kcal mol−1, depending on the indirectly derived reference value for proton solvation. (4) Importantly, our results reproduce the experimental differences between the solvation enthalpies of the ions both by four H2O molecules in clusters and in bulk solution. Both the experimental data and the present independent model show that the differences among bulk solvation enthalpies are reproduced by partial solvation by only four H2O molecules, although major continuum and hydrophobic solvation terms are absent in these small clusters. (5) Our model reproduces the experimental and quantum mechanical calculated stepwise solvation energies by one to four H2O molecules. (6) Our MD model reproduces the stepwise solvation enthalpies from the LD model of the ions by up to 10 000 water molecules. In particular, both models predict slow convergence to the bulk water vaporization enthalpy of 10.5 kcal mol−1 through smaller values of 8.5−10.5 kcal mol−1 up to the 10 000 size. However, the few experimental data show convergence of the stepwise solvation enthalpies to bulk water condensation enthalpy already at about 10−14 H2O molecules. (7) Our model thus demonstrates the validity of the assumptions used when developing the LD approach. As far as we know, this has never been shown previously. (8) The calculated stepwise solvation enthalpies of the ammonium ions and neat water converge to a common value independent of the presence and nature of the core ion, already after 50 water molecules. This is important because these data can be used in thermochemical cycles

IV. CONCLUSIONS We investigated theoretically the behavior of alkylammonium ions and of the monatomic cation K+ at the microscopic level, in water droplets and in bulk water, using a polarizable modeling approach, whose parameters were obtained from quantum calculations. The model is purely theoretical, and its 6231



Article

use of the Stampede supercomputing system at TACC/UT Austin, funded by NSF award OCI-1134872.

to derive single-ion solvation energies. The agreement of our extrapolated values with accepted single-ion solvation energies rules out a cluster-based suggestion that the accepted single-ion solvation energies are too small by tens of kcal mol−1 but allows alternative suggestions about the discrepancy, that the cluster binding energies are different from bulk values up to large clusters.3,15,18 Lastly, our model predicts properties of NH4+ and K+ solvated in water droplets or in the aqueous phase to be close, as experimentally reported. We have to remember here that, in our approach, the weights of the microscopic interactions (like repulsive, electrostatic, polarization, and dispersion ones) occurring in ion/water small clusters differ noticeably for the latter two ions (see Figure 8 and the Supporting Information). This supports our model parameter assignment strategy, as it allows us to account in a balanced manner for the latter microscopic interactions, which favor/ disfavor ion solvation. Similar to theoretical studies on the solvation of monatomic halide ions,41,53−56 we found an overall strong sensitivity of the ion/water potential of mean force at air/water interfaces to the model parameters, in particular, of methyl/water interactions. Nevertheless, all the other ion solvation properties discussed in the present study appear to be less sensitive to the model parameters, confirming the reliability of all our conclusions. In summary, the present simulation approach appears to be well suited to study organic ions that involve ion−solvent HBs and hydrophobic alkyl−water interactions. It applies sophisticated modeling approaches that were developed for monatomic ions to organic ions that have been less investigated theoretically so far. First, our approach investigates directly the structural and energetic properties of ions solvated in water nanodroplets. This is useful because, despite their importance in aerosol chemistry and thus in atmospheric pollution,57−59 there is a lack of experimental data on ions solvated in nanodroplets. Beyond the droplets, our objective is to develop methods that can extrapolate to bulk ion solvation in natural and biological environments, to obtain solvation thermochemistry that cannot be obtained experimentally. We demonstrated that our approach can yield the solvation energies of organic and biological ions, and identify and quantify the physical terms that contribute to the thermochemistry of solvation.

■

■

(1) Lehninger, A. Biochemistry; Worth Publishers: New York, 1970. (2) Taft, R. W. Protonic Acidities and Basicities in the Gas Phase and in Solution: Substituent and Solvent Effects. Prog. Phys. Org. Chem. 1983, 14, 247−350. (3) Meot-Ner (Mautner), M. The Ionic Hydrogen Bond. Chem. Rev. 2005, 105, 213−284. (4) Arnett, E. M.; Jones, F. M.; Taagepera, M.; Henderson, W. G.; Beauchamp, J. L.; Holtz, D.; Taft, R. W. Complete Thermodynamic Analysis of the “Anomalous Order” of Amine Basicities in Solution. J. Am. Chem. Soc. 1972, 94, 4724−4726. (5) Aue, D. H.; Webb, H. M.; Bowers, M. T. Quantitative Relative Gas-phase basicities of alkylamines. Correlation with solution basicity. J. Am. Chem. Soc. 1972, 94, 4726−4728. (6) Henderson, W. G.; Taagepera, M.; Holtz, D.; McIver, R. T.; Beauchamp, J. L.; Taft, R. W. Methyl Substituent Effects in Protonated Aliphatic Amines and Their Radical Cations. J. Am. Chem. Soc. 1972, 94, 4728−4729. (7) Aue, D. H.; Webb, H. M.; Bowers, M. T. A Thermodynamic Analysis of Solvation Effects on the Basicities of Alkylamines. An Electrostatic Analysis of Substituent Effects. J. Am. Chem. Soc. 1976, 98, 318−329. (8) Aue, D. H.; Webb, H. M.; Bowers, M. T. Quantitative Proton Affinities, Ionization Potentials, and Hydrogen Affinities of Alkylamines. J. Am. Chem. Soc. 1976, 98, 311−317. (9) Aue, D. H.; Webb, H. M.; Bowers, M. T.; Liotta, C. L.; Alexander, C. J.; Hopkins, H. P. A Quantitative Comparison of Gas- and SolutionPhase Basicities of Substituted Pyridines. J. Am. Chem. Soc. 1976, 98, 854−856. (10) Munson, M. S. B. Proton Affinities and the Methyl Inductive Effect1. J. Am. Chem. Soc. 1965, 87, 2332−2336. (11) Kebarle, P. Ion Thermochemistry and Solvation from Gas-Phase Ion Equilibria. Annu. Rev. Phys. Chem. 1977, 28, 445−476. (12) Meot-Ner, M.; Sieck, L. W. Proton Affinity Ladders from Variable-Temperature Equilibrium Measurements. 1. A Reevaluation of the Upper Proton Affinity Range. J. Am. Chem. Soc. 1991, 113, 4448−4460. (13) Szulejko, J. E.; McMahon, T. B. Progress toward an Absolute Gas-Phase Proton Affinity Scale. J. Am. Chem. Soc. 1993, 115, 7839− 7848. (14) Hunter, E. P.; Lias, S. G. Evaluated Gas Phase Basicities and Proton Affinities of Molecules: An Update. J. Phys. Chem. Ref. Data 1998, 27, 413−656. (15) Meot-Ner (Mautner), M. Update 1 of: Strong Ionic Hydrogen Bonds. Chem. Rev. 2012, 112, PR22−PR103. (16) Brauman, J. I.; Riveros, J. M.; Blair, L. K. Gas-Phase Basicities of Amines. J. Am. Chem. Soc. 1971, 93, 3914−3916. (17) Bartmess, J. E.; Scott, J. A.; McIver, R. T. Substituent and Solvation Effects on Gas-Phase Acidities. J. Am. Chem. Soc. 1979, 101, 6056−6063. (18) Meot-Ner, M. Heats of Hydration of Organic ions: Predictive Relations and Analysis of Solvation Factors Based on Ion Clustering. J. Phys. Chem. 1987, 91, 417−426. (19) Brugé, F.; Bernasconi, M.; Parrinello, M. Ab Initio Simulation of Rotational Dynamics of Solvated Ammonium Ion in Water. J. Am. Chem. Soc. 1999, 121, 10883−10888. (20) Hrobárik, T.; Vrbka, L.; Jungwirth, P. Selected Biologically Relevant Ions at the Air/Water Interface: A Comparative Molecular Dynamics Study. Biophys. Chem. 2006, 124, 238−242. (21) Gopalakrishnan, S.; Jungwirth, P.; Tobias, D. J.; Allen, H. C. Air/Liquid Interfaces of Aqueous Solutions Containing Ammonium and Sulfate: Spectroscopic and Molecular Dynamics Studies. J. Phys. Chem. B 2005, 109, 8861−8872. (22) Douady, J.; Calvo, F.; Spiegelman, F. Structure, Stability, and Infrared Spectroscopy of (H2O)nNH4+ clusters: A Theoretical Study at Zero and Finite Temperature. J. Chem. Phys. 2008, 129, 154501.

ASSOCIATED CONTENT

S Supporting Information *

Full description of the quantum computations, of the model parameter assignment strategy, of the molecular dynamics details, and of data mentioned in the text. This material is available free of charge via the Internet at http://pubs.acs.org.

■

REFERENCES

AUTHOR INFORMATION

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Othman Bouizi and Emmanuel Oseret (Exascale Computing Research Laboratory, a joint INTEL/CEA/UVSQ/ GENCI laboratory) for their help in optimizing the code POLARIS(MD). We acknowledge access to the supercomputing systems of the CCRT (Centre de Calcul et de Recherche Technologique) of the French Nuclear Agency (CEA) and the 6232



Article

(23) Babiaczyk, W. I.; Bonella, S.; Guidoni, L.; Ciccotti, G. Hydration Structure of the Quaternary Ammonium Cations. J. Phys. Chem. B 2010, 114, 15018−15028. (24) Chang, T. M.; Cooper, R. J.; Williams, E. R. Locating Protonated Amines in Clathrates. J. Am. Chem. Soc. 2013, 135, 14821−14830. (25) Thomson, J. J. Conduction of Electricity through Gases; C. J. Clay and Sons: London, 1903; pp 149−150. (26) Holland, P. M.; Castleman, A. W. Thomson Equation Revisited in Light of Ion-Clustering Experiments. J. Phys. Chem. 1982, 86, 4181−4188. (27) Peslherbe, G. H.; Ladanyi, B. M.; Hynes, J. T. Cluster Ion Thermodynamic Properties: The Liquid Drop Model Revisited. J. Phys. Chem. A 1999, 103, 2561−2571. (28) Yu, F. Modified Kelvin-Thomson Equation Considering IonDipole Interaction: Comparison with Observed Ion-Clustering Enthalpies and Entropies. J. Chem. Phys. 2005, 122, 084503. (29) Donald, W. A.; Williams, E. R. Evaluation of Different Implementations of the Thomson Liquid Drop Model: Comparison to Monovalent and Divalent Cluster Ion Experimental Data. J. Phys. Chem. A 2008, 112, 3515−3522. (30) Bragg, A. E.; Verlet, J. R. R.; Kammrath, A.; Cheshnovsky, O.; Neumark, D. M. Hydrated Electron Dynamics: From Clusters to Bulk. Science 2004, 306, 669−675. (31) Donald, W. A.; Leib, R. D.; Demireva, M.; O’Brien, J. T.; Prell, J. S.; Williams, E. R. Directly Relating Reduction Energies of Gaseous Eu(H2O)n3+, n = 55−140, to Aqueous Solution: The Absolute SHE Potential and Real Proton Solvation Energy. J. Am. Chem. Soc. 2009, 131, 13328−13337. (32) Donald, W. A.; Leib, R. D.; Demireva, M.; Williams, E. R. Ions in Size-Selected Aqueous Nanodrops: Sequential Water Molecule Binding Energies and Effects of Water on Ion Fluorescence. J. Am. Chem. Soc. 2011, 133, 18940−18949. (33) O’Brien, J. T.; Williams, E. R. Effects of Ions on HydrogenBonding Water Networks in Large Aqueous Nanodrops. J. Am. Chem. Soc. 2012, 134, 10228−10236. (34) Fujii, A.; Mizuse, K. Infrared Spectroscopic Studies on Hydrogen-Bonded Water Networks in Gas Phase Clusters. Int. Rev. Phys. Chem. 2013, 32, 266−307. (35) Horvárth, L.; Beu, T.; Manghi, M.; Palmeri, J. The Vapor-Liquid Interface Potential of (Multi)Polar Fluids and its Influence on Ion Solvation. J. Chem. Phys. 2013, 138, 154702. (36) Réal, F.; Vallet, V.; Flament, J.-P.; Masella, M. Revisiting a Many-Body Model for Water Based on a Single Polarizable Site. From Gas Phase Clusters to Liquid and Air/Liquid Water Systems. J. Chem. Phys. 2013, 139, 114502. (37) Thole, B. Molecular Polarizabilities Calculated with a Modified Dipole Interaction. Chem. Phys. 1981, 59, 341−350. (38) Lee, H. M.; Kim, J.; Lee, S.; Mhin, B. J.; Kim, K. S. Aqua/ Potassium(I) Complexes: Ab Initio Study. J. Chem. Phys. 1999, 111, 3995−4004. (39) Lee, H. M.; Tarakeshwar, P.; Park, J.; Kolaski, M. R.; Yoon, Y. J.; Yi, H.-B.; Kim, W. Y.; Kim, K. S. Insights into the Structures, Energetics, and Vibrations of Monovalent Cation/(Water)1−6 Clusters. J. Phys. Chem. A 2004, 108, 2949−2958. (40) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: London, 2002. (41) Caleman, C.; Hub, J.; van Maaren, P.; van der Spoel, D. Atomistic Simulation of Ion Solvation in Water Explains Surface Preference of Halides. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 6838− 6842. (42) Warren, G. L.; Patel, S. Hydration Free Energies of Monovalent Ions in Transferable In- termolecular Potential Four Point Fluctuating Charge Water: An Assessment of Simulation Methodology and Force Field Performance and Transferability. J. Chem. Phys. 2007, 127, 064509. (43) Yu, H.; Whitfield, T. W.; Harder, E.; Lamoureux, G.; Vorobyov, I.; Anisimov, V. M.; MacKerell, A. D.; Roux, B. Simulating Monovalent

and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field. J. Chem. Theory Comput. 2010, 6, 774−786. (44) Dang, L. X.; Chang, T.-M. Molecular Mechanism of Ion Binding to the Liquid/Vapor Interface of Water. J. Phys. Chem. B 2002, 106, 235−238. (45) Tian, C.; Byrnes, S. J.; Han, H.-L.; Shen, Y. R. Surface Propensities of Atmospherically Relevant Ions in Salt Solutions Revealed by Phase-Sensitive Sum Frequency Vibrational Spectroscopy. J. Phys. Chem. Lett. 2011, 2, 1946−1949. (46) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Aqueous Solvation Free Energies of Ions and Ion: Water Clusters Based on an Accurate Value for the Absolute Aqueous Solvation Free Energy of the Proton. J. Phys. Chem. B 2006, 110, 16066−16081. (47) Fan, Y.; Chen, X.; Yang, L.; Cremer, P. S.; Gao, Y. Q. On the Structure of Water at the Aqueous/Air Interface. J. Phys. Chem. B 2009, 113, 11672−11679. (48) in ’t Veld, P. J.; Ismail, A. E.; Grest, G. S. Application of Ewald Summations to Long-Range Dispersion Forces. J. Chem. Phys. 2007, 127, 144711. (49) Tuttle, T. R.; Malaxos, S.; Coe, J. V. A New Cluster Pair Method of Determining Absolute Single Ion Solvation Energies Demonstrated in Water and Applied to Ammonia. J. Phys. Chem. A 2002, 106, 925− 932. (50) Meot-Ner (Mautner), M.; Lias, S. G. In Binding Energies Between Ions and Molecules, and the Thermochemistry of Cluster Ions; Linstrom, P. J., Mallard, W. G., Eds.; National Institute of Standards and Technology; http://webbook.nist.gov (retrieved Dec 11, 2013). (51) Shi, Z.; Ford, J. V.; Wei, S.; Castleman, A. W. Water clusters: Contributions of Binding Energy and Entropy to Stability. J. Chem. Phys. 1993, 99, 8009−8015. (52) Kelterbaum, R.; Kochanski, E. Behavior and Evolution of the First 28 Protonated Hydrates from Monte Carlo Studies. J. Phys. Chem. 1995, 99, 12493−12500. (53) Dang, L. X. Computational Study of Ion Binding to the Liquid Interface of Water. J. Phys. Chem. B 2002, 106, 10388−10394. (54) Wick, C. D. Electrostatic Dampening Dampens the Anion Propensity for the Air-Water Interface. J. Chem. Phys. 2009, 131, 084715. (55) Wick, C. D.; Cummings, O. T. Understanding the Factors that Contribute to Ion Interfacial Behavior. Chem. Phys. Lett. 2011, 513, 161−186. (56) Netz, R. R.; Horinek, D. Progress in Modeling of Ion Effects at the Vapor/Water Interface. Annu. Rev. Phys. Chem. 2012, 63, 401−418. (57) Finlayson-Pitts, B. J. Reactions at Surfaces in the Atmosphere: Integration of Experiments and Theory as Necessary (but not Necessarily Sufficient) for Predicting the Physical Chemistry of Aerosols. Phys. Chem. Chem. Phys. 2009, 11, 7760−7779. (58) Rubasinghege, G.; Elzey, S.; Baltrusaitis, J.; Jayaweera, P. M.; Grassian, V. H. Reactions on Atmospheric Dust Particles: Surface Photochemistry and Size-Dependent Nanoscale Redox Chemistry. J. Phys. Chem. Lett. 2010, 1, 1729−1737. (59) Abbatt, J. P. D.; Lee, A. K. Y.; Thornton, J. A. Quantifying Trace Gas Uptake to Tropospheric Aerosol: Recent Advances and Remaining Challenges. Chem. Soc. Rev. 2012, 41, 6555−6581. (60) Klots, C. E. Solubility of Protons in Water. J. Phys. Chem. 1981, 85, 3585−3588.

6233


Simulated Solvation of Organic Ions: Protonated Methylamines in

Recommend Documents