Molecular Dynamics Simulations to Examine Structure, Energetics

Oct 2, 2012 - Molecular Dynamics Simulations to Examine Structure, Energetics, and Evaporation/Condensation Dynamics in Small Charged Clusters of ...
0 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/JPCA

Molecular Dynamics Simulations to Examine Structure, Energetics, and Evaporation/Condensation Dynamics in Small Charged Clusters of Water or Methanol Containing a Single Monatomic Ion Christopher D. Daub and Natalie M. Cann* Department of Chemistry, Queen’s University, 90 Bader Lane, Kingston, Ontario, Canada K7L3N6 ABSTRACT: We study small clusters of water or methanol containing a single Ca2+, Na+, or Cl− ion with classical molecular dynamics simulations, using models that incorporate polarizability via the Drude oscillator framework. Evaporation and condensation of solvent from these clusters is examined in two systems, (1) for isolated clusters initially prepared at different temperatures and (2) those with a surrounding inert (Ar) gas of varying temperature. We examine these clusters over a range of sizes, from almost bare ions up to 40 solvent molecules. We report data on the evaporation and condensation of solvent from the clusters and argue that the observed temperature dependence of evaporation in the smallest clusters demonstrates that the presence of heated gas alone cannot, in most cases, solely account for bare ion production in electrospray ionization (ESI), neglecting the key contribution of the electric field. We also present our findings on the structure and energetics of the clusters as a function of size. Our data agree well with the abundant literature on hydrated ion clusters and offer some novel insight into the structure of methanol and ion clusters, especially those with a Cl− anion, where we observe the presence of chain-like structures of methanol molecules. Finally, we provide some data on the reparameterizations necessary to simulate ions in methanol using the separately developed Drude oscillator models for methanol and for ions in water.

1. INTRODUCTION

We have now continued this work by leaving for the moment the effect of the electric field and focusing instead on the evaporation and condensation of solvent molecules from small clusters (up to 40 solvent molecules + 1 ion). We consider both water and methanol to make the results more relevant to ESI, where the most common solvent is a water/methanol mixture. We simulate both a cluster in vacuum as well as the additional effect of a surrounding gas. By heating either the gas or the initial cluster above room temperature, we can simulate the effects of heating either the curtain gas17 or the ion capillary tube,18 as is commonly done in modern ESI-MS instrumentation.19,20 There are two main goals for this study. First, we wish to quantify the ability of increased temperature, either in the charged cluster or in the surrounding gas, to remove solvent from the cluster in the absence of electric field. Our new findings reinforce the conclusions of our initial study,12 demonstrating that gas alone, even quite hot gas, cannot be expected to remove all of the solvent in most of the systems that we studied, and therefore, the electric field must play a role as well. Second, we use long simulations of clusters of varying sizes that are allowed to evaporate and condense solvent as a convenient means to examine the energetics and structure of a broad range of cluster sizes, from nearly bare ions up to 40 solvent molecules. We are particularly interested in studying the

1−3

Electrospray ionization (ESI) is a popular method for transferring large ions into the gas phase for analysis by mass spectrometry (MS). The process entails the efficient loss of solvent, by both droplet fission and solvent evaporation, until a bare or nearly bare ion can be obtained. Several interesting ESIfocused molecular simulation studies have considered the intermediate stages of the process, with hundreds4,5 or thousands6−11 of remaining solvent molecules. Our interest is the removal of the final few solvent molecules in the last stage of ESI. In a recent publication,12 we reported classical molecular dynamics (MD) simulations of small clusters of hydrated ions accelerated by an electric field through a gas, the first time gas molecules have been explicitly considered in an atomistic simulation of ESI. Our findings stressed the importance of both the electric field and gas on the removal of the final solvent molecules from the ion, suggesting that collision-induced solvent loss from the accelerated cluster interacting with surrounding gas must be considered in new models for ESI mechanisms that go beyond the debate8,11,13,14 between charged-residue models1 and ion evaporation models,2,15 at least where simple monatomic ions are concerned. There has already been some experimental support for our ideas in recent ESI-MS experiments on brine solutions, where it was found that the magnitude of the electric potential, that is, so-called soft or hard ionization conditions, could make a difference in allowing the complete removal of the solvent from the ion.16 © 2012 American Chemical Society

Received: August 17, 2012 Revised: September 25, 2012 Published: October 2, 2012 10488

dx.doi.org/10.1021/jp308217q | J. Phys. Chem. A 2012, 116, 10488−10495

The Journal of Physical Chemistry A

Article

