J. Phys. Chem. 1995,99, 9974-9985
9974
Transfer of Small Ions across the Water/l,2-Dichloroethane Interface Karl J. Schweighofer and Ilan Benjamin* Department of Chemistry, University of Califomia, Santa Cruz, Califomia 95064 Received: January 24, 1995; In Final Form: April 10, 1995@
Molecular dynamics equilibrium and nonequilibrium simulations are used to examine the mechanism and dynamics of transfer of small ions across the water/l,2-dichloroethaneinterface. The equilibrium calculations provide information about the potential of mean force governing the ion transfer. We show that a twodimensional free energy surface, which is a function of both ion position and the degree of solvation by one of the liquids, is necessary for properly describing the transfer process. Extensive nonequilibrium calculations for the transfer of ions from the organic phase to water, as well as from water to the organic phase, under conditions of different external electric fields, are described. These calculations show that the transfer into the aqueous phase is an activated process that is facilitated by surface roughness and that the activation energy is associated with the breakup of interfacial water hydrogen bonds. The transfer into the organic phase involves significant perturbation to the interface structure and dragging of the first hydration shell of the ion into the organic phase.
I. Introduction The transfer of ions across the interface between two immiscible liquids is a fundamental process in many areas of chemistry and biology. Ion transfer across the interface between two electrolyte solutions, together with electron transfer, is one of the most fundamental charge transfer processes in electrochemistry.’ The transfer of inorganic ions from the aqueous phase to the organic phase is an essential step in the important technique of phase transfer catalysis.* Many cellular functions are regulated by ions, and the interface between two immiscible liquids provides a simplified model of the cell membrane which can be used to study ion transfer and processes linked to ion t r a n ~ f e r . ~Finally, -~ the problem of drug delivery in pharmacology is intimately connected with the process of ion tran~fer.~.’ The fundamental importance of the ion transfer process in electrochemistry has been the impetus for numerous experimental and theoretical studies dating back to the beginning of the century with the pioneering work of Nemst and Riesenfield.8 In recent years the applications of modern electrochemical methods, starting with the work of Gavach and Henry,9 have significantly contributed to our knowledge of the kinetics of the transfer of many types of ions in many different liquid/ liquid systems under a variety of conditions.’Xl0 Despite the impressive progress made and the extensive knowledge collected, our understanding of the mechanism of ion transfer can be described as incomplete at best. The fundamental problems are the lack of knowledge about the microscopic structure of the neat interface and, even more critical, the lack of information about the structure of the ionsolvent complexes at the polarized interface. This situation is beginning to change with the advances of new, molecularly sensitive experimental techniques. Although most of the progress has been limited to the study of other interfacial systems such as the electrolytelelectrode interface,’ I some promising work has been recently reported on both the neat structure and electrochemical processes at the liquidlliquid interface.I2-l6 In the future, techniques are likely to significantly contribute to our understanding of the structure and dynamics of the liquidniquid interface. @Abstractpublished in Advance ACS Abstracts, May 15, 1995.
Insight about the mechanism of ion transfer has usually been obtained by fitting measured current-voltage relationships to theoretical expressions. These expressions are obtained from models which are based on several assumptions regarding the structure of the interface, the distribution of ions at the interface, and the associated inner potential change across the interface. The approaches taken in the literature usually fall into the following three categories presented here in an approximate chronological order: 1. The approximate similarity between the measured chronopotentiograms for ion transfer and those for electron transfer at a metal electrode have led to the use of the standard Marcus theory of electron transfer at a solutiodmetal electrode in deriving a Butler-Volmer type equation for the current-voltage relati~nship.”-’~Although reasonable fit with experimental data could often be achieved by assuming some model for the potential drop across the interface and choosing a value for the activation energy, this approach fails to give any insight into the mechanism of the transfer. In particular, it is hard to see how the concept of the reorganization free energy, which is the central element in Marcus’ theory, has any relevance to the transfer heavy ions across the interface. Later, it was found from surface tension and capacity measurements that the inner potential drop across the interface between water and 1,2dichloroethane (DCE) is very small and independent of the applied potential,20and thus the model lost much of its initial appeal. 2. A more recent approach views the transfer of ions as a diffusion across a mixed solvent layer such that the ion electrochemical potential monotonically varies across the interface.’Sl0 Although reasonable agreement with experimental data can be obtained, the model is based on an unproven assumption (mixed solvent region) and is not free from unknown parameters (free energy of change across the interface) which must be empirically determined. 3. Attempts have also been made to treat the ion transfer process as a chemical reaction assuming some reaction free energy curve as a function of ion position (which is the “reaction coordinate”) and employing transition state theory.2’,22 The fundamental problem with this approach is that the shape of the curve is essentially unknown and must be treated as a free parameter, and more importantly, it is hard to reconcile the
0022-365419512099-9974$09.00/0 0 1995 American Chemical Society
Transfer of Ions across WaterDCE Interface assumption of a well-defined transition state whose location is fixed with the modern view of liquid interfaces (see below). It is clear that what is required in order to resolve some of the issues discussed above is a microscopic molecular description of the liquidlliquid interface and of the process of ion transfer across the interface. Significant advances in integral equation and density functional theories of inhomogeneous provide one potential route to the understanding of the structure of the liquidlliquid interface. However, these approaches have been mainly used in the treatment of neat liquid interfaces (mostly liquidvapor and liquidsolid interfaces), and they are not yet at the level which is necessary to begin addressing the complex issues involved in the behavior of ions at an interface between water and an organic liquid. Some progress has also been made in a statistical mechanical description of the electric double layer,24-26and although this represents an important contribution to one aspect of the ion transfer problem, it lacks information about the solvated ion complex in the two-phase region which is critical to the proper description of the transfer process. Similarly, the modification of continuum models by the use of a nonlocal dielectric response represents an interesting way of incorporating liquid structure.27 However, because it is essentially impossible to experimentally distinguish the nonlocal effects from any specific solvent-solute interactions, the nonlocal response parameter of the media must be treated phenomenologically and cannot replace a rigorous microscopic description. It is, then, not surprising that most of our current knowledge about the microscopic structure of the liquiaiquid interface comes from Monte Carlo and molecular dynamics computer simulations. There have been quite a number of studies on the neat liquidlliquid interface^^^-^^ and on adsorption and reactions of solute molecules at the i n t e r f a ~ e , ~ ~but - ~ Oonly a few treat the problem of solute t r a n ~ f e # ~ .and ~ * ion t r a n ~ f e i ' ~across .~~ the liquidlliquid interface and of water transport through a lipid membrane .45 We have previously considered a simplified molecular dynamics trajectory calculation of the ion transfer rate across the interface between a nonpolar liquid and a polar liquid43and have shown that the rate calculated from the molecular dynamics trajectories agrees with the rate calculated using a solution to a one-dimensional diffusion under the influence of an external field, provided that this external field (and the diffusion constant) are calculated using the microscopic model. The simple model used for the two liquids (dumbbell diatomic molecules with Lennard-Jones intermolecular potentials) made it possible to accurately determine the parameters required in the diffusion model, but it did not give much insight into the mechanism of ion transfer at realistic electrochemical liquiaiquid interfaces. In this paper, we use a recently developed microscopically realistic model of the waterDCE interface32to address some of the fundamental questions about the mechanism of ion transfer across the interface between water and an organic phase. In particular, we focus on the issues of the potential of mean force for the ion transfer, the existence and the nature of the barrier to the transfer, the dynamics of the solvation shell exchange when the ion crosses the interface, and the role of the interface structure. We consider the transfer of small ions from the DCE to the water phase, and vice versa. (We will postpone for future publication the more complex issues associated with the transfer of large ions and of ion pairs.) A preliminary report of the work described in detail below has been published.44 The rest of the paper is organized as follows: In section 11, we present a brief discussion of the neat liquidlliquid interface
J. Phys. Chem., Vol. 99, No. 24, 1995 9975 and the potentials and methods employed. In section 111, we consider the free energy profile for transfer of C1- ion as a function of the distance along the interface normal and introduce the concept of solvation coordinate. The main idea is that the transfer of the ion should be viewed as a motion on a twodimensional free energy surface, which is a function of the distance of the ion along the interface normal and the solvation coordinate. In section IV, we present results of dynamic trajectory calculations which show the utility of this view. This section concentrates on the microscopic picture of the ion transfer, the structure of the solvent shell, and the local perturbation to the interface structure for both the transfer from the water to the organic phase and the reverse process. We conclude with a summary and an outlook for future work in section V.
11. Systems and Methods The liquidlliquid interfacial system consists of two adjacent liquid slabs of 343 water molecules and 108 DCE molecules arranged in a rectangular box of size 21.7 x 21.7 x 70.0 A. The interface between the two liquids is contained in the XY plane at z o A. The water molecules occupy approximately the region Z < 0, and the DCE molecules are mostly in the Z > 0 region. The location of the interface is usually identified as the Gibbs dividing surface which, for the liquidvapor interface of a one-component system, is the planar surface drawn such that the decrease in mass on the bulk side of this surface (because of reduced density) is exactly balanced by the excess mass on the other side of the surface. For two-component systems, the location of the Gibbs dividing surface will depend on the reference liquid chosen. However, for the waterDCE system, the two possible choices result in almost the same location, which is approximately where the density of each liquid is near 50% of the bulk density. This system is created by first equilibrating the water and DCE in separate boxes, with periodic boundary conditions applied in all directions, and letting the box length vary to reproduce the correct bulk density for each liquid. The two boxes are joined along one face, while the total length of the box along the normal to this face is increased to a value that results in only a single liquidlliquid interface. This method of creating the interface results in a water liquid-vapor interface at e-20 A and a DCE liquid-vapor interface at e 2 5 A. These interfaces do little to perturb the bulk liquids or the liquidliquid interface, as discussed in detail elsewhere.32 The system is equilibrated for 200 ps at T = 300 K. Several properties of the system, including the temperature of different modes, the interaction energy between the two liquids, and the density profiles, are monitored to ensure that the system reaches an equilibrium before data collection begins. The time-averaged density profiles over a 100 ps trajectory are shown in Figure 1. These profiles are calculated by dividing the simulation box into 1.0 8, slabs parallel to the X-Y plane. The position of each slab is given with respect to the system center of mass. The number of molecules of a certain type are counted in each slab. This is done for each time, step, and the results are averaged. Each number density is converted to g/cm3 by the appropriate conversion factor. In the interfacial region, the density of each liquid varies smoothly; however, in the bulk there are oscillations about the expected average values. The short simulation time combined with the relatively small system size may account for these oscillations to some extent, especially in the DCE, which has only 108 molecules. Another choice of periodic boundary conditions that avoids the establishment of the two liquidvapor interfaces results in two liquidlliquid interfaces separated by a slab of water and a
Schweighofer and Benjamin
9976 J. Phys. Chem., Vol. 99, No. 24, 1995
I
-30
-15
0
15
z (A,
30
Figure 1. Density profiles of water (solid line) and of 1,2-dichloroethane (dotted line) in contact at 300 K. The liquidkquid interface is near Z = 0 A. Note the two liquid-vapor interfaces at the two ends of the simulation box due to the choice of the periodic boundary conditions, as explained in the text.
TABLE 1: Parameters for the Solvent-Solute Potentialp i-i
0-C1-
(A)
H-C1-
C-C1-
2.14 0.359 -0.410
3.96 0.308 -0.221
E,]
3.55 0.359 0.820
a
U,]= 4~,,[(~,,/r,)’* - (a~l~r,l)61 + qulru.
(A]
(kcaymol) qll ( a 4
C1-C1’ 3.55 0.645 0.221
slab of DCE. This choice is not a desirable one here (unless much larger bulk regions are used), as the transfer of the ions across one interface produces large density fluctuations and could cause the two interfaces to “interact”, thus washing out the bulk region. The intermolecular and intramolecular potentials for the two liquids are standard multisite models which give reasonable properties of the two liquids. For the water, we use a version of the flexible SPC where the intramolecular potential is described by Kuchitsu and M ~ r i n o . ~For ’ the DCE molecule, we use a four-site model in which each one of the two CH2 groups is replaced by a united atom of atomic mass 14. The intermolecular potential is given by a pairwise additive sitesite Lennard-Jones plus Coulomb terms. The intramolecular potential for each DCE molecule is fully flexible and includes one torsion, two bending, and three stretching modes. The parameters for all the terms described above are given elsewhere.32 The water-DCE interaction is calculated by using the standard combining rules for mixtures. It is important to stress that no artificial constraints are placed on the system to ensure solvent immiscibility. The properties of each solvent ensure that the two fluids do not mix and result in an interfacial region that is sharply defined. Density fluctuations normal to the interface give rise to capillary waves that broaden the interfacial region without mixing the solvents. The low mutual solubility of the liquids is largely a result of the strong hydrogen-bonding forces in the water.48 The potentials governing the ion-solvent interactions are a combination of Lennard-Jones and Coulombic terms, and although we have performed simulations on up to five different ions, we mainly focus on the results for a chloride ion whose Lennard-Jones parameters are given in Table 1. These potential energy functions (as well as the liquid-liquid potentials) approximate the many-body nature of the potentials by effective pair interactions. Although there are quite a number of models which explicitly take into account the polarizable nature of the ~ o l v e n t ~ and ~ - ~the ’ ionic s o l ~ t eand ~ ~give , ~ better ~ agreement with experiments (they also have more parameters), it is not clear whether their application will change the basic physical
picture of the ion transfer process. For example, there is evidence that the structure of the interfacial water is only mildly effected by such nonadditive f0rces.5~Nevertheless, the neglect of the polarizable nature of the liquids and the solute is a major approximation, and future studies are planned where the effect of polarizability will be considered. All intermolecular interactions are smoothly switched to zero over an interval of 1 A, at a distance of half the box width, using as a switching function a quintic polynomial with continuous first and second derivatives. This helps minimize the error due to the truncation of the long-range Coulomb interactions. There are well-established methods for including the effect of the long-range Coulomb interactions in simulations of bulk and these methods have been extended to simulations of a liquid near a flat wa1156-57and to simulations of thin membranes5* It is not clear how to adapt these methods to the case of the liquidiquid interface and what the effect would be of including infinite replicas of the interface inhomogeneous charge distribution on the surface density correlations. The integration of the equations of motion is done using the velocity version of the Verlet algorithm55 with a time step of 0.5 fs. The results reported below correspond to canonical ensemble averages at T = 300 K generated by running many constant-energy trajectories with initial conditions that are selected from a Boltzmann distribution. Experimental studies of ion transfer usually involve driving ions across the interface with an extemal potential, and a significant part of the calculations reported below are done in the presence of these fields. The effects of electric fields on the structure and dynamics of the neat interface have been characterized e l ~ e w h e r ewith , ~ ~ the general result that the local structure of the solvents remains unchanged by the fields we employ here. The magnitude of the fields we use, 0.1-0.3 VIA, are close to those that are commonly employed in voltammetric experiments, assuming that the potential drop is across an interface that is on the order of 10 8, (the average thickness of our interface). The transfer of a single ion between two liquids considered in this paper is a nonequilibrium process, because there is no finite ionic concentration to allow for the establishment of constant electrochemical potential throughout the system. One type of calculations reported below represents nonequilibrium molecular dynamics trajectories that enable us to explore the mechanism of the ion transfer process. In these calculations, which are schematically depicted in Figure 1, the ion is equilibrated at the center of the bulk phase and then allowed to cross to the other phase, usually with the help of a static homogeneous electric field. Each trajectory is integrated for 50 ps (or until the ion reaches the other bulk phase) and repeated 10-20 times by selecting different initial conditions. This is done both for the transfer from the DCE to the aqueous phase (the ion starts in the bulk of the DCE phase) and for the reverse process (the ions starts in the bulk of the water phase) and for different values of the extemal electric field. The total length of the calculations is about 5 ns. Some of these calculations were repeated for other ions (Na+, K+, and F-). Much insight can also be gained from equilibrium ensemble averages of local thermodynamics properties, which may be obtained by restricting the ion to a given region (in Z). In this second type of calculations, the region between the two bulk regions is divided into 11 3 8, thick slabs parallel to the interface, whose centers are 2 8, apart (so that there is 1 8, overlap between neighboring slabs). The ion is constrained to a slab by a “window” potential, and statistics about the ion position, the average force on the
Transfer of Ions across WaterDCE Interface ion and its solvation energy, and the structure of the solvation shell are collected from a 400 ps trajectory in each slab. The second type of calculations are done in order to generate enough statistics at every point along the Z coordinate so that the potential of mean force may be obtained. We discuss these equilibrium calculations first in section 111 and the nonequilibrium trajectories in section IV.
111. Equilibrium Calculations 111.1. The Free Energy Profile. We define the local free energy of a single ion in the inhomogeneous system as follows:
where ,8 = l/kT, zi is the position of the ion along the interface normal, 6 is the Dirac delta function, and the ensemble average is over all positions of the fluid molecules subject to having the ion position fixed. A(Z1) - A(Z2) represents the reversible work necessary to bring the ion from point 2 to point 1 . It is important to distinguish this property from the chemical potential of the ion in a system where an equilibrium distribution of ion concentration profile is maintained. In such a case, the chemical potential is constant throughout the system. Under the condition of no extemal field, the quantity defined in eq 1 is simply related to the singlet direct correlation function C(I)(Z).~ The difference A ( Z B I )- A(&), where ZBI and zB2 are positions in the bulk region of the two liquids, is the free energy of transfer of a single ion from the bulk of one liquid to the other. This is a quantity of fundamental importance to diverse areas such as organic synthesis and protein folding. Since, experimentally, one can measure only the free energy of transfer of pairs of ions, certain assumptions are needed in order to split the contribution to individual ions.61 There have also been a few attempts to compute the net free energy of transfer by electrostatic6’ nonlocal continuum dielectric models27 and statistical mechanical analytical methods.61 All of these approaches successfully avoid the actual physical transfer across the interface by taking advantage of the path-independent nature of the free energy of transfer. This, of course, does not apply to the free energy projZe, A(Z), but the net free energy of transfer is a useful check on the validity of the particular model used to compute the profile. Knowledge of the full free energy profile A(Z) is crucial, however, when one considers the actual physical transfer of ions across the interface. Despite its importance, this quantity is basically unknown for realistic systems, in part, because it is experimentally inaccessible. In most of the treatments where this quantity is needed, a mathematical form is assumed that includes free parameters, which are then fitted to reproduce experimental data or estimated using some models. A(Z) may be computed directly from eq 1 by sampling zi in overlapping windows with the help of a biasing potential in order to increase the sampling statistics in the high-free energy region^.^*,^^ This method has been used to study the free energy profile of solutes at different types of interface^.^^.^-^^ As will be discussed in detail below (section 111.2), the wandering nature of the water/ DCE interface and the existence of capillary fluctuations, both of which get more pronounced when ions are present at the interface, present some technical, as well as fundamental, difficulties in implementing this technique for the water/DCE interface. We therefore discuss first a simpler, approximate route to obtaining A(Z). This method is also less computationally demanding but accurate enough to include all the important physical features.
J. Phys. Chem., Vol. 99, No. 24, 1995 9977 The solvation free energy of an ion in a polar liquid may be partitioned into ideal gas, cavity, and electrostatic contributions:
Since the ideal gas contribution (Aldeal)depends only on the temperature, volume, and other constants (independent of the position relative to the interface), it will be omitted in what follows. The partition in eq 2 is based on a thermodynamic path where one first creates the uncharged “ion” with the final value of the Lennard-Jones size (0)and dispersion interaction parameter ( E ) (representing A,,,) and then slowly increases its charge to the final value (which yields Aelec). The cavity term, Acav,represents the reversible work necessary to create a cavity of a given size at a position Z. In creating this cavity, work must be done to change the equilibrium structure of the liquids, although some is gained due to the dispersion interactions between the liquid and the uncharged ion. For small ions in polar liquids, this term is usually much smaller than the electrostatic contribution. The free energy for creating small cavities (diameter less than 3 A) in a dense liquid can be easily and accurately found by determining the size distribution of naturally occurring cavities using molecular dynamics or Monte Carlo simulation^.^^ Somewhat larger cavities appropriate for our need here could still be computed with reasonable accuracy in bulk liquids, but the method becomes increasingly more computationally demanding when one seeks the free energy of atomic size cavity formation as a function of Z-the distance along the interface normal. Instead, we estimate the cavity work by generalizing the scaled particle the~ry~O-’~ to the interfacial system. In the scale particle theory, the actual liquid is modeled using an equivalent hard-sphere liquid of diameter d, and the chemical potential of a hard-sphere solute of diameter D is calculated by computing the probability for a successful insertion of this hard sphere into the hard-sphere l i q ~ i d . ~ ’Our generalization is based on assuming a mathematically sharp interface separating two hardsphere liquids of diameters dl and d2. The probability of inserting a hard sphere of diameter D at a distance Z from the interface is approximated by the product of inserting the appropriate fractions of the cavity in each liquid. The result for the free energy of cavity formation is
(3) where p ~ , ~ ” (=i )-1n(l - vi)
(1 - vi?
+ 3 7 1 , ~+; ~1)/(1 ~ - vi)+ 1/29vi2xi2/
+ viu + vi + vi2 - vi3)xi3/(1 - vi31 (4)
= Di/dj, and v i = (n/6)~idi3.Di is the diameter of a sphere whose volume equals the fractional volume of the cavity of diameter D in liquid i, which depends on the location of the cavity: Di = [‘/2(D- Z)(D 2Z)2]1’3.Z = I/2D corresponds to the cavity being in liquid i and Z = -I/2D in the other liquid. Because we are using a hard-sphere model to represent a Lennard-Jones molecular liquid, the parameters di are effective sizes that can be selected by a variety of ways.73 We have chosen the location of the peak of the liquid’s radial distribution function as the diameter di. Fortunately, variations in di values within a reasonable range do not effect the qualitative shape of A(Z), to be discussed below. The cavity profile for transfer of C1- from water to DCE is shown as a thin line in panel A of Figure 2. The general shape of the curve is in agreement with the qualitative expectation that it is harder to create a large cavity in water than in a chlorinated hydrocarbon fluid. The larger Xi
+
Schweighofer and Benjamin
9978 J. Phys. Chem., Vol. 99, No. 24, 1995
15
’\_
‘
’
’
’
‘
B
0 -20
..-10
0
z (4
fl
I
10
.
20
Figure 2. Free energy of a chloride ion as a function of the distance from the interface. Panel A presents the cavity contribution (dotted line) and the electrostatic contribution calculated using either the molecular dynamics data (solid line) or from a solution of the Poisson equation (dashed line). Panel B shows the sum of the cavity term and the electrostatic term (solid line, molecular dynamics data: dashed line,
Poisson equation). solubility of inert gases in hydrocarbon fluids and in DCE, compared with that in water,6’ is a direct manifestation of this effect.69 The electrostatic contribution to the free energy of the ion near the interface may be obtained by a standard free energy perturbation calculation, where the charge on the ion is increased by a small fraction Aq = Q/N reaching the full charge Q in N steps. The total electrostatic contribution is N- I
where each term in the sum is the electrostatic free energy change at the nth step. U(Z)is the electrostatic potential induced by the liquid(s) at the location Z of the ion, and ( . . J n represents the ensemble average with the ion having a charge ndq. Typically, one needs to use n > 20 with a few tens of thousands of configurations (a few tens of picoseconds) per term (a total of a few hundred picoseconds for each value of Z) to get reasonably convergent results. This does not present a major difficulty, as long as the ion is restricted to be in a window far from the interface. However, when the ion is at the interface, this type of calculation suffers from the same problem mentioned above (and discussed further in section 111.2). To overcome this, we use the fact that, to a very good approximation, U is a Gaussian random variable, and so the ensemble average could be replaced by the integral
u,
where is the average electrostatic potential at the location of the ion when the charge of the ion is ndq. k,, is related to the width of the Gaussian: k, = ‘/2(v- .,)?i Substituting eqs 6 and 7 in eq 5 gives
The advantage of this approximate formula is that short equilibrium calculations are sufficient for achieving convergence of the first and second moments of the distribution P(U). Levy and co-workers discussed this method for a multiple-site solute and show that good results can be obtained even for n = 2. (Thus, the ensemble averages needed involve the initial and final charges only.74) We have used this approach previously to compute the free energy profile for an ion in a simple polar/ nonpolar liquid-liquid interface.43 In that system, it was possible to directly compare this approximate method with the results of umbrella sampling, and good agreement was found. The results of the calculations obtained from a short 50 ps run in each of the windows (corresponding to different Z values) is shown by the thick solid line in panel A of Figure 2. The electrostatic contribution monotonically decreases as the ion is moved from the DCE to water. This feature is well understood as being a direct consequence of the different Born charging energy in the different dielectric media. We also show in panel A of Figure 2, by the dashed line, the numerical solution of the Poisson equation for a charged sphere of diameter D at different locations relative to the sharp interface between two dielectric media whose dielectric constant is determined from independent molecular dynamics calculations of the neat individual liquids (82.5 for flexible SPC and 10 for flexible DCE). The continuum dielectric model results are qualitatively similar to the molecular dynamics calculations. Note that, in a typical application of the continuum electrostatic model, the radius of the charged sphere is a parameter that is normally fit to reproduce the solvation free energy of the ion in the bulk liquid. This procedure will result in different ion sizes in different liquids. The results shown here could be somewhat improved by making the ion’s size dependent on its location. However, the characteristic (and unrealistic) sharp change in the free energy as the ion crosses the interface is a direct result of the discontinuity in the dielectric constant in this model. By choosing a distancedependent dielectric constant E ( Z ) (about which there is no experiment or theory), the continuum model could be made to fit the molecular dynamics data. Panel B of Figure 2 presents the total A , & ) + Aelec(Z)free energy profile. The net free energy of transfer from the molecular dynamics data is about 15 kcallmol, in reasonable agreement with the experimental value of 12.4 kcallmol. The discrepancy reflects mainly the inaccuracy of the simple cavity model. As discussed previously, if one is only interested in the net free energy of transfer, one could calculate it by a thermodynamic cycle which avoids the physical transfer across the interface and obtain more accurate result^.'^ The interesting feature of Figure 2B is a significant barrier of about 3 kcallmol in the “downhill” direction in free energy. Because the existence of this barrier is inferred from the individual approximate calculations of the cavity and the electrostatic contributions, the exact value is somewhat uncertain; but we may assume that this barrier is the result of the higher cost of creating a cavity for the ion in the water phase, which eventually is “paid off’ by the lower (more negative) electrostatic contribution. Note that the dashed line in Figure 2B, which represents the cavity plus the continuum electrostatic contribution, does not show a barrier. This is because the electrostatic contribution from the solution of the Poisson equation changes too fast at the interface (reflecting the discontinuity in the dielectric model). Further support for the fact that the structure of the water is an important factor in creating such a barrier is provided by an earlier study
Transfer of Ions across WaterDCE Interface
J. Phys. Chem., Vol. 99, No. 24, 1995 9979
of ion transfer across the interface between a polar and nonpolar simple In this case, an umbrella sampling calculation of the potential of mean force for the ion transfer shows no barrier in a system where the cavity work contribution is the same in the two liquids. The existence of a barrier to ion transfer between two different liquids has been speculated about before based on rate measurements.*’ The barrier has been attributed to short-range intermolecular forces involved in the process of the switching of the members of the ion’s solvation shell-a necessary step in the transfer process (which we will discuss in detail below). However, it must be stressed that exchange of solvent molecules in the ion’s coordination shell is a process that occurs continuously when an ion moves in bulk homogeneous liquid, yet we do not associate a well-defined potential of mean force and a barrier with this process. (Indeed, in bulk homogeneous media, A(Z) must be a constant because of translation invariance.) Taking into account the fact that it costs more work to break away a water molecule than a DCE molecule from the C1- first solvation shell does not change this conclusion, because we observe a barrier for the “downhill” process Cl-(DCE) Cl-(water). Thus, although the process of solvation shell exchange is extremely important in the transfer process, it is not clear how such a process by itself could produce a welldefined barrier along Z. Our nonequilibrium trajectory calculations described in section IV allow us to resolve these issues, at least for the system discussed in this paper. 111.2. Interfacial Fluctuations. In calculating any equilibrium property of a solute using molecular dynamics computer simulations, one assumes that the ensemble average includes most of the frequently occurring solvent configurations. This presents a significant technical problem in a system like the water/DCE interface because of slowly relaxing surface density fluctuation^.^^ These fluctuations can be described as protrusions of water molecules into the DCE phase, creating water “fingers” whose amplitudes are on average 6 A long, but can be as long as 10 A. Similarly, but to a lesser extent, there are DCE “fingers” which are smaller. Estimates of the time scale for the decay of these fluctuations (in the system of the size considered here) have been given,32and they are in the range of tens of picoseconds. Thus, the 50 ps time scale we use to calculate the electrostatic component of the free energy in the previous section seems reasonable. There is, however, a much more serious technical and fundamental problem, which we discovered by extending the time scale of the simulation while the ion is constrained to be in a given window. During these simulations, a water “finger”, which in the neat waterDCE system would retract, may interact with the ion and become stabilized. Because the ion is constrained in a window relative to the center of mass of the system, the location of the ion relative to the average location of the interface is fixed within the specified range, as it should be. However, the structure of the interface has now changed and become broader. Worse, the electrostatic contribution to the ion’s free energy has now also changed, and if one waits long enough, it becomes almost Z-independent as water molecules “leak” into the region where the ion is located and are able to solvate it! The surface fluctuations are greatest when the ion is nearest the interface and only occur when the ion is on the DCE side of the interface. This is demonstrated in Figure 3, where the 200 ps average density profiles of water and of DCE, in a system where the C1- ion is restricted to a given window, is shown. In panels A and D, no significant effect is evident because the ion is far enough from the interface (in bulk DCE, in panel A, or in bulk water, in panel D). However, when
-
-30 -15 Oe 15 30
z (A)
-30 -15 On 15 30 z (A)
Figure 3. Density profiles of water (solid lines) and of DCE (dotted lines) in a system where a C1- ion is restricted to be in a 3 8, thick slab centered at different locations relative toethe interface. Paneis A-D refer to the window center located at 8.5 A (bulk DCE), 6.5 A (bulk DCE), 2.5 8, (interface), and -7.5 8, (bulk water), respectively. Each density profile is calculated from a 200 ps trajectory.
”
0
100
200 300 (PS)
400
t
Figure 4. Panel A: number of water molecules in the first solvation shell of the C1- as a function of time. The first solvation shell is determined from the first minimum in the ion-oxygen radial distribution function (3.8 A). The ion is confined to a window centered at 10.5 8, (thick solid line), at 8.5 8, (dotted line), and at 6.5 8, (thin line). Panel B: ion-water (0)radial distribution functions for each window averaged over the first 200 ps. Panel C: time-dependent ionwater (0)radial distribution functions for the ion located at the window centered at 10.5 8, as a function of time starting at 200 ps. Each line corresponds to a time average over a 20 ps period.
the ion is constrained to be in the DCE phase at 6.5 A away from the interface (panel B) or at 2.5 A from the average interface (panel C), there i s a visible spread in the water density profile. Figure 4 demonstrates the problem even more clearly by looking at the ion’s coordination shell and for longer times. In panel A, the number of water molecules in the first coordination shell is shown as a function of time for the ion restricted to three different windows. Note that after 400 ps, in all the windows, the first coordination shell is completely replaced by water molecules! In panel B, we show the ionoxygen radial distribution function averaged over the first 200 ps for all the windows. The functions in the region Z -= 0 are
Schweighofer and Benjamin
9980 J. Phys. Chem., Vol. 99, No. 24, 1995 essentially those of C1- in bulk water. Those in the region Z > 0 (DCE phase) show that, as one gets closer to the interface, water molecules begin to form a first and second solvation shell. Note that, as is clear from panels A and B, during the first 200 ps, the water did not reach the ion located in the farthest window (10.5 A into the DCE phase), but the first solvation shell for this case gradually builds to essentially a bulk water-like structure by 400 ps (panel C). An animated “movie” of this trajectory shows a small cluster of water molecules surrounding the ion in bulk DCE, with a chain of several water molecules extending from the bulk water into the DCE. This is one of the most dramatic examples we have seen of a water “finger”. Besides the technical difficulties that the above phenomena present when one wants to use an umbrella sampling to calculate the free energy profile, there is a fundamental issue that must be addressed. One first should realize that the interactions of interfacial water with ions located up to 10 .,& away should be viewed as an artifact. There are experimental data supporting a water profile that is up to 10 .,& thick in a waterloil system,I2 and small ions should be able to strongly interact with such interfacial water molecules. The results presented in Figures 3 and 4 suggest that the characterization of the state of the ion at the interface, in terms of a single coordinate Z, is incomplete and can be misleading when the ion is close to the interface. Constraining the ion to have a certain value of Z is not sufficient to describe its state at the interface. For example, if we wait long enough, two different ion positions relative to the average location of the interface will have almost the same solvation shell structure and very similar interaction energy with the two liquids. The proper description of the state of the ion is, of course, a key ingredient in any attempt to understand the mechanism of the ion transfer, and we must quantify this state better than just specifying a value of Z. To explicitly take into account the fact that different solvation environments could be realized almost independently of Z (at least near the interface), we use the electrostatic potential induced by each liquid at the location of the ion as a quantitative measure of its solvation by the corresponding liquid. This enables us to define a solvation coordinate as follows:
(9)
where Vi is the electrostatic potential induced by liquid i at the location of the ion, and upq is the equilibrium value of this quantity when the ion is solvated in the bulk of liquid i. These electrostatic potentials and thus the solvation coordinate are generalized solvent coordinates because they are functions of the positions of all solvent molecules in the system. We note that these two coordinates may not be independent because we expect that, as the ion is desolvated by one liquid, it is more likely to be solvated by the second liquid. Still, it will be useful to consider both parameters as a measure of the solvation state of the ion. In the rest of the paper, reference to s should be taken to be either one of these two parameters unless stated otherwise. This solvation coordinate is similar to the one used in simulations of heterogeneous electron transfer reaction^^^-^^ where a proper account of solvent configurations different from the equilibrium one is of prime importance.80 111.3. A Two-DimensionalFree Energy Surface. In view of the discussion in the previous section, the free energy profile shown in Figure 2B should be considered a projection of a twodimensional free energy surface on the Z coordinate, the other coordinate being the degree by which the ion is solvated by
mm
5 4 $3 2 1 \ I 0 -0.2 0.2 0.6 1 1.4 -0.2 0 2 0.6
1
1.4
DCE Figure 5. Free energy associated with the fluctuations in the solvent coordinate for an ion located at different positions relative to the interface. The panels on the left (A-D) are for water, and those on the Pght (E-H) are for DCE. Panels A and E are for theion in the 8.5 A window; panels B and F, 6.5 8, window; C and G , 4.5 A window; D and H (solid lines), 2.5 8, window; D (dotted line), -7.5 8, window (bulk water).
one of the liquids. Each point in Figure 2B gives the free energy of the ion located at a distance Z from the interface, while its solvation shell structure approximately reflects a sharp interface between the two liquids. On the other hand, the calculations described in Figures 3 and 4 are done for an ion constrained to be at a given Z, while the system is allowed to explore all values of s and find a minimum (in s). In this case, the projection along Z of the free energy will look very different. It is not clear which surface one should use if one wants to interpret the process of the ion transfer as a motion on a one-dimensional free energy surface (or potential of mean force), nor it is obvious if such a reduction of the many-body dynamics is at all possible. An important issue to consider in this respect is the time scale of the ion transfer in relation to the time scale that governs surface fluctuations. This will be addressed later when we discuss the nonequilibrium ion transfer trajectories. Regardless of the suitability of the free energy surface in interpreting the dynamics, it is useful at this point to consider its shape along the s coordinate. Figure 5 presents the free energy associated with fluctuations in the solvent coordinates S H ~ Oand SDCE for different locations of the ion. The free energy is calculated from the equilibrium distribution of the coordinates SDCE and s H 2 0 . For each location (Z) of the ion, this free energy represents the minimum work necessary to change the configuration o f the solvent from the equilibrium configuration around the ion. In general, as the ion is coming closer to the interface (top to bottom panels), the location of the DCE minima shifts from 1 (bulk DCE) to a smaller value, while that of the water shifts from small values to near 1 . This simply reflects the desolvation by the DCE and the simultaneous solvation by water. When the ion is in bulk DCE, the DCE contribution (panel E) is nearly a parabola, which is consistent with the Gaussian assumption discussed earlier. In this case, the water-induced electrostatic potential at the location of the ion is randomly fluctuating around zero (panel A). In panels B and F, the ion is near the interface, and the free energy curve exhibits a double minimum representing two relatively stable and different DCE and water configurations
Transfer of Ions across WaterDCE Interface
J. Phys. Chem., Vol. 99, No. 24, 1995 9981
15
$
-5 -15
v 0.5
0
-25
0
$ -15-5 -25
-55 -15
-25
50
K]ma$ CI
0
0 15
$
25
25
50
0
25
50
$l:
2
v 0.5
0
25 t (PSI
50
0
25
50
t (PS)
Figure 6. Time-dependent ensemble averages (20 independent trajectories of 50 ps each) of the ion position (left panels) and the solvent coordinate (right panels) for “downhill” transfer of CI- from DCE to water under the influence of different constant electric fields, as indicated on the right margin. The solvent coordinate for water and DCE is defined in the text.
around the ion. The barrier in s separating them needs to be surmounted during the ion transfer. This can be done by a proper solvent fluctuation at fixed Z or by moving the ion closer to the interface (panels C and G), where the barrier is much smaller. As the ion approaches the water phase, the DCE curves shift to zero while those of water attain nearly the bulk behavior (panels D and F). In panel D, the dashed line is for C1- in bulk water. The free energy cuts shown in Figure 5 will be useful in the interpretation of the nonequilibrium trajectories we discuss next.
IV. Nonequilibrium Calculations We turn now to the examination of the nonequilibrium ion transfer trajectories. In this case, there is no external potential constraining the ion to a slab, and so it is free to diffuse and cross the interface. Our goal in conducting these calculations is to provide a molecular picture of the events that accompany the transfer in terms of the location of the ion, its interaction with the liquids, and the structure of the solvation shell. These nonequilibrium calculations can thus provide clues to the mechanism of the transfer across the interface. This section is divided into the discussion of the transfer from the DCE to the water phase (“downhill”) and the transfer from the water to the DCE phase (“uphill”). In each section, we describe the trajectories for different external electric fields, focusing first on the time-dependent ion location in the r-s plane and then on the time-dependent structure of the solvation shell. IV.l. Ion Transfer from DCE to Water. This process is examined using three different external electric fields: 0 (no field), 0.1, and 0.2 VIA. The average location of the ion and the average value of the normalized solvent (water and DCE) coordinates as a function of time are shown in Figure 6. Under the conditions of no driving force due to an external field (top two panels), the average location of the ion is almost time-independent, reflecting a simple diffusion process with a very small drift. Of course, during the time interval of 50 ps, individual members of the ensemble do change location. Some get close to the interface quite rapidly, but then spend the rest of the time near the interface, and some diffuse away from the
Figure 7. Time-dependent ensemble averages for several structural properties of the ion hydration shell during a “downhill” transfer (DCE -water) with ]El = 0.1 V/A. Upper left: coordination number. Upper right: CI--water oxygen radial distribution function. Lower left: the probability distribution for the angle (19)between the water dipole @) and the CI--water oxygen vector (d). Lower right: the probability distribution for the angle (4) between the water H-H vector and the plane defined by the vectors p and d. Each time slice in the 3-dimensional plots represents an average over 5 ps interval and over all the trajectories. Note that the ion begins to significantly interact with the water only toward the end of the 50 ps period, so only the last 20 ps of the 50 ps period is shown.
interface. None of the trajectories in the zero field case cross the interface. This situation is independent of the initial position of the ion as long as it starts farther than 5 8, away from the average location of the Gibbs dividing surface. This result supports the existence of a barrier in the “downhill” direction, in agreement with the conclusion drawn from Figure 2. Without such a barrier, 50 ps would be more than enough to bring the ion across the interface or at least show a drift due to the significant ( x 10 kcaymol) driving force. The two middle panels in Figure 6 show that an external electric field of 0.1 Vl8, is still not enough to drive the C1- ion from bulk DCE to water on the 50 ps time scale of the trajectories. Most of the trajectories reach the outer region of the interface (6 A), but only one trajectory makes it across the aperoximate location of the barrier. This suggests that the 0.1 VIA is almost enough to “flatten” the barrier. This could be used to put a lower limit of e4 f 1 kcaVmol on the height of the barrier. The bottom two panels of Figure 6, on the other hand, demonstrate that doubling the field strength to 0.2 VIA increases the probability of crossing to 100%. Complete solvation by water occurs in an average time of 30 ps. This can be used to put an upper limit of 8 % 2 kcal/mol on the barrier height. It is interesting to note that although the DCE solvation coordinate decreases monotonically to zero, the water solvation coordinate shows values greater than 1 when the ions begin to lose the last trace of interaction with the organic phase, indicating a tighter hydration shell than the equilibrium one attained in bulk water (see below). Notewithstanding this, one may note that S H ~ O SDCE sz 1. This suggests that the ion tries to maximize its solvation environment as it leaves one liquid and enters the other liquid. More insight into the transfer process is obtained by examining the structure of the C1- solvation shell for the trajectories described above. Since no interaction with water is evident in the zero-field case, we present results for the 0.1 VIA (Figure 7) and 0.2 VIA (Figure 8). The top left two-dimensional plot
+
Schweighofer and Benjamin
9982 J. Phys. Chem., Vol. 99, No. 24, 1995
0
25
50
-
Figure 8. Same as in Figure 7 for a “downhill”transfer (DCE water) with IEl = 0.2 V/A. Note that after 30 ps the system reached equilibrium, so, for greater clarity, only the first 30 ps of data are shown in the 3-D plots.
shows the time-dependent coordination number. It closely resembles the plot of the solvation coordinate in Figure 6, showing that the full coordination number of 6 for C1- in water is formed under the influence of the higher driving force (Figure 8). Note that, during the time interval in which the water solvation coordinate is larger than 1 (mentioned at the end of last paragraph), there is evidence of a seventh water molecule in the first coordination shell. The top right panel shows the time-dependent oxygen (water)-Cl- radial distribution function. In the lower field case (Figure 7; note that we show only the last 20 ps of the trajectories), the first peak of this function grows monotonically up to a value that is about 15% of the bulk value at a distance r, which is slightly shorter than the one for C1- in bulk water, and no second peak is apparent. This plot and an examination of the individual trajectories indicate that in one or two of the trajectories a full solvation shell is formed, but the ion plus its hydration shell is surrounded by DCE molecules (and so there is no second peak). Examination of the animated trajectories shows that the ion is initially interacting with a water “finger” that extends out of the water phase into the DCE phase. (A picture of this was provided in ref 44.) For the upper field case (Figure 8), we show only the first 25 ps of the 50 ps trajectories. The full first and second peaks are evident by this time. It is interesting to note that while the size of the first peak grows monotonically with time, the location and width of the first peak are nearly time-independent. This suggests that the time scale over which a water molecule joins the first solvation shell is much shorter than the 5 ps time interval we use to average the radial distribution function. So, as soon as one water molecule interacts with the ion, the rest follow very quickly. This “Trojan horse” picture is seen directly in the animated trajectories. The total time for formation of the full solvation shell is about 15 ps, which can be viewed as the lifetime of the ion with a mixed solvation shell. The bottom left panel presents the time-dependent probability distribution for the angle between the 0-C1 vector and the water dipole. The equilibrium distribution for C1- in bulk water is broad and centered at about 45”. The result for the low-field case (Figure 7) shows that, during the time that the solvation shell is formed around the single ion that is able to cross the interface, the distribution peak shifts from near zero starting at t = 30 ps to the final value. This indicates that the initial approach of the first water molecule to the ion is along the water
dipole, with both hydrogens pointing toward the ion. Examination of the animated trajectory shows that this water molecule is hydrogen bonded to only one water molecule (compared with an average number of 2-2.5 for water at the neat water/DCE i n t e r f a ~ e ~ Additional ~). water molecules that enter the solvation shell bring the most probable orientational angle to the equilibrium value. This behavior is also apparent with the higher field (Figure 8). Note that because the orientational probability distribution is normalized by the number of the water molecules in the first shell, the total area under the curve is constant: J:P(8) sin 8 d 8 = 1. Thus, as the peak width increases (reflecting a less tight distribution), its height decreases, but the effect is not large. In the bottom right panel of Figures 7 and 8 we consider the time-dependent probability distribution of the angle q5 between the H-H vector and the normal to the plane defined by the water dipole and the 0-C1- vectors. For C1- in bulk water, the distribution peaks at 90”. Taken together with the distribution of the angle 8, this indicates that, in bulk water, the O-HC1- are only 7” off the linear configuration, and so only one hydrogen bond is “used up” when solvating the ion. The timedependent orientational distribution for the low external field case shows (Figure 7) that initially the H-H vector is broadly distributed between 0” and 180”. This is consistent with the observation that initially the C1- is interacting with both hydrogens. (Rotation of the H-H vector around the water dipole when the C1- ion is along this dipole does not change the H-C1 distance.) This is also the situation for the highfield case. The detailed picture of the molecular events that accompany the ion transfer suggests that the rate-limiting step is the ability of the C1- ion to find a water molecule which, although attached to the water phase, has one or both of its hydrogen atoms available for “bonding” with the ion. Loosely connected water molecules are spontaneously formed due to thermal motion and constitute a roughening of the surface, or capillary fluctuations. We suggest that the barrier to the ion transfer is a combination of the entropic factor associated with matching the ion with a water “finger” plus an energetic component which comes from the need to break some of the water-hydrogen bonds and to create such a “finger”. (Using macroscopic terminology, this process requires free energy because work must be done to increase the water surface area.) We showed that the actual process of ion solvation shell replacement is quite rapid and thus not significantly contributing to this barrier. It is interesting to contrast these results with the results of Hayoun, Meyer, and Turq,4* who found a barrier for a transfer of an atomic solute between two identical Lennard-Jones liquids. They attributed this barrier to the process of solvent shell exchange. However, in that work, the two identical liquids were made immiscible by repulsive pair interactions between atoms which belong to “different” liquids. Thus, it is likely that the barrier in their case results from the instability of the mixed solvation shell due to the repulsion between its members. Indeed, no such barrier was found in our earlier study where two Lennard-Jones diatomic liquids were kept from mixing by an external field.43 It is not clear whether instability of a mixed solvation shell due to unfavorable interactions between its members is an important effect in real systems. In the water/DCE system (and several others watedorganic liquid systems), the water-DCE interactions are more favorable (by 5 kcal/mol) than the DCE-DCE interactions, and thus the model of repulsive (or reduced attractive) interactions between the two liquids is not appropriate. IV.2. Ion Transfer from Water to DCE. The transfer of C1- from water to DCE must be done with the help of an
Transfer of Ions across WaterDCE Interface
J. Phys. Chem., Vol. 99, No. 24, 1995 9983
15
1
0.5 -15
0
25
50
n 0
25
50
-
Figure 10. Same as in Figure 8 for an “uphill” transfer (water DCE) with /El = -0.1 V/A. Note that the ion remains solvated in water t (PSI Figure 9. Same as in Figure 6 for the “uphill” transfer from water to DCE with the help of an electric field of different intensities, as indicated on the right margin.
external electric field of an opposite polarity to the one used in the previous section, in order that the intrinsic large difference of the free energy of the ion in DCE and in water be overcome. Figure 9 presents the average trajectory of the ion under the influence of an electric field of intensities -0.1, -0.2, and 0.3 VIA. The top two pairs of panels of Figure 9 shows that the two lowest field intensities are insufficient to drive the ion to the DCE phase on the time scale of the simulations. After a period of linear increase in (z), the ion gets stuck at the interface. The external electric field is able to modify the electrochemical potential of the ion in an expected way: the interface becomes the lowest free energy region. This is easy to see by adding a linearly decreasing function of 2 with a slope -eJEl to the potential of mean force of Figure 2B. The -0.3 VIA electric field is able to drive the ion to the DCE phase, as is clear by the value of the solvation coordinate SDCE or S H ~ O . Note that a specification of the location Z alone is not sufficient to describe the state of the ion. It is also increasing to note that here, as in the transfer “downhill”, the sum SDCE f S H ~ Ois nearly constant. Figures 10-12 describe the structure of the ion solvation shell as it attempts to cross the interface. They are arranged in the same way as the plots for the transfer “downhill”. As before, the coordination numbers mirror the behavior of the solvation coordinate and provide a better gauge to the state of the ion than does its position. The behavior of the 0-C1- radial distribution functions and the orientational probability distribution functions at the two lowest fields (Figures 10 and 11) reflects the fact that the ion’s first hydration shell is essentially unaffected, as expected. The radial distribution function shows that the ion’s second shell is, however, significantly diminished. Figure 12 (the highest-field case) shows that, as the ion begins to lose its hydration shell (the first peak of the radial distribution function gets smaller), the location of the peak is unchanged, and the orientational probability distributions are almost unaffected. The ion keeps the structure of its aqueous solvation shell as intact as possible. Note that the probability distribution for the angle a shows a slight narrowing down, reflecting a tighter hydration shell. Examination of the animated trajectories for the high-field case shows that the external field drives the ion to the bulk of
even though its average position Z (see Figure 9) indicates that it is at the interface. Each time slice in the 3-dimensional plots represents an average over 10 ps interval and over all the trajectories.
0
10
20
Figure 11. Same as in Fi ure 10 for an “uphill” transfer (water DCE) with /El = -0.2 V / f .
0
10
-
20
Figure 12. Same as in Figure 11 for an “uphill” transfer (water DCE) with IEl = -0.3 V/A.
-
the DCE while most of the first hydration shell is still intact. The interface gets wider, and there is a water “finger” connecting the ion to the bulk water. This is almost a playback of the
9984 J. Phys. Chem., Vol. 99, No. 24, 1995 microscopic events that occur in the “downhill” transfer, and it is consistent with the reversible nature of the transfer process, a fact well documented experimentally.8’ Unfortunately, because of the small size of the system, as the ion is further pushed into the DCE phase and the water “finger” grows at the expense of the size of the bulk water, the interface becomes unstable. So simulations in a larger system (which are under way) are necessary to observe the complete dehydration of the ion in the high-external field case. Discussion of the size effects and stability of the neat water/DCE system to electric fields of varying intensities can be found e l s e ~ h e r e . ~ ~
V. Conclusions The equilibrium and nonequilibrium molecular dynamics calculations discussed in this paper help shed light on several outstanding problems regarding the process of electrochemical ion transfer across the interface between two immiscible liquids. We have demonstrated that, although the transfer of a single C1- ion from the organic phase to water has a significant intrinsic driving force, the kinetics of the process suggests that a barrier of a few kcal/mol exists in this “downhill” direction. Approximate free energy calculations support this conclusion. We presented evidence that this barrier is directly related to the unique water structure at the interface; specifically, the first water molecule which begins the hydration process has most of its interfacial hydrogen bonds broken. We also found that the transfer process is facilitated by interface roughness. Water capillaries, or water “fingers”, naturally have loosely connected water molecules, and they are found to be a precursor for the transfer process. Examination of the reverse process is also possible with the halp of an external electric field driving the ion “uphill”. We show that the ion tends to keep its first hydration shell during the crossing of the interface and that the interface gets perturbed. The existence of a water “finger” connecting the ion to the bulk aqueous phase is consistent with the reversible nature of the process. In attempting to describe the transfer of the ion as a motion on a one-dimensional reaction coordinate (position of the ion along the interface normal), one must note that the time scale for the ion solvation shell reorganization is comparable to the time scale of interface capillary fluctuations and to the time it takes the ion to cross the interface. Since the potential of mean force at a given Z is significantly dependent on how many molecules of each solvent are present in the solvation shell, it is not clear that a one-dimensional, position-dependent reaction coordinate is a useful reduction of the many-body dynamics. We have shown the utility of a solvation coordinate, and thus a two-dimensional description of the potential of mean force, as a way to better reflect the state of the ion while it is in the interface region. Although we did not show the data, some of the calculations reported here have also been done for Na+, K+, and F- in order to examine the effect of charge and size on the transfer process. The major conclusions stated above about the role of surface roughness and the ability of the small ions to keep their first hydration shell as they are driven to the organic phase hold for these other ions. Although the first step in the hydration of positive ions involves a water molecule with the oxygen pointing toward the cation, only one of the two hydrogen bonds exists, which is similar to the situation found for the hydration of the anions. Unfortunately the statistics collected to date are not accurate enough to determine the effect of ion charge on the ion transfer rate. More detailed examinations of the effect of ion charge and size are needed. We conclude by mentioning two other topics of importance to the majority of experiments on electrochemical ion transfer,
Schweighofer and Benjamin which have not been discussed in this paper but are currently under investigation. First, many of the experimentally studied systems involve transfer of ions that are bulkier than the ones we considered thus far (e.g., tetramethylammonium, N(CH&+). The transfer of these ions involves little or almost no intrinsic driving force. Thus, longer simulations on larger liquidhquid systems are required. In such systems, the ion will spend significantly more time at the interface region, and some aspects of the mechanism discussed above will most likely be different. For example, the binding of water molecules to the ion will not be as dramatically favored over that of the organic molcules, as has been the case in the systems discussed in this paper. Second, the transfer of small ions (such as Na+ and C1-) with large free energies of transfer could be facilitated by the formation of ion pairs (NaCl), especially at moderate to high ionic concentrations. The mechanism of ion association at the water/DCE interface, which is an important part of the transfer process, is currently under investigation.
Acknowledgment. This work was supported in part by NSF Grant CHE92-21580 and the Senate of the University of California, Santa Cruz. References and Notes (1) Girault, H. H. J.: Schiffrin, D. J. In Electroanalytical Chemistry; Bard, A. J., Ed.; Dekker: New York, 1989; p 1. (2) Starks, C. M.; Liotta, C. L.;Halpem, M. Phase Transfer Catalysis; Chapman & Hall: New York, 1994. (3) Honig, B. H.; Hubbeli, W. L.; Flewelling, R. F. Annu. Rev. Biophys. Biophys. Chem. 1986, 15, 163. (4) Gennis, R. B. Biomembranes; Springer: New York, 1989. ( 5 ) Shao, Y.; Girault, H. H. J . Electroanal. Chem. 1990, 282, 59. (6) Knepp, V. M.; Guy, R. H. J . Phys. Chem. 1989, 93, 6817. (7) Arai, K.; Ohsawa, M.; Kusu, F.: Takamura, K. Bioelecrrochem. Bioenerg. 1993, 31, 65. (8) Nemst, W.; Riesenfield, E. H. Ann. Phys. 1902, 8, 600. (9) Gavach, C.; Henry, F. C.R. Acad. Sei. 1972, C274, 1545. (10) Hanna, G. J.; Noble, R. D. Chem. Rev. 1985, 85, 583. (1 1) Spectroscopic and Diffraction Techniques in Interfacial Electrochemistry; Gutierrez, C., Melendres, C., Eds.; Kluwer Academic Publishers: Dordrecht, 1990; NATO AS1 Series Vol. 320. (12) Lee, L. T.; Langevin, D.; Farnoux, B. Phys. Rev. Lett. 1991, 67, 2678. (13) Wirth, M. J.; Burbage, J. D. J. Phys. Chem. 1992, 96, 9022. (14) Higgins, D. A.; Com, R. M. J. Phys. Chem. 1993, 97, 489. (15) Kott, K. L.;Higgins, D. A,; McMahon, R. J.; Com, R. M. J . Am. Chem. Soc. 1993, 115, 5342. (16) Conboy, J. C.; Daschbach, J. L.;Richmond, G. L. J . Phys. Chem. 1994, 98, 9688. (17) Gavach, C. J . Chim. Phys. 1973, 70, 1478. (18) Gavach, C.; D’Epenoux, B. J. Electroanal. Chem. 1974, 55, 59. (19) Kuznetsov, A. M.: Kharkats, Y. I. In The Interface Structure and Electrochemical Processes at the Boundary Between Two Immiscible Liquids; Kazarinov, V. E., Eds.; Springer: Berlin, 1987; p 11. (20) Kakiuchi, T.; Senda, M. Bull. Chem. Soc. Jpn. 1983, 56, 1753. (21) Gurevich, Y. Y.; Kharkats, Y. I. J. ElecrroanaL Chem. 1986, 200, 3. Samec, Z.; Kharkats, Y. I.; Gurevich, Y. Y. Ibid. 1986, 204, 257. (22) Simonin, J. P.; Moreau, M. Mol. Phys. 1990, 70, 265. (23) Fundamentals of Inhomogeneous Fluids; Henderson, D., Ed.: Marcel Dekker: New York, 1992. (24) Carnie, S. L.: Torrie, G. M. Adv. Chem. Phys. 1984, 56, 141. (25) Torrie, G. M.; Kusalik, P. G.; Patey. G. N. J . Chem. Phys. 1989, 91, 6367 and references therein. (26) Blum, L. Adv. Chem. Phps. 1990, 78, 171. (27) Kornyshev, A. A,: Volkov, A. G. J . Electroanal. Chem. 1984, 180. 363. (28) Linse, P. J. Chem. Phps. 1987, 86, 4177. (29) Gao, J.; Jorgensen, W. L. J . Phys. Chem. 1988, 92, 5813. (30) Meyer, M.; Mareschal, M.: Hayoun, M. J . Chem. Phys. 1988, 89, 1067. (31) Carpenter, I. L.; Hehre, W. J. J . Phy.7. Chem. 1990, 94. 531. (32) Benjamin. I. J. Chem. Phys. 1992, 97, 1432. (33) Vanburren, A. R.; Manink, S. J.: Berendsen, H. J. C. J . Phys. Chem. 1993, 97, 9206. (34) Smit, B. Phys. Rev. A 1988. 37, 3431.
Transfer of Ions across WaterDCE Interface (35) Smit, B.; Schlijper, A. G.; Rupert, L. A. M.; van Os, N. M. J . Phys. Chem. 1990, 94, 6933. (36) Benjamin, I. J . Phys. Chem. 1991, 95, 6675. (37) Benjamin, I. In Sfructure and Reactivity in Aqueous Solution: ACS Symposium Series 568; Cramer, C. J., Truhlar, D. G., Eds.; American Chemical Society: Washington, DC, 1994; p 409. (38) Benjamin, I. Chem. Phys. 1994, 180, 287. (39) Benjamin, I. In Reaction Dynamics in Clusters and Condensed Phases; Jortner, J., Levine, R. D., Pullman, B., Eds.; Kluwer: Dordrecht, The Netherlands, 1994; p 179. (40) Michael, D.; Benjamin, I. J. Phys. Chem. 1995, 99, 1530. (41) Hayoun, M.; Meyer, M.; Mareschal, M. In Chemical Reactivity in Liquids; Ciccotti, G., Turq, P., Eds.; Plenum: New York, 1987. (42) Hayoun, M.; Meyer, M.; Turq, P. J . Phys. Chem. 1994, 98, 6626. (43) Benjamin, I. J . Chem. Phys. 1992, 96, 577. (44) Benjamin, I. Science 1993, 261, 1558. (45) Marrink, S.-J.; Berendsen, H. J. C. J. Phys. Chem. 1994, 98,4155. (46) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hennans, J. In Intermolecular Forces; Pullman, B., Ed.; D. Reidel: Dordrecht, 1981; p 331. (47) Kuchitsu, K.; Morino, Y. Bull. Chem. Soc. Jpn. 1965, 38, 814. (48) This can easily be shown by artificially breaking the hydrogen bonds within the water using a strong electric field, which causes the liquids to start mixing. When the electric field is removed, the interface is restored. (49) Sprik, M.; Klein, M. L. J . Chem. Phys. 1988, 89, 7556. (50) Ahlstrom, P.; Wallqvist, A,; Engstrom, S.; Jonsson, B. Mol. Phys. 1989, 68, 563. (51) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141. (52) Dang, L. X.; Rice, J. E.; Caldwell, J.; Kollman, P. A. J . Am. Chem. Soc. 1991, 113, 2481. (53) Dang, L. X. J . Chem. Phys. 1992, 97, 2659. (54) Motakabbir, K.: Berkowitz, M. Chem. Phys. Lett. 1991, 176, 61. ( 5 5 ) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987. (56) Fincham, D.; Heyes, D. M. Chem. Phys. 1983, 78, 425.
J. Phys. Chem., Vol. 99, No. 24, 1995 9985 (57) Spohr, E. J . Phys. Chem. 1989, 93, 6171. (58) Hautman, J.; Klein, M. L. Mol. Phys. 1992, 75, 379. (59) Schweighofer, K. J,; Benjamin, I. J . Electroanal. Chem., in press. (60) The proper definition of local thermodynamic functions in an inhomogeneous fluid is somewhat problematic. See: Rowlinson, J. S.; Widom, B. Molecular Theory of Cnpillariry; Clarendon: Oxford, 1982; Chapter 4. (61) Marcus, Y. Ion Solvation; Wiley: New York, 1985; Chapter 6. (62) Torrie, G. M.; Valleau, J. P. J . Comput. Phys. 1997, 23, 187. (63) Chandler, D. Introduction to Modem Statistical Mechanics: Oxford University Press: Oxford, 1987; p 169. (64) Wilson, M. A,; Pohorille, A. J . Chem. Phys. 1991, 95, 6005. (65) Pohorille, A,; Benjamin, I. J . Chem. Phys. 1991, 94, 5599. (66) Rose, D. A,; Benjamin, I. J . Chem. Phys. 1991, 95, 6856. (67) Benjamin, I. J. Chem. Phys. 1991, 95, 3698. (68) Pohorille, A,; Benjamin, I. J . Phys. Chem. 1993, 97, 2664. (69) Pohorille, A,; Pratt, L. R. J . Am. Chem. Soc. 1990, 112, 5066. (70) Reiss, H.; Frisch, H. L.; Lebowitz, J. L. J . Chem. Phys. 1959, 31, 369. (71) Stillinger, F. H. J . Solution Chem. 1973, 2, 141. (72) Pierotti, R. A. Chem. Rev. 1976, 76, 717. (73) Andersen, H. C.; Chandler, D.; Weeks, J. D. Adv. Chem. Phys. 1976, 34, 105. (74) Levy, R. M.; Belhadj, M.; Kitchen, D. B. J . Chem. Phys. 1991, 95, 3627. (75) Schweighofer, K. J.; Benjamin, I. Manuscript in preparation. (76) Straus, J. B.; Voth, G. A. J . Phys. Chem. 1993, 97, 7388. (77) Rose, D. A,; Benjamin, I. J . Chem. Phys. 1994, 100, 3545. (78) Xia, X.; Berkowitz, M. L. Chem. Phys. Lett. 1994, 227, 561. (79) Rose, D. A,; Benjamin, I. Chem. Phys. Lett. 1995, 234, 209. (80) Marcus, R. A. J . Chem. Phys. 1956, 24, 979. (81) For example, the location of the peaks in the cyclic voltammogram for a typical ion transfer process is independent of the sweep rate, and the heights of the anodic and cathodic peaks are equal.
JP950232Y
