J. Phys. Chem. 1994,98, 6626-6632
6626
Molecular Dynamics Study of a Solute-Transfer Reaction across a Liquid-Liquid Interface Marc Hayow* and Madeleine Meyer Laboratoire des Solides Irradih, CEA-CEREM-DTM, CNRS URA No. 1380, Ecole Polytechnique, 91 128 Palaiseau Cedex, France
Pierre Turq Laboratoire d’Electrochimie, Universitb Pierre et Marie Curie, 8 rue Cuuier, 75005 Paris, France Received: December 14, 1993; In Final Form: April 21, 1994’
The aim of this study is to model a solute transfer through a liquid-liquid interface (LLI) by an activated process. This transfer has been investigated by molecular dynamics simulation (MD). Modified LennardJones potential functions have been used to model the solute and the two solvents which compose the LLI. The transfer reaction investigated is the change of solvation of the solute particle when it crosses the LLI. The full rate constant is calculated as the product of the transition-state theory (TST) rate constant and the transmission coefficient (dynamical contribution). The free energy barrier, encountered by the solute crossing the interface, is derived from the mean force acting on the solute. The TST rate constant of the transfer is deduced from the height and the shape of this barrier. The transmission coefficient is obtained from the plateau value of the normalized reactive flux. The prediction of the Kramers theory is compared with the MD results and discussed with respect to the specificity of our model.
I. Introduction The solvent extraction, or liquid-liquid extraction,I-* is a technical process which is frequently used to retrieve components from mixtures. The principle of this technique generally consists in the transfer of a solute from an aqueous solution to an organic solvent immisciblewith water. Thestudy of the extractionkinetics is rather complicated since the transfer involves many steps such as chemical reactions in the bulk or at the interface and transport phenomena. Among the numerous phenomena controlling the extraction rate, it is difficult (i) to determine which step controls the overall process and (ii) to estimate the rate constant of the interfacial reaction. Experimental studiess7 concerning solvent extraction have shown that the transfer rate can be limited by a chemical reaction occurring at the liquid-liquid interface (LLI) between the solvents. In this paper, we have chosen to model such a reaction occurring during the transfer of a solute particle. It is an elementary reaction consisting of a change of solvation of the solute. It is also important to recall that the phenomenological description generally used to clarify the experimental results is not able to provide reliable information on the mechanism occurring at the molecular scale. In this context, it is interesting to use a model approach which allows one to investigate directly at a microscopic scale the rate of the transfer reaction occurring at the LLI. Moreover, with our model it is possible to study only the interfacial chemical reaction without the transport contributions of the reactants and products. In order to achieve that, we have used the molecular dynamics (MD) simulation technique”” already employed in previous studies12-14 to investigatethe generic properties of planar LLIs. The interface model selected for the present study has already proved suitable to generate a stable interface.14 This system contains two atomic liquids A and B in which the interatomic forcesare described via modified LennardJones (LJ) potentials. In order to model the chemical reaction describing the change of solvation, we introduced, in the LLI model, an atomic solute interacting with a similar potential. For such a particle there is no steric effect and this choice simplifies the calculations. This model is simple but suitable to describe the change of solvation. Abstract published in Advance ACS Absrracrs, June 1, 1994.
0022-3654f 94/2098-6626$04.50/0
Indeed, the transfer reaction is investigated here from a generic point of view. Nevertheless, the change of solvation is a process occurring in every transfer reaction. The activation energy of the solvation reaction is obtained by implementing the model in such a way that the solute particle encounters an energy barrier in the interfacial region. The heterogeneous chemical reaction of the solute S going from a solvated state in liquid A to a solvated state in liquid B is given by
with k the rate constant of the transfer from liquid A to liquid B. The interaction potentials of the solute with the solvents have been chosen in order to obtain an energy barrier at the interface and a high solubility of the solute in both solvents. At this stage, we would like to emphasizethat no assumption is madeconcerning the shape of this barrier. This energy profile is obtained by calculating the potential of mean force acting on the solute as a function of the solute position with respect to the LLI. In this model, thechoiceof thereaction coordinatedescribingthe transfer is simple to define, since a good approximation is given by the solute position. An estimate, kmT, of the exact transfer rate constant, k,can be calculated from the transition-state (TST). One of the main approximationsof the TST is that when the activation state is reached, the system evolves necessarilyto give the products of the reaction. The r e c r o ~ s i n g s ~of~the J ~ activation ~~ barrier are thus neglected. These recrossings can be taken into account by analyzing the trajectories of the solute in the interfacial region with the reactiveflux method,22-26 which provides the transmission coefficient K . This coefficient corrects the TST prediction as follows:
There exist some theoretical approaches(e.g. Kramers theor~2~-29) which allow one to estimate the value of the transmission coefficient. MD simulations will provide a direct calculation of K and allow the comparison with the Kramers value calculated using the parameters yielded by the simulations. 0 1994 American Chemical Society
The Journal of Physical Chemistry, Vol. 98, No. 26, 1994 6627
Solute Transfer across a Liquid-Liquid Interface
TABLE 1: MD Simulation of the UP
I " " " " " ' , .
'
0
I
\
I
-*...-
I
I
u(r)
I
/
, \
\
0.6
-
10.3~ 35.0A
, I
u,
/
\
- 6 ' .
L, N PB T PNb 20.5~ 1728 0.845d 0.83alkg 1.35a/U3 70.0A 21.4 nm-3 99.0 K 0.56kbar a L,, L,, and Lz are the dimensions of the box. N,T,f", and pg are the particle number, the temperature, the normal component of the pressure tensor," and the bulk density of the two liquids, respectively. u = 3.405 A and c = 119.8kg;k g is the Boltzmann constant. Without long-range correction. L x = L,
I
f
us,
I
. . I . . . \ /
0.8
I
1
I . .
1
.
1.2
.
=Use I
I
a
1.4
I
. ' 1.6
r/o Figure 1. Interaction potentials. A, B, and S refer to particles of liquid A, liquid B, and solute S,respectively. See eqs 3 and 4. u = 3.405 A and c = 119.8kg;kg is the Boltzmann constant.
In section 11, we describe the manner of modeling the interactions of the solute particle with the solvents and we give the main features of the liquid-liquid interface. Section I11 is devoted to the methodology of calculation of the mean force potential and of the reactive flux. All the relevant details of the MD computations are also reported. The MD results for the transmission coefficient and the rate constant are given and discussed in section IV. The TST rate constant is deduced from the mean force potential, while the transmission coefficient is obtained from the plateau value in time of the normalized reactive flux. Next, we compare the Kramers prediction with the MD results. Then, our main conclusions are summarized in the last section.
11. Model A. The Liquid-Liquid Interface. The main features of the model interface we have used, in this study, are recalled hereafter, while ref 14 should be consulted for details. The model system contains two identical liquids A and B where the atoms interact via LJ potential functions U I J (IJ = AA, BB, AB) with the usual parameters u an c for argon (see relation 3). An extra parameter, CYIJ, has been introduced in order to obtain an interface between the two liquids. Although the two liquids are identical, they do not mix because the attractive part of UAB is reduced by a factor CYABequal to 0.3, whereas CXM and LYBB are simply equal to 1. These interactions are described in a condensed form by
uIJ(r)= 0 for r > rc
(3)
where rc = 3u is the cutoff radius. These potentials are displayed in Figure 1 . The system obtained with these potentials is symmetric with respect to the average plane of the interface. The simulation cell is a parallelepipedic box containing 1728 particles (864 particles of each liquid). The center of the box coincides with the origin of the coordinates; liquid A is on the side of negative z whereas liquid B is on the other side. The dimensions of this cell together with the thermodynamic conditions of the computationsare given in Table 1. Since our system is symmetric, the average plane of the interface, parallel to the xy plane, is located at z = 0. Periodic boundary conditions have been applied in order to (i) simulate an infinite size system in the x and y directions and (ii) to stabilize the average positions of the LLI. These conditions generate a second interface which is parallel to the first one and located at the boundaries of the simulation box. These two interfaces do not interact, because they are separated by regions which have the properties of bulk liquids A and B (same pressure at the same density and temperature).
The characterization of the interface properties has included the density and pressure profiles.14 The average position of the LLI is stable over the time scale of the simulations. The interface thickness is about three atomic diameters. One can emphasize an important point for the present study: the existence of a well in the density profile in the interfacial region. B. Modeling the Solute. In order to simulate reaction 1, we introduced in the previous model an atomic solute which has the same mass as the solvent particles. The selection of the potential parameters used to describe the interactions of this solute with the solvent is based on the following criteria. The solute must have a high solubilityin both liquid and must be strongly solvated. So, when it will move from one solvent to the other one, it will encounter an energy barrier corresponding to the replacement of the solvation shell composed of A by the solvation shell composed of B. In fact, the activated state should coincide with the bottom of the well of the total density14which is located at the position of the average plane of the LLI (z = 0). Thus, the solute will be less solvated in the interfacial region (see section 1I.C) and will have a higher energy. The potential function chosen to describe the interactions between a solute particle and a particle of liquid A or B is the following:
USA(')
= us&) = 0 for r > rc
(4)
where (YSA is equal to 2. This potential is more strongly attractive than U M and UAB, as shown by Figure 1. We do not need any solute-solute potential since we have introducedonly twosolute particles (seesection 111) in thesystem, one for each interface (see section 1I.A). In fact, these two particles are always separated by a distance greater than the cutoff radius of the potential. This corresponds to a low concentration of solute equal to 0.039 mol/L. At such a concentration, the interfacial tension is not significantly modified with respect to the value of 33 f 6 erg/cm2 obtained in the system without solute. Indeed, we have obtained 34 f 3 erg/cm2 with the two solute particles fixed at the average position of the interfaces by holonomic constraints (see section 1II.A). The properties of the solute particle are symmetric with respect to the average plane of the interface. The environment of the solute particle is the same in each liquid, and the distribution coefficient of the solute between the two liquids is equal to unity. This symmetry results in a simplification for the study of the transfer kinetics because (i) we have to deal with only half of the interfacial region to determine the energy barrier and (ii) we have to compute only one rate constant in reaction 1 since the forward and backward rates are the same. C. Solute-Solvent Radial Distribution Function. The computation of the transmission coefficient requires the analysis of the solute first neighbors shell. We need to know the radius, rs, and the number of neighbors, ns, of this shell. These values are obtained from the solute-solvent radial distribution function&). This distribution function has been computed in a bulk region of one of the two solvents and at the interface. The two distribution functions, plotted in Figure 2, are representativeof a small particle,
6628 The Journal of Physical Chemistry, Vol. 98, No. 26, 1994 7
"
"
?ss
6
#
~
~
~
"
"
~
~
interface
~
I -
5
bulk
3 2
---
0
' .
0.5
'
jb
'
' .
0.75
.
____--v
1
rs
t
1.25
.
1.5
'
'
'
'
'
'
1.75
'
'
2
rlo Figure 2. Solute-solvent (A and B) radial distribution functions, g(r), obtained at the interface and in a bulk region of one of the two solvents. rs is the radius of the first neighborsshell. The peak of the g(r) calculated at the interface is larger than in the bulk. However, the number of
neighbors is higher in the bulk (see text); this is due to the fact that the density is lower at the interface. u = 3.405 A. and the first neighbors peaks are very sharp, showing a strong solvation. The integration of these functions, up to rs = 1.1 u, gives an average value of the number of first neighbors, ns, in the bulk and at the interface. The number ns is smaller in the interfacial region (5.3) than in the bulkone (6.4). The difference ofsolvationof thesolutebetween the bulkand theinterfacecomes from the decrease of the total density of the system in this region. The solute is less solvated at the interface, and we expect an activated state correspondingto the change of solvationof reaction 1.
111. Computational Methodology All the MD simulations have been carried out at a temperature of 99 K, either in the microcanonical ensemble (NVE) or in the canonical ensemble (NVT) using the Verlet's algorithm3&)'with a time step 6t equal to s. In the (NVT) ensemble, the Noses scheme32 is applied with predicted velocitie~.~' MD simulations have been performed in order to calculate the rate constant of the transfer reaction given by relation 2. These calculations involve two main steps: (i) the computation of the free energy barrier via the calculation of the mean force potential (section 1II.A) and (ii) the determination of the transmission coefficient K using the reactive flux method (section 1II.B). As already stated in the introduction, there is no steric effect with our atomic solute, and we can consider that there is only a single degree of freedom involved in the transfer reaction. The corresponding reaction coordinate is thus simple to define by considering the position of the solute particle, ZS, with respect to the average position of the interface (z= 0; see section 1I.A). The reaction coordinate should depend on all the degrees of freedom; nevertheless, a very good approximation is provided by the Cartesian coordinate zs. A. Mean Force Potential. Our aim is to calculate the free energy profile associated with the transfer of the solute. In order to achieve that, it is sufficient to compute the difference, AF(zs), between the free energy, F(zs), of the system for a given value of zs and the free energy, P,of a reference state. P is chosen as the free energy of the system with the solute particle located in a bulk region of liquid A. The mean force potential, acting on a solute particle located at ZS, is the Helmholtz free energy difference, AF(zs):
is the zcomponent of the total force acting on the solute particle, and the brackets indicate a canonical average. ZQ is chosen in a region of the system for which the free energy is constant. This is verified by the fact that (fi(zs = ZO)) is close to zero. Relation 5 can be demonstrated by using the thermodynamic integration method.34
fi
TABLE 2
~
"
Computation of the Mean Force (f,)
~
"
~
"
"
~
purpose ensemble step 1 new configuration NVT NVE step 2 introduction of the solute NVT step 3 quilibration NVT step 4 computation offi 0 6 t = 10-14 S.
4
1
Hayoun et al. ~
~
t/6P
500 31000
2500 3000
~
~
constraints none along z alongz alongz
The choice of the method used to calculate the mean force acting on the solute takes into account the structural and dynamical properties of the interface. On average, the interface is located at z = 0. This mean position remains stable throughout the simulation: it is stabilized by the periodic boundary condit i o n ~ . 'The ~ instantaneous configurations of the LLI are not flat at the atomic scale (see Figure 8 of ref 14). The position of the interface determined locally on an area Ax Ay = u2 may slightly deviate from its average plane (z = 0). This deviation depends on the xy position and is autocorrelated over about 10 ps.14 Nevertheless, the mean position of each part of the LLI is zero, provided the average is computed over a sufficiently long period of time. On the other hand, the instantaneous position of the interface averaged in the x and y directions, on the whole area of the LLI, is also z = 0 at any time. Taking into account these properties, the position of the interface is z = 0 when we perform either a spatial or a time average. In these conditions the soluteinterface distance can be defined as Izsl. We can thus compute the mean force acting on the solute located at a given distance in two different ways. In the first case, the average is calculated on a long trajectory. In the second case, the averaging is performed on several short trajectories starting at different independent configurations. These configurations correspond to different sets of x,y values of the solute position in the plane z = 2s. These two sampling methods provide equivalent results, but the second one is more efficient. The calculation of the mean force at zs involves four steps defined by Table 2. A new equilibrium configuration of the interface is obtained at the end of the first step. Then, at the beginning of the second step, one particle of liquid A and one particle of liquid B are selected in a range Az equal to 0 . 0 5 ~ around a coordinate z* and replaced by two solute particles. This notation z* of the coordinate corresponds to the fact that we consider in the same simulation two solute particles which are at the same position with respect to the average position of the first and the second interfaces, respectively. The change of the nature of the two solvent particles is simply a replacement of the interaction potential: relations 3 are replaced by relations 4. Since the beginning of step 2, the coordinate zs of each solute particle is fixed at a given position z* by an holonomic constraint.35 In order to obtain the independent configurations mentioned above, we have repeated the whole simulation (four steps) nine times with different initial xy coordinates of the solute. B. Reactive Flux. The transmission coefficientcan be obtained via the calculation of the normalized reactive flux k(t).22-26 The direct calculation of the transmission coefficient K requires the generation of a large number of trajectories going through the barrier top. These spontaneous trajectories are rare events which are not attainable in a reasonable computational time. But the Bennett's method3638 is an alternative way to obtain them. In this technique, the solute particle is put in the close vicinity of the barrier top, at t = 0, and the evolution of the trajectory is followed forward and backward in time. The five steps, summarized in Table 3, describe how we proceed to obtain these uncorrelated trajectories. The solute particle is generated at zs = 0 f 0 . 0 5 ~and constrained to remain at this location. The average position of the barrier top is located at zs = 0, and its instantaneous equilibrium position fluctuates around this average at the same time as the LLI. So we also apply constraints along x and y in order to avoid the escape of the solute particle from the barrier
Solute Transfer across a Liquid-Liquid Interface
TABLE 3
The Journal of Physical Chemistry, Vol. 98, No. 26, 1994 6629
Computation of the Transmission Coefficienta purpose
t/6tb
step 1 new configuration step 2 introduction of the solute particles step 3 equilibration step 4 evolutionlt > 0 step 5 evolution/t < 0
constraints runs trajectories
lo00 none =1900 along x, y , z
108 108
2500 alongx,y,z SO00 none 5000 none
51 51 51
102 81 81
Note that the number of trajectories selected at step 4 is not equal to twice the number of runs because the trajectories which do not have a suitable environment (see text) have been rejected. 6t = l&I4 s.
top by moving in the xy plane. Then, after the equilibration step, the constraints are released and we verify that the solute is indeed in the close vicinity of the barrier top. Its environment at f = 0 must be roughly composed of half of solvent A and half of solvent B. If the environment does not fulfill this condition, the configuration will not be used to generate a full trajectory. This rejection of some configurationsprovides a simple way to obtain instantaneous initial configurations of the solute at the barrier top. The selection of configurations is not in strict agreement with the definition of the reaction coordinate,but this is not really a significant point. At the beginning of step 4 (defining t = 0), the constraints are removed. The z component of the velocity is given a random valuedrawn from a Gaussian distribution weighted by and referred to as Anderson’s one: i
~
~
~
Figure 3. Mean force acting on the solute particle (fz) plotted as a function of the solute position IS with respect to the average plane of the LLI. The dots refer to the computed data, whereas the continuous line is the result of a smoothing camed out over five points. u = 3.405 A and c = 119.8&~; &e is the Boltzmann constant. 5
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A
,
~
~
9
~
1
~
-
0 -3
~
~
-2
-1
0
1
2
3
0 The backward in time trajectory is obtained in step 5 with the same initial configuration (t = 0) as in step 4 by taking a negative value for the time step in the Verlet’s algorithm. These trajectories are then used to calculate the transmission coefficient of transfer reaction 1. The reactant state corresponds to negativevaluesofzs whereas, in the product state, zsis positive. The activated state, which is at the top of the free energy barrier, corresponds to zs = 0. The transmission coefficient is obtained from the reactiveflux, & ( f ) , which is a function of the reaction coordinate zs and of its velocity is:
k(r) =
(iS(0) 6(ZS(O)> e ( z S ( 0 ) ) (iS(0) 6(ZS(O)) W,tO)))
(7)
where (0) stands for t = 0 and 6 and 0 are the Dirac 6 function and the Heaviside step function, respectively. Because of the symmetry of the system, relation 7 is simplified (see expression A7 of the appendix) and the forward in time trajectories are sufficient to perform the computation. However, the backward trajectories allow good insight to be obtained on the whole process of transfer; this is why they have also been computed.
IV. Results and Discussion A. Free Energy Barrier and TST Rate Constant. The mean force acting on the solute particle has been calculatedas a function of the solute position, ZS. The results, computed for values of zs in the range [0,3.2u],are plotted in Figure 3. Notwithstanding the number of independent configurations taken into account to calculate the forces, there is a significant scatter of the computed values. The mean force obtained for zs = 0 is very close to zero (O.OSc/a), and indeed the expected equilibrium value must be zero since the interface is symmetric. This result justifies a posteriori the sampling method used to calculate the mean force (section 1II.A) since it shows that ksl defines satisfactorily the distance to the LLI. We also note that the mean force is nearly zero for a distance to the average position of the LLI greater than
Figure 4. Freeenergy, AF,and potential energy, AU,barriers calculated with a reference state chosen as 9J= - 2 . 8 ~( u = 3.405 A). The left parts of the barriers have been calculated, whereas the right parts are obtained by symmetry with respect to the plane zs = 0.
2 . 8 ~ .The mean force potential is thus constant beyond this distance. So, the free energy profile can becalculated, according to relation 5 , taking as a reference state the value I” defined by the solute located at zs ,< zO = 2 . 8 ~ Only . half of the values of the mean force potential have been computed (negative 2s). The whole profile plotted in Figure 4 takes into account the symmetry of the system. The smooth aspect of this profile is due to the integration (relation 5 ) of the mean force. Thisfree energy profile has the shape of an energy barrier. The model has thus produced the activationbarrier required by reaction 1. This result confirms the pertinence of the choice of the potential parameters used to describe the solute interactions. The barrier height hFy = 4.2keT is sufficiently high to model an activated process. It gives the activation free energies of the forward and backward reactions. The TST estimate for the rate constant, k, of reaction 1 is given byz5
where the brackets denote a canonical average. A simple calculation leads to
Relation 9 shows that kmT depends not only on the height of the activation barrier but also on the shape of this barrier through the integral of the denominator. The frequency at which the system leaves the activated state is given by the first term of
6630 The Journal of Physical Chemistry, Vol. 98, No. 26, 19I94
Hayoun et al.
1
0.8
ht)
1
0.6 0.4
-1
-2 0
10
20 30 Time (ps)
40
50
Figwe 5. Normalized reactive flux plotted as a function of the time. The
transmission coefficient corresponds to the plateau value reached after 10 ps.
expression 9, calculated by numerical integration, and is equal to 46 ns-1. The value of the TST rate constant is 0.71 ns-l. Since the computation of this barrier requires a free energy calculation, it is always time consuming. It is interesting to check if the potential energy barrier could approximate the free energy barrier. To this purpose, we have calculated the potential energy, U(zs),of the system as a function of the solute position. U(ZS) remains constant when the distance to the LLI is greater than 2 . 8 ~ . This defines the reference state used to calculate the potential energy difference AU. This quantity is plotted in Figure 4 together with AF. The comparison of the two profiles shows that the potential energy barrier is rather close to the free energy one. This potential energy profile yields a barrierbeight AV equal to 3 . 9 k ~ Tand a corresponding TST rate constant equal to 0.90 ns-1. With our model, the potential energy gives a good approximation of the TST rate constant. This result is clearly due to the fact that the free energy and potential energy profiles are close. It means of course that the entropy variation due to the crossing of the interface by the solute is not very important. The potential energy profile could be used in other systems to calculate the TST rate constant. This approximation is only valid if it is possible to prove otherwise that the entropy contribution is negligible. B. Transmission Coefficient. The reactive flux has been calculated with 81 full trajectories (forward and backward in time) for which the initial environment of the solute particle was composed, on average, of 2.7f 0.7 particles of solvent A and 2.8 0.7 particles of solvent B. The average total number of solvent particles (Aand B) is equal to 5.5 0.7and agrees very well with the value of the number of first neighbors, ns * 5.3, obtained from the radial distribution function at the interface. If a macroscopic rate law exists, the transmission coefficient, K , will correspond to the plateau value in the time variation of the normalized reactive flux. The normalized reactive flux has been computed according to relation A7. The result obtained is given in Figure 5. The reactive flux decreases from the TST value, k(t = O+), equal to 1 at t = 0+, and reaches a true plateau. The existence of this true plateau proves26 the validity of the macroscopic rate law corresponding to reaction 1.
*
*
where C ~ and A C ~ are B the solute concentrations in solvents A and B,respectively. Since the system is symmetric, the forward and backward rate constants are identical and are both referred to as k. The time needed to reach the plateau (10 ps) is an intermediate one. It is smaller than the average time between two reactions but greater than the time required for the equilibration of the solvent with the solute.
-3 Time (ps)
Figure!6. Solutetrajectories, generated by Bennett's method,illustrating a solute transfer (dotted line) and an unsuccessful attempt (continuous line). The diffusivebehavior is clearly visible, particularly at the barrier top ( i s = 0). It is characterized by a large lifetime of the activated state associated with multiple recrossings.
The transmission coefficient, K , represents the dynamical contribution to the full rate constant and accounts for the recrossingsof the barrier top induced by the solvent. It has been obtained from the plateau value of the normalized reactive flux which has been estimated by averaging between 10 and 50 ps. We obtained K = 0.15 0.05. This value, significantly smaller than 1, indicates a large deviation from the transition-state theory, which assumes that no recrossing occurs ( K = 1). This result can be qualitatively understood by looking at the solute trajectories crossing the interface (Figure 6). These trajectories exhibit a diffusive behavior even in the vicinity of the barrier top, where many recrossingsof the activated state occur. In theseconditions, the correlation with the initial positive velocity is lost and consequently K is small. Another consequence of this diffusive behavior associated with the multiple recrossings is the large lifetime (about 20 ps on average) of the activated state. C. Full Rate Constant. The full rate constant, k, obtained as the product of the TST rate constant (0.71 ns-1) with the transmission coefficient (0.15)is equal to 0.11 ns-1. The kinetics of the solute-transfer reaction is heterogeneous, and therefore therateconstant kdependson theratioAlVbetween the interface area and the volume of the system. So, in order to obtain a rate constant independent of the system size, we define k'by
*
k' = k-V A The value of the full rate constant, corresponding to the transfer of the solute through the energy barrier located at the interface, is k' = 37 cm s-l. This rate constant is large with respect to the experimental values. The reaction rate is fast, but we have only investigated a chemical reaction at the interface which is not always the limiting step of the whole experimental process. There is no easy comparison between the value of this rate constant computed for an elementary step of the liquid-liquid extraction, namely the change of solvation of the solute at the interface, and the experimental results. If the chemical reaction at the interface is the limiting step in the experimental process, the activation energy can be determined. The experimental energies, thus obtained, are significantly higher than our barrier height. This difference is explained by the following remarks. The interfacial reaction may involve not only a solvation change but also other phenomena. Even if the interfacial reaction is simply a solvation change, the experimental energy barrier will be higher since in our model the solute does not have a rigid solvation structure. D. Application of Kramers Theory. We can now discuss the possibilityof calculating an approximate valueof the transmission coefficient. Such an approximation can be obtained with the Kramers theory*'-29 (KT), which is a stochastic approach. In
Solute Transfer across a Liquid-Liquid Interface
The Journal of Physical Chemistry, Vol. 98, No. 26, 1994 6631
this theory, thesolventis takenintoaccount by meansofaconstant friction coefficient, 5: This coefficient can be estimated from the diffusion coefficient, D,of the solute in the solvent, via Einstein's relation. The low solute concentration (0.039 mol/L; two solute particles in the system) makes it difficult to obtain an accurate determination of D from the usual equilibrium methods: mean square displacement or velocity autocorrelation function. So we have used a nonequilibrium molecular dynamics4M2J0 method, which is more efficient in such conditions. The statistics of the result are improved by using the subtraction technique.4345The diffusion coefficient thus calculated is introduced in Einstein's relation in order to estimate the friction coefficient. We obtained 127 f 4 reduced units or 5.7 f 0.2 X 10-10 g/s for S: The Kramers transmission coefficient, KKT, is given by19928
with ubthecurvatureofthebarrier top. Thiscurvatureisobtained by fitting the equation of a parabola (13)
on the top region of the free energy barrier (Figure 4). Ob has been found equal to 1.0 ps-1. Considering the numerical values of the barrier curvature and the friction coefficient, it is worth pointing out that the diffusive limit of Smoluchowski~J7is reached, since S" >> 1. In this 2wb case the value of the transmission coefficient is given by (14)
The transmission coefficient, calculated in this way, is equal to 0.12. This prediction is close to the MD result. This statement is not so surprising since, as already stated, the behavior of the solute during the whole transfer process is diffusive. In order to obtain quantitative arguments explaining why Kramers theory gives a good estimate for K, even though it contains the simplifying assumption of the zero frequency friction coefficient, we have investigated relevant characteristic times. The appropriate comparison is between the time evolution of the reactive flux and the correlation time of the friction forces exerted on the solute. In order to investigate the friction characteristics, we have calculated the autocorrelation function of the force exerted on the solute by the solvent in the activated state. Although only the z component of the force is required for our purpose, we have also studied the x and y components to have a reference nonperturbedby the local equilibriumfluctuationsof the interface. For the z component, the analysis of the correlation function is not so obvious as for the x and y components. The same decorrelation time is observed for the three directions, but the long time tail does not fluctuate exactly around zero for the z component. Nevertheless, it oscillates between positive and negative values, and the departure from zero is small (less than 1%). This behavior could be explained by the fact that although the solute is constrained to remain at z = 0, we are not sure it stays at the barrier top (local fluctuations of the LLI). After some picoseconds, the constrained position of the solute could perhaps induce a local deformation of the interface. However, it is clear that the decorrelation occurs in 2 ps; the characteristic time of the friction forces is thus significantly smaller than the time required for the reactive flux to reach its plateau value (10 ps). In these conditions, the Kramers theory is applicable and we do not observe any inconsistency as reported by Ciccotti et a1.47 It is not necessary to perform a complete analysis to decide whether or not the KT approximation in the Smoluchowski limit
is valid to calculate the transmission coefficient. A qualitative criterion can be suggested. It is simply necessary to calculate some full trajectories of the solute crossing the interface. If these trajectories exhibit a diffusive behavior during the crossing of the barrier, the Smoluchowski approach is suitable. V. Concluding Remarks
In this paper, we have investigated the kinetics of the transfer of a solute across a liquid-liquid interface. The transfer is a thermally activated process consisting of the change of solvation of the solute occurring at the interface. This elementary chemical reaction is of course a very simple case of transfer, compared to experimental ones. But nevertheless, the change of solvation occurs in any solute transfer and can be considered as a generic feature of the transfer. Since our approach is performed at the atomic scale, the separation of the contribution of the chemical reaction from the transport contributions (to and from the LLI) is easily done and the kinetics of the reaction can be quantitatively studied. We have shown, for a simple model, that it is possible to calculate, by MD simulation, the full rate constant of the solute transfer. This study has also allowed us to compare for the same model the activation energy and the transmission coefficienttoapproximatevalues of these quantities. Some recipes explaining when it is possible to use these approximations are also given. The full rate constant has been determined from both the activation and the dynamical contributions. This prediction is made without any approximation for the shape of the activation barrier. The free energy barrier encountered by the solute is a mean force potential calculated by averaging over the configurations of the solvent. It provides the activation contribution to the rate constant and leads to the TST estimate. The activation free energy is compared to the approximated value obtained from the potential energy profile. The analysis of the reactive transfer trajectories has shown that the dynamics of the solvent induces recrossings at the barrier. Therefore, the TST prediction must be corrected with the transmission coefficient. This coefficient has been obtained from the plateau value of the normalized reactive flux and compared to the Kramers theory estimate. The two values are in good agreement, and this is explained by the diffusive behavior of the solute. In the present state of the model, there is no direct comparison with the experimental results. In fact, the aim of this study was the computer modeling of the transfer by an activated process. The advantage of the MD simulation is the ability to investigate separately the different steps of the whole process. The counterpart of such an approach is the necessity to begin with a simple model. It is possible to foresee developmentswith a more complex model involving, for instance, a more rigid solvation. A possible consequence of such a solvation could be a higher energy barrier. A more realisticdescription of the liquid-liquid extraction could be obtained with polar solvents and a molecular solute.
Acknowledgment. The authors thank D. Borgis, G. Ciccotti, and V.Pontikis for helpful discussions. The scientific direction of SPM in CNRS is also acknowledged for the allocation of computing time. Appendix. Calculation of the Reactive Flux In our system, the coordinate of the activated state is located at zs = 0. The B(zs(t)) term of relation 7 characterizes the reaction: its value is 1 when the solute is located between the activated and the product states, whereas it is equal to zero when the solute is on the reactant side. In order to compute k(t), relation 7 is expressed as a function of B(z~(0)) and B(-ks(O)) by introducing O(iS(0)) + O(-&(O)) = 1 into the numerator. Since we are dealing with a system at equilibrium, we also use the following relation:
6632 The Journal of Physical Chemistry, Vol. 98, No. 26, 1994
(iS(0)6(ZS(O))) = 0
(AI)
to modify the denominator. We obtainz5
Hayoun et al. sign of the time, the velocity remains positive and points toward the product state. These conditions are the right ones required to computeexpression A7. Finally. the average has been achieved over 162 halitrajectories (forward or bac&ard ones).
References and Notes (1) Handbook of Solvent Extraction; Lo, T. C., Baird, M. H. I., Hanson, C., Eds.; Wiley: New York, 1983. (2) Science and Technology of Tributyl Phosphate; Schultz, W. W., Navratil, J. D., Talbot, A. E., Eds.; CRC Press: Boca Raton, FL, 1984; Vol. I. (3) Lewis. J. B. Chem. Enp. Sci. 1954. 3. 248. (4j Horner, D.; Mailen, J.;~iel,S.;Scdtt,’T.;Yates,R.Ind.Eng. Chem.
which can be symbolically written as
k(t) = V ( z s W ) ) +
(A3)
Let us now describethe details of the computation of k(t). We could perform the calculationdirectly since we have forward and backward trajectories, but we will take advantageof the symmetry of the system to simplify relation A3. In what follows, we will only discuss in detail the case of the forward in time trajectories. It is clear that (i) because of the 6 functions, the trajectories for which the solute position zs is not strictly equal to zero at t = 0 do not contribute to the averages. So, we should consider only forward trajectories with zs = 0 at the time origin. The reaction coordinate we have chosen at t = 0 is close to zero (f0.05~) but not exactlyequal to. So we made a small approximation: instead of computing O(zs(r)), we computed B(zs(t)- zs(0)) for t close to zero. (ii) Rather than sampling the time origin velocities from a Maxwell distribution (density of probability) and weighting the averages by is, it is more efficient (we need less trajectories to obtain the same statistics) to compute simple averages after thevelocities have been sampled from the Anderson’s distribution (relation 6). The time origin velocities of the solute particle used in step 4 have been randomly drawn from this distribution. If we consider the first term of relation A2, we see that, because of the B functions of the velocity, there is no contribution of the forward in time trajectories with a negative velocity at t = 0. This, added to remarks i and ii, leads to (e(zS(0) )+
= (e(zS(0))
I
r,(O)=O
(A41
iS’0 Anderson
The second term of expression A2 requires the contribution of negative velocities: (e(zs(t)) )-
= (e(zs(r)))
I
(A5)
; : % = O n
but since our system is symmetric, we may exchange the role of the reactant and product states. If we deal with the same forward trajectories with such an exchange,the velocity will point toward the reactant. That corresponds to a negative velocity for the reaction coordinate. Notwithstanding the exchangebetween the reactant and the product states, the same forward in time trajectories are used; this leads to (~(ZS(0))-
= 1 - (e(zs(t))>+
(-46)
and expression A3 becomes
= 2ms(t)))+ - 1
(A71
Forward in time trajectories are suitable to calculate relation A7. Nevertheless, backward trajectories can also be used. The negative time involved is dealt with as if it were positive, provided that the reactant and product states are inverted. Since the direction of the reaction coordinate is reversed together with the
--.
Fundam. 1980. -19., 103. .--.
-(e(zs(t))>-
(5) Murdoch, R.; Pratt, H. R. C. Trans. Inst. Chem. Eng. 1953,31,307. (6) Lewis, J. B. Nature 1956, 178, 274. (7) Lewis, J. B. Chem. Eng. Sci. 1958,8, 295. (8) MolecularDynamics Simulation of Statistical-Mechanical Systems; Ciccotti. G., Hoover, W. G., Eds.; North-Holland Amsterdam, 1986. (9) Simulationof LiquidsandSolids-Molecular Dynamicsand Monte CarloMethods inStatisticalMechanics; Ciccotti,G., Frenkel, D., McDonald, I. R., Ed%;North-Holland Amsterdam, 1987. (IO) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (1 1) ComputerSimulation in Materials Science;Meyer, M., Pontikis, V., Eds.; NATO AS1 Series: Kluwer Academic Publishers: Dordrecht, The Netherlands, 1991. (12) Hayoun, M.; Meyer, M.; Mareachal, M.; Ciccotti, G.; Turq, P. In Chemical Reactivity in Liquids, Fundamental Aspects; Moreau, M., Turq, P., Eds.; Plenum: New York, 1988; p 265. (13) Hayoun, M.; Meyer, M.; Turq, P. Chem. Phys. Lett. 1988,147,203. (14) Meyer, M.; Mareschal, M.; Hayoun, M. J. Chem. Phys. 1988.89, 1067. (15) Eyring, H. J. Chem. Phys. 1935,3, 107. (16) Evans, M. 0.; Polanyi, M. Trans. Faraday Soc. 1935,31,875. (17) Glasstone,S.;Laidler,K. J.; Eyring, - - H. The TheorvofRateProcesses: McGraw-Hill: New York, 1941. (18) Truhlar, D. G.; Hase, W. L.; Hynes. J. T. J. Chem. Phys. 1983,87, 2664. (19) Hynes, J. T. In Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC Press: Boca Raton, FL, 1985; Vol. IV, p 171. Wigner, E. J. Chem. Phys. 1939, 7, 616. (20) Hirschfelder, J. 0.; (21) Anderson, J. B. J. Chem. Phys. 1973, 58,4684. (22) Chandler, D. J. Chem. Phys. 1978,68,2959. (23) Rosenberg, R. 0.;Beme, B. J.;Chandler, D. Chem.Phys. Lett. 1980, 75, 162. (24) Berne, B. J. In Multiple Time Scales; Brackbill, J. U.,Cohen, B. I., Eds.; Academic Press: New York, 1985; p 419. (25) Berne, B. J.; Borkovec, M.; Straub, J. E. J. Phys. Chem. 1988,92, 3711. (26) Chandler, D.; Kuharski, R. A. Faraday Discuss. Chem. Soc. 1988, 85, 329. (27) Kramers, H. A. Physica 1940, 7,284. (28) Chandrasekhar, S. Reo. Mod.Phys. 1943, 15, 1. (29) Hinggi, P.; Talkner, P.; Borkovec, M. Reu. Mod. Phys. 1990,62, 251. (30) Verlet, L. Phys. Reo. 1967, 159, 98. (31) Berendsen, H. J. C.; van Gunsteren, W. F. In ref 8, p 43. (32) Nost, S. J. Chem. Phys. 1984,81, 511. (33) Ferrario, M.; Ryckaert, J. P. Mol. Phys. 1985,54, 587. (34) Frenkel, D. In ref 8, p 151. (35) Ryckaert, J. P. In ref 8, p 329. (36) Bennett, C. H. In Diffusion in Solids; Nowick, A. S., Burton, J. J., Eds.; Academic Press: New York, 1975; p 73. (37) Bennett, C. H. In Algorithms for Chemical Computation; Christofferson, R. E., Ed.; ACS Symposium Series 4; American Chemical Society: Washington, DC, 1977; p 63. (38) Gillan, M. J.; Harding, J. H.; Tarento, R. J. J. Phys. C 1987, 20, 2331. (39) WecanobtainvaluesofXdistributedaccordingtoP(X) BXedX’I2 by using the following relation: X = (-(2/B)ln(l - Y))V*, with Y values randomly drawn from a uniform distribution. (40) Hoover, W. G. Physica A 1983,118, 11 1. (41) Hoover, W. G. Annu. Reo. Phys. Chem. 1983, 34, 103. (42) Evans, D. J.; Hoover, W. G.; Failor, B. H.; Moran, B.; Ladd, A. J. C. Phys. Reu. A 1983, 28, 1016. (43) Ciccotti, G.; Jacucci, G. Phys. Rev. Lert. 1975,35, 789. (44) Ciccotti, G.; Jacucci, G.; McDonald, 1. R. Phys. Rev. A 1976, 13, 426. (45) Ciccotti, G.; Jacucci, G.; McDonald, I. R. J. Srat. Phys. 1979, 21,
__
-
1.
(46) Smoluchowski, M. Z . Chem. Phys. 1918, 92, 129. (47) Ciccotti, G.; Ferrario, M.; Hynes, J. T.; Kapral, R. J. Phys. C h m . 1990, 93,1137.