freedom ndof in addition to three degrees of freedom for the Drude particle. The Drude particle is kept cooled with a Nose− Hoover chain thermostat43 to 1 K and thus can be neglected in conversions between the kinetic energy and temperature. Minor technical differences38,42 between the charge on spring and Drude oscillator methodologies developed by the two groups did not affect our results in any noticeable way. The same ion models used with the NDP water model were used as a base for the methanol simulations but required some reparameterization. These modifications have been detailed in section 3.1. Before proceeding further, we take a moment to define some of the quantities that we compute. The total cluster kinetic energy will be denoted by K. It is computed from the instantaneous velocities of all of the molecules in the cluster, including the ion, with membership in the cluster determined by the cluster analysis described below. We also subtract the cluster center-of-mass velocity so that conversions between K and the cluster temperature Tc are accurate. This is especially important in the smallest clusters where recoil from the condensing or evaporating solvent molecule can be significant. The total cluster potential energy will be denoted by U. It is computed from the total configurational energy of the system, including pairwise electrostatic and van der Waals terms plus the spring energy in the Drude oscillators. K is determined only for the cluster, and U is determined for the whole system. Because the residual potential energy of the evaporated molecules is negligible, these energies can be usefully compared. We detect strongly interacting solvent-only pairs, triples, and so forth during our cluster analysis, and we do not include these configurations when we compute U. The initial cluster configuration was obtained by first equilibrating a bulk configuration of solvent and dilute ions using Monte Carlo moves in the absence of electrostatic forces. Then, the initial cluster was formed by selecting the closest ion to the center of the simulation box and the nearest solvent molecules up to the desired number. The simulation box was then greatly expanded by a factor of 10 to avoid any effects from interactions of the cluster with its periodic images but still allow evaporated solvent molecules to interact with the cluster at a later time. The final cubic simulation box had sides of length 234.8 Å in the water clusters and 238.8 Å in the methanol case. In the first 50 ps, the temperature was gradually increased from 150 K up to the desired initial temperature. A radial confining potential was applied to hold the cluster together while this temperature equilibration took place. All of our cluster simulations computed all of the minimum-image pairwise intermolecular interactions, both electrostatic and Lennard-Jones. Only our test simulations for parametrization purposes in section 3.1 used Ewald summation. We have mainly done two types of simulations. First, we have considered clusters initially in vacuum, with different numbers of solvent atoms and starting temperatures ranging from below room temperature (275 K) up to very hot temperatures necessary in some cases to sample the smallest clusters (∼2000 K). After the initial thermalization, the thermostat was turned off (NVE ensemble; but note again that the thermostat on the Drude particle remained at 1 K). Depending on the initial temperature, some number of solvent molecules evaporated until an equilibrium with the remaining cluster near room temperature was reached. After this point, we tracked the evaporation and condensation events and computed the average kinetic energy of the cluster, as well as the total

changes that accompany these evaporation and condensation events. We have introduced polarizable models for both the solvent and ions. This should allow us to accurately represent the reciprocal polarization between the ion and the surrounding solvent. There is a general consensus that polarizable models are more accurate at describing anions at aqueous interfaces,21−23 although there is still some debate as whether the polarizability per se is the main factor determining the interfacial ion behavior.24 The question of whether ions, especially large anions, have a preference for the liquid−vapor interface rather than the bulk has been a topic of intense study. There is little we can add to this discussion regarding anions in water, where it is fairly well established that anions other than fluoride seek the surface of the air−liquid interface, both in bulk salt solutions21−28 and in single-ion clusters,29 with a transition in the clusters from surface to interior solvation with increasing temperature30 or increasing cluster size.31 There is much less literature on solutions with other polar solvents. Simulations of solvents without hydrophobic groups (e.g., formamide, ammonia) showed a surface preference for anions similar to that seen at the aqueous liquid-vapor interface;32 however, it appears in methanol that the methyl group’s own preference for the liquid−vapor interface is strong enough to keep the anion in the solution interior.33,34 The specific question regarding anions in methanol that we can address is whether the same behavior observed at the liquid−vapor interface of a salt solution carries over to a charged cluster. The only other papers to address this reported surface solvation in small (three and four methanol molecules) Cl− clusters.35,36 Our own answer to the question is nuanced; instead of strictly arguing for surface or interior solvation, we stress the importance of extreme nonsphericity in the methanol−Cl− cluster and the presence of long chains of methanol molecules, long known to be a dominant motif in liquid methanol.37 The rest of the paper is organized as follows. In section 2, we describe the potential models and simulation methods that we used. Section 3 details the results of our study. In section 3.1, we present the adjustments required to reparameterize the polarizable ion models for use with methanol solvent. Section 3.2 concerns our analysis of evaporation and condensation events and focuses on the cluster temperature changes accompanying these events. In section 3.3, we turn to the structure and energetics of the clusters. Finally, in section 4, we summarize, draw some general conclusions, and make some suggestions for future work.

2. SIMULATION METHODS We perform classical MD simulations of clusters of H2O or CH3OH molecules and one Na+, Ca2+, or Cl− ion. The initial clusters have either 10, 20, 30, or 40 solvent molecules. In our prevous work,12 we did not incorporate polarizability in the potential models. Now that we wish to include anions in our simulations, it is important to include polarization effects explicitly. For the water and ions, we have chosen to use the Drude oscillator models developed by Lamoureux et al.38−40 and have updated our in-house code to incorporate this methodology. The specific models chosen were the negative Drude particle (NDP) water model40 and associated ion force field parameters developed by Yu et al.41 For methanol, we have adopted the so-called charge on spring methanol (COS/ M) model, also developed by Yu with different co-workers.42 Both of these models have a total of six molecular degrees of 10489

dx.doi.org/10.1021/jp308217q | J. Phys. Chem. A 2012, 116, 10488−10495

The Journal of Physical Chemistry A

Article

height and the solvent−ion coordination number Nc, with Nc defined here and throughout the paper as the average number of solvent molecules separated from the ion by less than the radius of the first minimum in the ion−solvent oxygen radial distribution function giO(r) measured in the bulk solution. We have done the required modifications by altering the mixing rules for the Lennard-Jones size parameter σ and the products of charges. The advantage of this strategy is that it will allow us in the future to mix together water and methanol in the same simulation with minimal difficulty. Surprisingly, the Ca2+ model could be used without alteration. The modifications used for all ions are listed in Table 1, and the resulting bulk ion−oxygen radial distribution functions are graphed in Figure 1.

