J. Phys. Chem. 1995,99, 4293-4305
4293
Nonequilibrium Solvation: The Mutual Influence of Solute and Solvent Dynamics Manuel A. Aguilar and Antonio Hidalgo* Chemical Physics Department, University of Extremadura, Badajoz, Spain Received: September 29, 1994; In Final Form: December 29, I994@
A new approach for the study of the mutual influence of solute and solvent dynamics is presented. This method allows one to specify the solute-solvent Hamiltonian using the in vacuo solute potential and the nonequilibrium solvation free energy. A criterion to discard or to take into account the different solvent coordinates, depending on their influence on the solute dynamics, is established. Variational results are obtained using the global (two-dimensional) Hamiltonian and the two one-dimensional limits, SCF and BO. Two typical dynamic processes have been studied, the vibrational motion of a diatomic molecule and a proton transfer reaction in solution.
I. Introduction There exist a great number of phenomena occurring in solution for which a static description of the solvent is insufficient. Electronic transitions, vibrational spectra, or chemical reactions in solution are only a few examples that can be mentioned. In all these cases it is necessary to describe the solvent in nonequilibrium situations. Because of the difficulty of nonequilibrium problems, the study of the effects of solvent dynamics has not been carried out at the same level of sophistication as equilibrium solvation and very simple models must be employed. For nonequilibrium phenomena, two theoretical treatments have been used in recent years: continuum models' (CM) and the mean spherical approximation (MSA)? The CM focus on the long-range solute-solvent nonspecific interactions and predict that solvation should proceed on a time scale that is approximately equal to the longitudinal relaxation time of the solvent. The MSA goes beyond a continuum representation by recognizing the molecular nature of the solvent and predicts much more complex dynamics. Although the MSA is a more sophisticated theory, in most cases, it is found that the average solvation times measured in experiments are close to the continuum prediction. Taking this into account, Maroncelli et aL3 have proposed a continuum model that permits one to calculate the solvent frequency. Their results showed close agreement with the predictions of computer simulations. Prior to the study of solvent dynamic effects, the number and nature of the solvent coordinates must be determined. The solvent coordinates are independent scalar variables that describe the degree of deviation of the solvent polarization from equilibrium. The problem is then to select from all possible solvent polarizations those which are strongly coupled to the solute dynamics. Several approximations have been proposed in the literature. In electron transfer reactions, for instance, the solvent polarization is constrained to be intermediate between the equilibrium polarization field for the donor and acceptor diabatic charge-localized states. 'b,4 An extension of this approximation to the case of chemical reactions has been proposed by Hynes and c o - w o r k e r ~ .In ~ ~a ~ proton transfer reaction, the diabatic states either are taken to be the neutral and ionic forms of valence bond states for a hydrogen bonded complex5 or are replaced by hypothetical charge distributions in which the charges appropriate to the
* To whom the correspondence @
should be addressed. Abstract published in Advance ACS Abstracts, February 15, 1995.
0022-3654/95/2099-4293$09.0QIO
reactant or product state are held at the geometries corresponding to the different values of the reactive coordinate ( S Nand ~ S N ~ reactions) .6,7 For other authors,8s9the set of solvent polarizations considered consists of the set of equilibrium solvent polarizations for other solute geometries. In other words, the nonequilibrium solvent polarization experienced by the solute at the geometry x is modeled by the equilibrium solvent polarization when the solute is at another geometry x'. The above definitions present some problems. The election of the diabatic states is complicated, and the associated solvent coordinate is only valid if the solute wave function may be written as a linear combination of the two diabatic states. If more complex wave functions are used (CI, for instance), a larger set of solvent coordinates must be introduced. In this case it is necessary to consider as many solvent coordinates as there are electronic configurations.I0 The election of the equilibrium solvent polarization to obtain the nonequilibrium values is a convenient way to represent the effect of the solvent delay, but clearly, other solvent polarizations (associated with thermal fluctuations, for instance) that affect the solute dynamics are not taken into account. In section 11, we present a method of selecting solvent coordinates that, without making any reference to diabatic states, allows one to include all the solvent polarizations relevant to the process under study. In our approximation, the change of the solute charge distribution during the course of the process determines the solvent coordinates that are to be included in the dynamics. More specifically, the solvent coordinates will be determined by the solute multipoles. A very similar approximation for the case of a dipolar solute with fixed spatial orientation has been proposed by Berezhk~vskii.~~ This author, however, neglects completely the effect of the self-polarization of the surface charges. In our model this effect is taken into account. Section I1 also introduces explicit expressions for the nonequilibrium free energy. The Hamiltonian we use is discussed in section III. It is a special case of the type of Hamiltonian extensively employed in the study of chemical reactions in s ~ l u t i o n . ~The - ~ Hamiltonian couples the solvent coordinate above indicated with a solute coordinate. Its nature depends on the process under consideration. We have also studied two limit situations usually treated in the bibliography: the Born-Oppenheimer (BO) and the self-consistent-field (SCF) limits. Both give place to onedimensional Hamiltonians. The BO limit corresponds to a slow 0 1995 American Chemical Society
4294 J. Phys. Chem., Vol. 99, No. 12, 1995
Aguilar and Hidalgo
solvent response regime, and SCF, to a fast solvent response regime (i.e., equilibrium solvation). Section IV describes the numerical method used to solve quantum-mechanically the two-dimensional Hamiltonian and the two limits BO and SCF. Applications to molecular vibrations and proton transfer reactions are described in section V. Finally, section VI presents our conclusions.
AG
h
11. Nonequilibrium Free Energy and Solvent Coordinates We first present the expressions for the electrostatic free energy for a solute with dipole moment ,&ac immersed in a medium with orientational polarization in equilibrium with an arbitrary dipole moment ji,. Two treatments, classical and quantum, are possible: (i) In the classical case, the nonequilibrium electrostatic free energy for a polarizable dipole is (see the Appendix)
where p ( c ) is
and
K = g(Eo) - g(Ew)
(3)
The form of g(E), the factor of the reaction field, depends on the model and cavity chosen. In the present work, CM and a spherical cavity are used. The extension of the formalism to MSAI2 is immediate. The use of elliptical or general cavities is also possible, but then the expression for the free energy and the solvent coordinate definition are somewhat more complicated due to the tensorial nature of g(E). If the cavity is a sphere, g(c) is the Onsager reaction factor:I3
2 ( E - 1) 1 g(E) = 2 E 1 a3
+
(4)
with a the cavity radius. The free energy for a nonpolarizable dipole moment can be obtained as a particular case of eq 1 by setting a = 0.
AG = - 1 / 2 g ( ~ o ) P v a c 2
+ '/&Gs
- PvaJ2
(5)
(ii) In the quantum treatment, the electronic wave functions of the solute are obtained via the effective Schrodinger equation
(2+ P2)IY) = WlY)
(6)
A3
eq
I
I
I "BE
Ps
Figure 1. Free energy versus different solvent configurations characterized by a dipole moment &). The equilibrium situation corresponds to the minimum when ,us= ,uvac.
Throughout this paper we take ji, as a dynamic variable associated with the solvent. It represents the actual orientational polarization of the solvent. Due to the vectorial character of ji,, three scalar coordinates (the three Cartesian components, psr,,us,,p,, for instance) will be necessary to represent completely the degrees of freedom of the solvent. Evidently the sofvent has many more degrees of freedom, but they are not relevant to the dynamics of the solute system because they are not coupled to the solute coordinate. We will retum to this point below. Tuming to eq 1, this only considers the electrostatic component. It would be interesting to include the dispersion component-the energy due to the correlation between the electronic density fluctuation in the solute and in the solvent. This component may be calculated by making use of an expression due to Linder:I4
where WI and (02 are the average electronic transition frequencies of the solute and solvent, respectively. Strictly, the dispersion component should also be included in the electronic effective Hamiltonian (eq 6). However, it has been ~ h o w n ' ~that * ' ~this component has a very small effect on the solute wave function. By contrast, it can have a major effect on the free energy and, thus, be included as an additional term in the free energy (eq 1). As we have above indicated, three solvent coordinates are needed to characterize completely the orientational polarization of the solvent in the most general case. However, in many interesting cases, a smaller number of coordinates may be used. To see this, we shall write eq 1 in a more convenient form:
where is the electric field generated by the solvent in the position occupied by the solute, in our case, the cavity center. W is the electronic Hamiltonian of the solute system in vacuo. The free energy of the system is now
+
AG = -1/2g(Eo)p2 1 / 2 K G s F)'
+ (Yl$IY)
- E"
(7)
withp = (Y1,fiIY)and IY) the solute wave function in solution.
Eo is the in vacuo solute energy. In eqs 1,5, and 7, the first term is the equilibrium free energy of an ideal dipole and the second is the distortion or reorganization energy, ;1 in Figure 1 (the free energy difference between an arbitrary solvent configuration and the thermally equilibrated solvent configuration). This last term describes the degree of deviation of the actual orientational polarization, characterized by p,, from the equilibrium one.
The number of necessary solvent coordinates is then given by the number of coupling terms that appear in eq 9. Three cases can be considered: (1) &ac = pxT. The solute dipole varies in magnitude but not in direction. Only the pSxcomponent of ,E, is coupled to the solute coordinate; the psy and psz components do not influence the solute movement. Equation 1 now becomes
J. Phys. Chem., Vol. 99, No. 12, 1995 4295
Mutual Influence of Solute and Solvent Dynamics ~
y P;
+ 2m2 +
H = - - Q(r) 2m,
with s = p, and where only scalar functions have been used. This case is found in diatomic vibrations and several chemical reactions that include electron transfer, proton transfer, or sN2 reactions and is the one we shall study in the rest of the paper. The extension of eq 10 to ellipsoidal cavities is immediate if the solute dipole varies along one of the major axes of the ellipsoidal cavity. Otherwise, a greater number of coordinates will be necessary. (2) ,Ev,, = p,Z p;. The solute dipole varies in direction and magnitude. This is the case found in most electronic transitions. If we take the dipole moment of the excited state to lie along the x direction, then the ground dipole moment will lie in the plane xy. Two solvent coordinates along the directions x and y are necessary to define the nonequilibrium solvation state and further dynamical solvent relaxation toward equilibrium of the solute excited state. (3) ,Ev,, = p,Z p; p&. This is the general case, the three components of p varying over the course of the process. Three solvent coordinates are necessary to characterize the solvent polarization. As one can see, in the above definition of the solvent coordinates, no reference has been made to the form of the solute wave function and very complicated wave functions could be used. The solvent coordinates are determined only by the number of solute dipole components that change appreciably during the process. This simplifies the problem with respect to the usual formulation where the number of solvent coordinates is related to the number of electronic configurations of the solute wave function.'O Finally, we must remark that it may be necessary to introduce new coordinates when the electrostatic free energy is a function of higher order multipoles, Qn.In this case, the free energy is"
+
+
+
where
In general, it will be necessary to introduce one solvent coordinate for every multipole that varies substantially over the course of the process. If a multipole is constant during the solute dynamics, the corresponding coupling term, derived from the expansion of eq 11, enters the free energy as a function only of the solvent coordinate and hence it does not modify the solute dynamics. 111. Hamiltonian
In what follows, we describe the solute motion with the aid of a single coordinate. Its nature depends on the process under consideration. It may describe a motion along the reaction coordinate or may be a vibrational coordinate. Extensions to systems with a greater number of coordinates will be the subject of a subsequent paper. The nuclear Hamiltonian of the system (solute solvent), expressed in internal Jacobi coordinates (after neglecting the angular dependence), is
+
+ AG(r,s)
(13)
where r is the internuclear separation, e(r) the potential energy function of the dynamic process (reaction or vibration) in vucuo, AG(r,s) the solvation free energy given in eq 10 but where the solute dipole moment function is expressed in units of distance, and g(c) is now
2(€ - 1 ) e2 2E 1 u3
g(€) = -
+
e being the electron charge to simply convert the reaction factor g(6) to energy/distance2 units.
A central role in what follows will be played by the solvent force constant Ks0l: d 2 ~ ~ ( ~ , s )
Ksol =
dS2
Its value is
When a = 0, the value of Ks0l is constant (Kso1= K), but in the general case, Ksol is a function of a, the solute polarizability, which usually depends on the solute coordinate, a a@). The precise expression that Ks0l takes depends on the definition of the solvent coordinate. In our case, and for convenience, the solvent coordinate is defined in such way that it has units of distance. Other definitions are possible; for instance, Fonseca et al." express the solvent coordinate in units of energy and find
where the influence of electronic polarization has been neglected. The precise form of Ksol does not affect the dynamics of the system because the mass associated with this motion also depends on Ksol. In physical terms, what matters is not the solvent force constant nor the solvent mass separately but the frequency of the solvent motion, and this is independent of the solvent coordinate definition. In eq 13, ml is the mass associated with the solute coordinate and m2 is the reduced mass of the solvent and the dynamic system:
where mT denotes the total mass of the solute and m, the solvent polarization mass. We are not able to assign a numerical value to m,,but we can calculate m2 directly from the relationship
where osolis the solvent frequency. This frequency can be obtained from computer simulations or calculated theoretically. Several theoretical methods have been proposed in the literature with the aim of obtaining wsol. In the present paper, we employ an extension to the dipole case of the method proposed by Maroncelli3 et al. The final result is
Aguilar and Hidalgo
4296 J. Phys. Chem., Vol. 99, No. 12, 1995
and hence
where psolis the dipole of the solvent molecules, e the solvent number density, I the moment of inertia of the solvent molecule, and (w2)the average squared rotational frequency of the solvent, determined entirely by the inertial properties of the solvent molecules. In what!ollows, eq 10 will be used together with the solute potential V(r) to construct potential surfaces with explicit inclusion of the solvent coordinate and the nuclear Hamiltonian (eq 13) will be used to study the dynamics of the solute-solvent system. As two particular cases of the complete Hamiltonian (eq 13) we have also studied the BO and SCF approximations. Both limit cases permit the transformation of the two-dimensional Hamiltonian into a one-dimensional Hamiltonian; they must however be used with caution because, as we will show below, they are only valid under restricted conditions. In the SCF limit, the solvent motion is much faster than the solute intemal movement. The solute atoms are virtually fixed in space while the solvent orientational polarization bends. For each value of the solute coordinate R, the solvent “sees” the solute as if it were frozen in time, and the solvent is always in the equilibrium position for the current displacement of the solute. This limit may be obtained by setting
and is valid when osol>> WR. In this case the solute Hamiltonian becomes
In the BO limit, the motion of the solute atoms is now faster than the solvent motion. The solvent dependence enters as a parameter, sB0,in the Hamiltonian of the solute: A
D
go= ly + ?(r) + AG(r,sBo) 2m1
The election of this parameter value is very important because it influences the dynamics of the solute. Different values imply different solute effective potentials. We considered the most appropriate value to be
the value obtained when the solvent is in equilibrium with the most probable position of the solute, although we also carried out calculations with other values of this coordinate. This limit corresponds to osol