potential energy of the system, before and after these events, for a total of 8 ns. This allowed us to build up statistics for the changes in the kinetic and potential energy in the cluster during these events as a function of the cluster size. To allow for the cluster size and temperature to settle into an equilibrium we did not initiate data collection until 1.5 ns had elapsed after the initial cluster was prepared. It is worth explaining why we simulate the clusters in a large simulation box with periodic boundaries rather than an isolated cluster in vacuum. One reason is quite simply that one of the main goals of this work is to study evaporation and recondensation, and this would be impossible to do with an isolated cluster. In addition, although at first it may seem that the isolated cluster is a better representation of the conditions in ESI experiments, we contend that during droplet breakup, there must be some concentration of solvent vapor present as well. Explicit addition of solvent vapor to ESI bath gas has been studied, albeit at concentrations on the order of 1 solvent molecule per 1000 gas molecules, somewhat lower than that in our simulations.14 For the second type of simulations, clusters with either 10 or 20 solvent molecules were prepared as described above, but we also included a surrounding argon gas at P = 1 atm. In these cases, the initial cluster started at Tc = 300 K. Some commercial ESI instruments introduce a gas with a temperature as high as 1050 K;20 we set the gas at one of four different initial temperatures ranging from Tg = 300 to 900 K. This allowed us to explore the effect of thermalization by a curtain gas. The 32 ns simulations were sufficient to reach the eventual equilibrium cluster size in all but a couple of our studied systems. The Ar Lennard-Jones parameters are from the classic simulation of Rahman.44 The same minimum-image convention was used for the gas and gas−cluster Lennard-Jones interactions as that used for the solvent and ion interactions. A cluster analysis45 was performed every 0.1 ps to track the size of the cluster around the ion. A molecule was considered to be part of the cluster as long as one of its atoms remained within 1.5σij = 0.75(σi + σj) of any atom in any other molecule in the cluster, for example, 4.34 Å for water and Na+. With this criterion most transient, subpicosecond hydrogen bond breaking and remaking events due to librations were ignored, but in some cases, the cluster size still fluctuated on 1−2 ps time scales due to hydrogen bonds. To eliminate such fluctuations, we smoothed the cluster size curves by averaging over 100 total determinations (10 ps) before and after the given time, and further, we accumulated the results of at least three separate simulations at each temperature to obtain reliable data. We also used the cluster analysis to track the presence of solvent-only clusters. As we will discuss further below, such solvent-only clusters only occurred with significant frequency in our simulations of the Cl− anion with methanol.

Table 1. Parameterization of Ion Models with Methanola ion

aσ

aq2

Nc

giOM(rpeak)

rpeak/Å

Ca Ca2+, CPMD48 Na+ Na+, CPMD47 Cl− Cl−, CPMD46

1.0

1.0

1.04

0.95

0.98

0.65

6.0 6.0 5.0 4.8−5.0 3.7 3.6

25 22 15.5 12 4.2 5.0

2.4 2.4 2.45 2.4−2.5 3.2 3.2

2+

a

aσ is a multiplicative factor applied to the mixed Lennard-Jones size parameters σixM = aσ(σi + σxM)/2. aq2 is a multiplicative factor applied to the products of the charges in the electrostatic interactions between the ion and methanol (Ecoul = aq2(qiqxM/r)). xM refers to all atomic centers in the methanol molecule, and both the atomic charges and the Drude oscillator charges are considered in scaling Ecoul. Nc is the coordination number of the ion, and giOM(rpeak) and rpeak are the height and position, respectively, of the first peak in giOM(r).

Figure 1. Bulk giOM (r) for simulations of 256 methanol molecules with a single ion, with potential parameters adjusted according to Table 1.

3.2. Evaporation and Condensation in Hot Clusters and Hot Gases. We start by considering a cluster initially in vacuum. In Figure 2, we show an example of the running cluster size nclust, cluster kinetic energy per cluster molecule K/nclust, and system potential energy per cluster molecule U/nclust over the course of a single 8 ns trajectory. Note that we will consistently define nclust as the total cluster size, including the ion, throughout. As the initial temperature of the cluster is high (Tc = 550 K), we observe about half of the initial cluster’s 20 solvent molecules very quickly evaporating. These vapor molecules are now available to be recondensed by the cluster. It can be seen that each evaporation and condensation event is accompanied by a resulting change in the K and U. By tracking the evaporation and condensation events for many different

3. RESULTS AND DISCUSSION 3.1. Parameterization of Ion Models with Methanol. The ion models introduced for use with the NDP water model41 required some reparameterization for use with the COS/M methanol model. For this purpose only, we use test simulations of a single solvated ion in bulk methanol (256 solvent, 1 ion, periodic boundary conditions, NVT ensemble, and Ewald summation). This system size is similar to that used originally by Yu et al.41 We compare with the available Car− Parrinello MD (CPMD) structural data for single ions in bulk methanol.46−48 We try to match the first peak position and 10490

dx.doi.org/10.1021/jp308217q | J. Phys. Chem. A 2012, 116, 10488−10495

The Journal of Physical Chemistry A

Article

Figure 2. Example of one trajectory for a simulation without curtain gas with an initial cluster of 1 Na+ + 20 H2O at an initial temperature of Tc = 550 K, showing the running cluster size nclust, the cluster kinetic energy K/nclust, and the negative of the potential energy U/nclust. Note that the potential energy is divided by 4 to fit it on the same graph. The blue dots indicate the times where an evaporation or condensation event occurred, along with the cluster kinetic energy averaged over the time since the previous event .

initial cluster sizes and temperatures, data can be averaged over many events, and thus, the whole range of cluster sizes from less than one solvation shell up to the largest initial clusters (40 solvent) can be scanned. The change in cluster temperature, ΔTc, accompanying each event can be obtained from the change in cluster ⟨K⟩.

Figure 3. Plots of the average ΔTc,cond, the temperature increase in the cluster accompanying a single-molecule condensation event, and the average −ΔTc,evap, the negative of the temperature loss in the cluster accompanying an evaporation event as a function of the cluster size after this event. Error bars are one standard deviation in the simulated data.

⎛ ⟨KN + 1⟩ ⟨KN ⟩ ⎞ 2 ΔTc,cond = ⎜ − ⎟ Nndof ⎠ kB ⎝ (N + 1)ndof ⎛ ⟨K ⟩ ⟨KN − 1⟩ ⎞ 2 −ΔTc,evap = ⎜ N − ⎟ (N − 1)ndof ⎠ kB ⎝ Nndof

larger. It is also remarkable how similar the curves are for Na+ in the two different solvents, while for Cl− and Ca2+, ΔT becomes quite a bit smaller in methanol at the smaller cluster sizes. In Table 2, we show ΔT for all ions and solvents at two selected cluster size transitions.

(1)

where N = nclust is the cluster size before evaporation or condensation, ndof = 6 are the molecular degrees of freedom, and kB is the Boltzmann constant. We show these results in Figure 3. The changes in ⟨K⟩ and ⟨U⟩ themselves are discussed later in section 3.3 and plotted in Figure 7, but for now, we wish to focus on the temperature changes. We include a negative sign for ΔTc,evap to allow direct comparison between evaporation and condensation. For a condensation (evaporation) event, ⟨KN⟩ is averaged over the time elapsed since the previous event, and ⟨KN+1⟩ (⟨KN−1⟩) is averaged over the time period up until the next event. As shown in Figure 2, to ensure that the system has equilibrated properly, we only start gathering data after 1.5 ns. In averaging the data for Figure 3, we also do not include transient events, which we define as occurring when the cluster size does not remain constant for at least 50 ps, nor do we include cases involving multiple solvent molecules participating in the same event. These multiple molecule events are very rare for the cations but more prevalent in simulations with Cl− and especially so in simulations with Cl− and methanol. The temperature changes during evaporation and condensation are, to within statistics, equal. Further, all three ions show a pronounced increase in ΔT at the smaller cluster sizes. Some of this increase is simply due to there being fewer molecules remaining in the cluster to absorb the loss (gain) in kinetic energy from the evaporating (condensing) molecule. However, especially in the case of the cations, much of the increase in ΔT comes about from the increased binding energy (cf. Figure 7), which must be overcome to remove the final few solvent molecules directly bound to the ion. Ca2+ in particular binds the solvent very tightly, and hence, the temperature changes are

Table 2. Temperature Change ΔT Accompanying Transitions between Two Sizes of Cluster Due to Evaporation or Condensationa

a

ion

solvent

ΔT(4 ↔ 5) (K)

ΔT(10 ↔ 11) (K)

Ca2+ Ca2+ Na+ Na+ Cl− Cl−

H2O MeOH H2O MeOH H2O MeOH

618 429 189 192 134 80

87 77 52 50 36 39

The values for evaporation and condensation have been averaged.

For our simulations of clusters in gas, all of the clusters are initially equilibrated at a temperature of Tc = 300 K, while the gas is separately thermalized at a different temperature Tg. This lets us observe the ability of a hot gas, as used in some ESI-MS instruments, to desolvate an ion. After the initial 50 ps equilibration period, the thermostats are removed (NVE ensemble except for the Drude thermostat). When the gas and cluster start at the same temperature, there is some system cooling due to a small amount of cluster evaporation. Further solvent loss does not occur due to the large temperature change that must accompany evaporation (cf. Figure 3). When the gas temperature is elevated, the cluster and gas exchange some kinetic energy, with the gas cooling somewhat and the cluster heating. This heating causes some additional solvent evaporation. Eventually, an equilibrium will be reached. In Figure 4, we show the average cluster size as a function of the simulation 10491

dx.doi.org/10.1021/jp308217q | J. Phys. Chem. A 2012, 116, 10488−10495

The Journal of Physical Chemistry A

Article

temperature of ΔTc(nclust) = 200 K can be expected to cause evaporation, subtract 1 for the evaporating molecule, and compare it with the final nclust for the cluster starting out with 10 solvent molecules in a gas initially at Tg = 500 K, that is, 300 + 200 K. The comparison is generally good, with the slightly larger cluster size sometimes seen in the gas simulations largely attributable to the lowered gas temperature due to equilibration between the gas and cluster, in addition to long equilibration times especially seen in the cluster with Ca2+, where even 32 ns was not sufficient to see the final cluster size. Deviations from equality would be larger for the clusters starting with 20 solvent molecules due to increased gas cooling and longer equilibration times. 3.3. Energetics and Structure of Clusters. In this section, we report some data for the energetics and structure of clusters initially prepared in vacuum. Similar data have been reported many times for hydrated ions;29,30 however, our results differ in two important ways, which make them a useful addition to the existing literature. First, our method for preparing clusters with a broad range of initial sizes and temperatures allows us to scan the whole range of clusters from the bare ion up to the largest size that we simulated (40 solvent molecules), rather than only clusters of a specific size. Second, we are not limited to studying clusters in minimum-energy configurations or at reduced temperatures; instead, we examine clusters at whatever temperature they happen to have in our simulation. This comes closer, we believe, to the real conditions experienced by clusters in ESI experiments. Our final equilibrated clusters tend to have temperatures between ∼270 and 300 K, but this can be much higher in some of our runs to examine the smaller clusters. In some cases, the starting temperature was as high as 2000 K, and the final cluster with less than a full solvation shell can be stable at temperatures much higher than 300 K. These conditions compare well with ESI experiments using heated gas, where gas temperatures as high as 1050 K are used to aid desolvation, but evaporative cooling results in a droplet with a steady-state temperature of around 317 K.20,49 We start by presenting some snapshots of the systems in Figure 5. Many of the numerical results that we will report can

Figure 4. Plot of cluster size as a function of the time in simulations with argon at P = 1 atm and different initial gas temperatures Tg. Solid lines are for an initial cluster size of 1 ion + 20 H2O, and dashed lines are for 1 ion + 10 H2O. The data are averaged over six trajectories for each initial cluster size and gas temperature. The initial cluster temperature is 300 K in all cases.

time for all of the simulations with gas. Even after long simulations (32 ns), in most cases, the last few solvent molecules remain bound to the ion, despite the presence of very hot gas. This supports the findings from our first publication,12 indicating that only collisions with fieldaccelerated clusters can be expected to quickly remove the last few tightly bound solvent molecules. The exception is the simulations with Cl−, where moderately hot gas (Tg ≈ 500 K in methanol, 900 K in water) is sufficient to remove all of the solvent. It is interesting to compare the data from clusters in vacuum in Figure 3 with the data for clusters in hot gas in Figure 4. Consider a cluster of size nclust that has a temperature Tc = 300 K. To a first approximation, for a gas molecule to cause an evaporation event from this cluster, it must be at a minimum temperature of Tg = Tc + ΔTc(nclust), with the ΔTc(nclust) read off of Figure 3. This is of course a rough argument, with many assumptions underlying it. Nevertheless, the data do appear to agree with this equation fairly well. For example, in Table 3, we show the lowest cluster size where a change in cluster Table 3. Comparison of Final Cluster Size (nclust − 1), Where ΔTc(nclust) = 200 K Can Be Expected to Cause Evaporation, Taken from Figure 3, with the Final nclust for the Clusters That Initially Contain 10 Solvent Molecules in a 500 K Gas, Taken from Figure 4 final cluster size ion

solvent

nclust − 1 (Figure 3)

nclust (Figure 4)

Ca2+ Ca2+ Na+ Na+ Cl− Cl−

H2O MeOH H2O MeOH H2O MeOH

6 5 4 4 3 1

9 8 5 4 5 1

Figure 5. Snapshots of representative clusters. The solvent is shown only as a wire frame. Colors are as follows: O, red; H, white; CH3, cyan; Na+, orange; Cl−, green; Ca2+, silver. 10492

dx.doi.org/10.1021/jp308217q | J. Phys. Chem. A 2012, 116, 10488−10495

The Journal of Physical Chemistry A

Article

difference in ⟨U⟩ at the two temperatures is very small; in fact, for such small clusters, the instantaneous fluctuations in temperature are large enough to overwhelm any significant temperature effect. For all three ions, the potential energy is lower for water than it is for methanol. As one would expect for the cations, the total potential energy per cluster molecule drops to a minimum at slightly less than a full solvation shell. Clusters smaller than this have less energy per molecule due to the fact that there are too few solvent molecules to allow for sufficient solvent−solvent hydrogen bonding interactions. For Cl−, there is no minimum apparent in either solvent. This surprising result can be explained for hydrated Cl− by the preference for the anion to remain on the surface of the cluster. Even in the smaller Cl− clusters, the first solvation shell is not filled. The situation is even more interesting in the larger Cl− ion plus methanol cluster. Here, the tendency for methanol to form long chains of molecules as seen in the snapshot above greatly reduces the binding energy, so that the potential energy is less negative than it would be in bulk methanol. Clusters of Cl− in water have more or less the same U per molecule as bulk water as long as nclust ≳ 10. Clusters with Na+ in both solvents have approximately the same U per molecule for nclust ≳ 20. The total potential energy for these larger Na+ clusters with water is very close to the value for the bulk solvent. Surprisingly, the potential energy for the larger Na+ clusters with methanol remains ∼3 kJ/mol more negative than the bulk value; this could be due to the methanol clusters having a lower average temperature compared to the simulations of bulk methanol. For Ca2+, even the largest clusters have significantly more binding energy per molecule than the bulk solvent as the doubly charged cation exerts proportionally more influence. In Figure 7, we show our data for the changes in cluster kinetic energy and system potential energy, ΔK and ΔU, accompanying evaporation and condensation events. The data collected for these quantities is subject to the same restrictions

be explained by referring to these snapshots. We show one smaller cluster and one larger cluster for each of the six combinations of ion and solvent. Snapshots of the smaller Na+ clusters were chosen to show configurations with four solvent molecules in the first solvation shell of the ion and one more molecule in transition between the first and second shell. Snapshots of the aqueous Cl− clusters consistently show the surface solvation state preferred by anions other than fluoride, demonstrating the ability of the Drude oscillator models to capture this important property. In methanol, even the smallest Cl− clusters show a number of methanol molecules that choose to bond to each other rather than the ion. The larger Cl− clusters with methanol show a very interesting structure. The coordination number of the anion remains very low (between 2 and 3); however, unlike the case in water, those methanol molecules that bind directly to the anion show no strong asymmetry. In our estimation, it would be more accurate to describe the structure as two or three more-or-less independent chains of methanol molecules, forming a highly nonspherical cluster. In such a case, arguing about interior or surface solvation states is rather ambiguous. This structure also explains why multiple molecule evaporation events are so prevalent in the methanol−Cl− clusters; with such extended chains, it should not be a surprise that some fluctuation could break one of the hydrogen bonds anywhere along the chain. Finally, the clusters with Ca2+ unsurprisingly have a very rigid first solvation shell with six water or methanol molecules in the larger clusters. Only in a very small range of nclust, shown here for the smaller methanol cluster, is there any observation of a solvent molecule choosing to bind to other solvent rather than becoming one of the six coordinating solvent molecules bound to the Ca2+ ion. In Figure 6, we show the average potential energy ⟨U⟩ per cluster molecule. We note again that configurations including

Figure 6. Plots of ⟨U⟩, the energy of the solvated ion cluster, divided by the cluster size. The horizontal green lines show the bulk values of ⟨U⟩ per molecule for water (⟨U⟩ = −41.5 kJ/mol from ref 40) and methanol (⟨U⟩ = −35.3 kJ/mol from ref 42). The other line colors indicate the type of ion in the cluster: black, Ca2+; red, Na+; blue, Cl−. The line styles indicate the cluster temperature: solid, lower Tc; dashed, higher Tc.

solvent-only clusters are not included in these averages; thus, we can neglect the small potential energy of nonbonded solvent vapor and assume that all of the system’s potential energy belongs to the cluster. We show data at two different values of Tc for each of the six combinations of ion and solvent. The

Figure 7. Plot of the changes in the cluster K and system U and their sum during evaporation and condensation events. Energy changes during evaporation have a negative sign to allow direct comparison with condensation. Error bars are one standard deviation. 10493

dx.doi.org/10.1021/jp308217q | J. Phys. Chem. A 2012, 116, 10488−10495

The Journal of Physical Chemistry A

Article

Nc ≳ 4 from the largest clusters down to clusters with a single solvation shell. As shown in the snapshots, there is considerable lability in the entire Na+ solvation shell, with frequent transitions between the first solvation shell and the rest of the cluster rather than the other possibility of four tightly bound solvent molecules, with another one or two molecules exchanging more frequently. Finally, Ca2+ clusters have a very tightly bound solvation shell, with six solvent molecules, down to very small cluster sizes. It is only at a cluster size below ∼13 for methanol and ∼7 for water that any reduction below Nc = 6 can be observed.

as that shown in Figure 3, that is, after 1.5 ns of equilibration and not including transient or multiple molecule events. Because the total energy must be constant before and after the event in these NVE simulations and assuming the condensing or evaporating solvent molecule has no potential energy, the sum of ΔK + ΔU can be considered a reasonable estimate of the kinetic energy of the evaporating or condensing solvent molecule, Kvap, also plotted in Figure 7. Note once again that we are subtracting the center of mass kinetic energy, which explains why the sum of ΔK + ΔU drops below zero in some cases in the smallest clusters. Neglecting the smallest clusters, the quantity Kvap is remarkably constant for all cluster sizes, averaging around 7 kJ/mol, but with considerable deviation between individual events. This energy is reasonable because 3RT is the average kinetic energy for a molecule with 6 degrees of freedom at ∼280 K. It is worth noting that the change in potential energy accompanying evaporation and condensation is much less than the average binding energy of molecules in the solvent (see Figure 6), even in the largest clusters. This can be explained by the fact that these events must by necessity select molecules at the edge of the cluster, and hence, these molecules do not have as much binding energy as molecules closer to the center. Another interesting piece of data that we have already mentioned is the ion coordination number Nc. These results are in Figure 8. There is a general trend of the Nc for methanol

4. CONCLUSIONS We have completed a long series of simulations of clusters of water and methanol containing a single Na+, Cl−, or Ca2+ ion. We implemented polarization with the Drude oscillator methodology. We simulated clusters both in an initial vacuum, at different initial temperatures, as well as in a gas of elevated temperature, with the cluster initially at room temperature. We computed the cooling (heating) of the cluster, which accompanied evaporation (condensation) events. When we compared these temperature changes in an isolated cluster of a given size with the final cluster size in the heated gas simulations, we obtained reasonable agreement. Results for both systems show that even at rather high temperatures of the surrounding gas or the initially isolated cluster, the final few solvent molecules remain. We also found that evaporation of these last few solvent molecules must be accompanied by a very sharp drop in the temperature of the solvent, especially for the cations Na+ and Ca2+. These data reinforce the conclusions of our first study in this area,12 that observation of bare ions in ESI-MS experiments cannot (in most cases) be explained solely on the basis of the action of either electric fields or gas, and now heated gas. Rather, it is the interplay of the two, with field-accelerated clusters colliding with gas, that provides enough collisional heating to remove the last few tightly bound solvent molecules. We note that some recent experimental data on the ESI-MS of brine solutions have appeared in the literature16 that support our findings, observing that the strength of the electric potential applied in the ESI apparatus (ie. soft versus hard ionization conditions) can be the difference between sending either bare ions or hydrated ions into the mass spectrometer. The ionization conditions also affect the ion flux, which must play a role as well. However, it is certainly true that the electric field will be larger in hard ionization conditions. Experimental setups that use heated gas also seem to agree with our data, where very hot gas (800 °C) aids desolvation of ions but only down to the first coordination shell.20 We look forward to further experimental data on these simple charged clusters appearing, which can be directly compared with simulation data. In particular, we would welcome experiments that systematically explore more of the parameters, including the path length from nozzle to the detector, the electric potential, different gases and different gas temperatures, and a greater variety of solvents. In addition to the relevance of our study for the development of new models for ESI-MS mechanisms, the method that we have devised to study such a broad range of cluster sizes provides a new way to look at the basic physical and chemical characteristics of charged ion and water clusters with varying sizes. We have presented a large amount of data on both the energetics and structure of the clusters. Our findings for water clusters support the existing literature on the aqueous solvation

Figure 8. Plot of the coordination number Nc as a function of the cluster size nclust. Horizontal dashed (water) or dotted (methanol) lines represent the bulk values (ref 41 for water; this work for methanol). The dashed orange line shows that the limiting line Nc = nclust − 1 is valid when all solvent molecules are in the coordination shell.

being significantly lower than that for water, a result of the larger size of and reduced hydrogen bonding in methanol. Even in the 40 molecule water clusters, the preference for Cl− to sit at the liquid−vapor interface results in Nc being much lower than the bulk value. We do not probe clusters large enough to see a transition from surface to interior solvation states.31 Turning now to the methanol clusters, it would be tempting to claim that the low Nc for Cl− in methanol, showing only two or three methanol molecules directly bound to the Cl− anion, demonstrates a surface preference for the ion. However, we argue that, as shown in the snapshot in Figure 5, a better interpretation is that instead, two or three independent chains of methanol molecules extend from the anion, with the positions of the anion-bound molecules more or less symmetric around the anion. The data for Na+ clusters are also notable, showing a very gradual decrease in coordination number until 10494

dx.doi.org/10.1021/jp308217q | J. Phys. Chem. A 2012, 116, 10488−10495

The Journal of Physical Chemistry A

Article

(18) Chowdhury, S. K.; Katta, V.; Chait, B. T. Rapid Commun. Mass Spectrom. 1990, 4, 81−87. (19) Manisali, I.; Chen, D. D. Y.; Schneider, B. B. Trends Anal. Chem. 2006, 25, 243−256. (20) Covey, T. R.; Thomson, B. A.; Schneider, B. B. Mass Spectrom. Rev. 2009, 28, 870−897. (21) Vrbka, L.; Mucha, M.; Minofar, B.; Jungwirth, P.; Brown, E. C.; Tobias, D. J. Curr. Opin. Colloid Interface Sci. 2004, 9, 67−73. (22) Chang, T.-M.; Dang, L. X. Chem. Rev. 2006, 106, 1305−1322. (23) Wick, C. D.; Dang, L. X. Chem. Phys. Lett. 2008, 458, 1−5. (24) Caleman, C.; Hub, J. S.; van Maaren, P. J.; van der Spoel, D. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 6838−6842. (25) Dang, L. X.; Chang, T.-M. J. Phys. Chem. B 2002, 106, 235−238. (26) Petersen, P. B.; Saykally, R. J. Annu. Rev. Phys. Chem. 2006, 57, 333−364. (27) Jungwirth, P. Faraday Discuss. 2009, 141, 9−30. (28) Bauer, B. A.; Ou, S.; Patel, S. Chem. Phys. Lett. 2012, 527, 22− 26. (29) Perera, L.; Berkowitz, M. L. J. Chem. Phys. 1991, 95, 1954− 1963. (30) Burnham, C. J.; Petersen, M. K.; Day, T. J.; Iyengar, S. S.; Voth, G. A. J. Chem. Phys. 2006, 124, 024327. (31) Koch, D. M.; Peslherbe, G. H. Chem. Phys. Lett. 2002, 359, 381−389. (32) Cwiklik, L.; Andersson, G.; Dang, L. X.; Jungwirth, P. ChemPhysChem 2007, 8, 1457−1463. (33) Dang, L. X. J. Phys. Chem. A 2004, 108, 9014−9017. (34) Sun, X.; Wick, C. D.; Dang, L. X. J. Phys. Chem. A 2011, 115, 5767−5773. (35) Cabarcos, O. M.; Weinheimer, C. J.; Martinez, T. J.; Lisy, J. M. J. Chem. Phys. 1999, 110, 9516−9526. (36) Pluhařová, E.; Jungwirth, P. Collect. Czech. Chem. Commun. 2008, 73, 733−744. (37) Haughney, M.; Ferrario, M.; McDonald, I. R. J. Phys. Chem. 1987, 91, 4934−4940. (38) Lamoureux, G.; Roux, B. J. Chem. Phys. 2003, 119, 3025−3039. (39) Lamoureux, G.; MacKerell, A. D.; Roux, B. J. Chem. Phys. 2003, 119, 5185−5197. (40) Lamoureux, G.; Harder, E.; Vorobyov, I. V.; Roux, B.; MacKerell, A. D., Jr. Chem. Phys. Lett. 2006, 418, 245−249. (41) Yu, H.; Whitfield, T. W.; Harder, E.; Lamoureux, G.; Vorobyov, I.; Anisimov, V. M.; MacKerell, A. D., Jr.; Roux, B. J. Chem. Theory Comput. 2010, 6, 774−786. (42) Yu, H.; Geerke, D. P.; Liu, H.; van Gunsteren, W. F. J. Comput. Chem. 2006, 27, 1494−1504. (43) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2635−2645. (44) Rahman, A. Phys. Rev. 1964, 136, A405−A411. (45) Soper, A. K.; Castner, E. W.; Luzar, A. Biophys. Chem. 2003, 105, 649−666. (46) Pagliai, M.; Cardini, G.; Schettino, V. J. Phys. Chem. B 2005, 109, 7475−7481. (47) Faralli, C.; Pagliai, M.; Cardini, G.; Schettino, V. Theor. Chem. Acc. 2007, 118, 417−423. (48) Faralli, C.; Pagliai, M.; Cardini, G.; Schettino, V. J. Chem. Theory Comput. 2008, 4, 156−163. (49) French, J. B.; Etkin, B.; Jong, R. Anal. Chem. 1994, 66, 685−691. (50) Cech, N. B.; Enke, C. G. Mass Spectrom. Rev. 2001, 20, 362− 387.

of ions, showing the preference of cations for the interior of the cluster and of anions for the surface. For methanol, our most important new finding concerns the solvation of the Cl− anion. These clusters are highly nonspherical and are dominated by long chains of methanol that emanate from the anion without any significant asymmetry. This chain structure also helps explain the comparative ease with which methanol can be removed from a cluster with Cl− and the prevalence of multiple molecule evaporation events in this system. Our future plans include the extension of our simulations to mixtures of water and methanol. It is well-known that the presence of an organic component like methanol in the solvent is a considerable aid in producing the desired bare ion in ESIMS.50 However, with the exception of the Cl− cluster, we did not observe many qualitative differences between our simulations of water and methanol. It may be necessary to combine the two solvents and observe their interactions directly if we want to understand the importance of different solvent compositions.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +1 (0)613 5332651. Fax: +1 (0)613 5336669. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Richard Oleschuk and Graham Gibson for helpful discussions of ESI experiments, NSERC for funding, and Compute Canada for computing resources. We also thank a referee from our first publication12 for suggesting that we take a closer look at the effects of evaporative cooling in small clusters.

■

REFERENCES

(1) Dole, M.; Mack, L. L.; Hines, R. L.; Mobley, R. C.; Ferguson, L. D.; Alice, M. B. J. Chem. Phys. 1968, 49, 2240−2249. (2) Iribarne, J. V.; Thomson, B. A. J. Chem. Phys. 1976, 64, 2287− 2294. (3) Fenn, J. B.; Mann, M.; Meng, C. K.; Wong, S. F.; Whitehouse, C. M. Science 1989, 246, 64−71. (4) Caleman, C.; van der Spoel, D. Phys. Chem. Chem. Phys. 2007, 9, 5105−5111. (5) Miller, M. A.; Bonhommeau, D. A.; Heard, C. J.; Shin, Y.; Spezia, R.; Gaigeot, M. P. J. Phys.: Condens. Matter 2012, 24, 284130. (6) Luedtke, W. D.; Landmann, U.; Chiu, Y.-H.; Levandier, D. J.; Dressler, R. A.; Sok, S.; Gordon, M. S. J. Phys. Chem. A 2008, 112, 9628−9649. (7) Consta, S. J. Phys. Chem. B 2010, 114, 5263−5268. (8) Chung, J. K.; Consta, S. J. Phys. Chem. B 2012, 116, 5777−5785. (9) Ahadi, E.; Konermann, L. J. Am. Chem. Soc. 2010, 132, 11270− 11277. (10) Ahadi, E.; Konermann, L. J. Am. Chem. Soc. 2011, 133, 9354− 9363. (11) Ahadi, E.; Konermann, L. J. Phys. Chem. B 2012, 116, 104−112. (12) Daub, C. D.; Cann, N. M. Anal. Chem. 2011, 83, 22393−22399. (13) Kebarle, P. J. Mass Spectrom. 2000, 35, 804−817. (14) Nguyen, S.; Fenn, J. B. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 1111−1117. (15) Thomson, B. A.; Iribarne, J. V. J. Chem. Phys. 1979, 71, 4451− 4463. (16) Schröder, D. Phys. Chem. Chem. Phys. 2012, 14, 6382−6390. (17) Li, L. Y. T.; Campbell, D. A.; Bennett, P. K.; Hanlon, J. Anal. Chem. 1996, 68, 3397−3404. 10495

dx.doi.org/10.1021/jp308217q | J. Phys. Chem. A 2012, 116, 10488−10495