1392
J. Phys. Chem. 1996, 100, 1392-1405
Molecular Dynamics Simulation Study of Polarizable Solute Solvation in Water. 1. Equilibrium Solvent Structure and Solute Rotational Dynamics Badry D. Bursulaya,† Dominic A. Zichi,‡ and Hyung J. Kim*,† Department of Chemistry, Carnegie Mellon UniVersity, 4400 Fifth AVenue, Pittsburgh, PennsylVania 15213-2683, and NeXstar Pharmaceuticals, Inc., 2860 Wilderness Place, Boulder, Colorado 80301 ReceiVed: August 8, 1995X
Equilibrium solvation for polarizable, polar, and nonpolar solutes in water is studied via molecular dynamics computer simulations. A valence-bond electronic description for the solute is employed to allow for the instantaneous solute dipole readjustment to the fluctuating solvent environment. The effects of the solute electronic structure variation with the solvent configuration on solvationsin particular, the solute-solvent structure and solute rotational dynamicssare examined. A general similarity between the effects of increasing polarizability and of growing dipole moment is noted. This is probably due to the preferential solvation of higher dipolar states of polarizable solutes, arising from the nonpairwise additivity. However, the details as to the way and extent polarizability influences a physical observable vary with the characteristics of the latter. Thus a mere increase of the solute dipole moment cannot reproduce many different aspects of diversified polarizability effects. As for the solute-solvent radial distribution functions, a polarizable solute tends to make a structure compared to that of a nonpolarizable one with the same average charge distribution. Most affected are the hydrogen radial distribution functions. The solute polarizability tends to enhance its hydrogenbonding ability. Thus even a nonpolar solute can form weak hydrogen bonds with water molecules, depending on its polarizability and short-range repulsive terms. Solute orientational dynamics are also found to be significantly influenced by the solute electronic structure variation. Unlike the equilibrium radial distribution functions, however, it is the polarizability anisotropy that is responsible for altering the solute reorientational dynamics. With increasing polarizability anisotropy, the solute rotational dynamics become slower and the corresponding rotational friction grows. This increasing friction trend is attributed to coupling between the different components of the solvent electric field, induced by the solute polarizability anisotropy.
1. Introduction One of the key features that make solution-phase processes unique is the solute electronic structure variation with the surrounding solvent environment. Its crucial role in connection with various charge shift and transfer reactions is now wellappreciated. Among the best-known examples is an SN1 reaction class RX f R+ + X-. The solvent stabilization of the ionic state R+X- allows the heterolytic bond cleavage to occursa reaction pathway strongly suppressed in the gas phase; as a result, the reaction becomes drastically accelerated with increasing solvent polarity.1 While its importance is well-appreciated for the reactive processes, the solute electronic structure variation has usually been neglected in many theoretical investigations of nonreactive processes in solution. For instance, in many molecular dynamics (MD) studies of solvation dynamics and related dynamic electronic spectroscopy, the probe solute molecules are often assumed to be nonpolarizable, so that their electronic charge distributions are fixed regardless of the surrounding solvent. While this treatment is reasonable for weakly polarizable solutes, many molecular probes for actual spectroscopic measurements are characterized by large transition dipole moments and thus their electronic structure variation with the solvent configuration may not be completely negligible. In addition, since the excited †
Carnegie Mellon University. NeXstar Pharmaceuticals, Inc. X Abstract published in AdVance ACS Abstracts, December 15, 1995. ‡
0022-3654/96/20100-1392$12.00/0
states are often highly polarizable, proper inclusion of the solute electronic polarizability is required for an accurate modeling of real systems. Recent MD simulation studies indeed show that solvation dynamics strongly depend on solute polarizability. Solvation dynamics become progressively slower with increasing solute polarizability.2,3 It is also found that solvation dynamics (both equilibrium and nonequilibrium) become dependent on the solute electronic state and that the validity of linear response becomes significantly limited in the presence of a polarizable solute.2 These are attributed to nonlinearity arising from the instantaneous solute charge adjustment of the fluctuating solvent configuration. This, for instance, reduces the solvent force constant in the presence of a polarizable solute compared to a nonpolarizable one, so that the amplitude of the solvent thermal fluctuation becomes larger and its frequency lower with the former solute. Another important and interesting issue that needs scrutiny is how and to what extent the solute electronic structure variation can influence the surrounding solvent structure. Since the electronic polarizability leads to nonpairwise additive interactions, it could modulate the solvation shell structure and thus change, e.g., the hydrogen-bonding character in a nontrivial manner. (Indeed there is some evidence that the electronic polarizability and resulting nonpairwise additive interactions are essential to the correct prediction of the solvent distribution function and coordination number.)4 This also has potential consequences for reaction free energetics involving charge shift (e.g., electron transfer and SN1/SN2 reactions) because their transition states are often highly polarizable due to the crossing © 1996 American Chemical Society
Polarizable Solute Solvation in Water of two different electronic states. In this paper, we present our initial theoretical analysis of the solute polarizability effects on the solvent radial distribution functions. As we will see below, the solute polarizability can exert a significant influence on the solvation shell structure, in particular, the hydrogen distribution around the solute. Another property that can be influenced by the solute electronic structure variation is its transport behavior, for example, translational and reorientational dynamics. Due to their fundamental importance in many different areas of chemistry, the transport phenomena have been under intensive theoretical and experimental scrutiny. Among others, a great deal of attention has been paid to the solvent effects, in particular, the modifications of the transport coefficients due to the long-range solute-solvent interactions; dielectric friction has provided a useful theoretical concept for their understanding.5-12 Recent experimental findings provide revealing evidence that the solute rotational dynamics are significantly modulated by the solute-solvent electrostatic interactions and thus become state-dependent.12 In many previous theoretical investigations of rotational dynamics (either via a dielectric continuum solvent description or via MD simulations), a nonpolarizable solute is posited to analyze dielectric friction; as mentioned above, this completely neglects, for instance, the solute dipole moment variation with the fluctuating solvent, which could yield considerable modulation of the torque on the solute and thus its rotational dynamics. In view of the marked solute polarizability effects on solvation dynamics reported earlier2 and also an approximate linear relationship between the time-dependent dielectric friction and solvation dynamics,8a one would thus expect a non-negligible influence on the solute rotational friction by the polarizability. As we will see below, the polarizability anisotropy indeed affects solute rotational dynamics considerably. In the present paper, we examine the role of solute electronic structure variation on its structure and dynamics in solvent water via MD simulations. Our main goal is to understand and quantify the effects of the induced solute dipole moment present over and above the permanent dipole moment. To clearly discern effects due to explicit polarization, the parameters have been carefully adjusted, so that the average charge distribution of the polarizable models agrees with the fixed charge models. As in our previous study,2 the solute will be assumed to be polarizable in a quantum mechanical way, while the solvent is nonpolarizable. This type of approximation is sometimes referred to as QM/MM in the literature, where the semiempirical methods for the solute electronic structure calculations are combined with the molecular mechanics algorithms for the classical solvents.13 In order to expose the polarizability effects in a transparent way, we employ a simple electronic description for the solute couched in a few valence-bond (VB) states14-16 rather than a more quantum chemically oriented representation with a large number of basis functions. Thus in the current formulation the solute polarizability is realized as the solventdependent mixing of the VB states of differing electronic character. We do not address here the issues associated with the solvent electronic structure variation. While there exist various classical models for polar and polarizable solvents,17 we nonetheless believe that a quantum mechanical description is needed for proper account of both the inductive and dispersive effects, including their dynamic aspects. In fact, a recent theoretical investigation employing a quantum dielectric continuum description for the solvent shows that the dynamic character
J. Phys. Chem., Vol. 100, No. 4, 1996 1393 associated with the solute-solvent electronic structure variation can significantly alter the free energy predictions based on classical polarizability.18 The implementation of such a quantum mechanical solvent description within the MD context will be postponed for future studies. The outline of this paper is as follows. In section 2, we present a general theoretical formulation for a VB solute electronic description tailored for MD studies. The solvent coordinate formulation of the electronic free energy surfaces in solution is briefly reviewed. The details about the theoretical models employed in our MD studies together with the simulation methods are described in section 3. The MD results on the solute electronic structure, solute-solvent radial distribution functions, and solute rotational dynamics are discussed in section 4, while section 5 concludes. The analysis of solvation structure and dynamics will be reported elsewhere.19 2. Theoretical Formulation We begin by considering the electronic Hamiltonian operator H ˆ el of the solute-solvent system, employed in our MD studies. As mentioned above, we assume that the solute is electronically polarizable, while the solvent is not. In this QM/MM approximation, H ˆ el becomes2
H ˆ el ) H ˆ ° + Vˆ int + Vsolv
(2.1)
where H ˆ ° is a q-number solute electronic Hamiltonian in vacuum, Vsolv is a c-number solvent potential energy, and Vˆ int is the interaction operator between the two. (Thus H ˆ el corresponds to the full Hamiltonian of the solute-solvent system minus the kinetic energies associated with the nuclear degrees of freedom.) As for the solute electronic structure, we employ a truncated basis set description couched in a few VB functions. For convenience, we use the vacuum adiabatic VB basis B° ) {ψ°1, ψ°2, ..., ψ°n} (n ) number of VB basis functions), in which H ˆ ° becomes a diagonal matrix with the energy levels E°1, E°2, ˆ ° with an eigenvalue E°i. ..., E°n; i.e., ψ°i is an eigenstate of H The q-number Vˆ int consists of short-range repulsion-dispersion terms and long-range electrostatic interactions. The classical solvent description employed here implicitly posits that the formersin particular, the dispersion interactionsis independent of the solute electronic state. By contrast, the solutesolvent Coulombic interactions vary with the solute electronic state. In the point dipole approximation for the solute charge distribution,20 the solute-solvent interaction becomes
b Vˆ int ) VLJ - pˆ ‚
(2.2)
where the c-number VLJ denotes the short-range solute-solvent terms that do not depend on the solute charge distribution, pˆ is the solute dipole operator, and b is the electric field at the solute arising from the surrounding solvent. In general, pˆ has nonvanishing off-diagonal elements b µ °ij ) 〈ψ°i|pˆ |ψ°j〉 (i * j), which correspond to the transition dipole moments between vacuum adiabatic VB states ψ°i and ψ°j. Thus the off-diagonal components of pˆ ‚ b introduce solvent-induced electronic coupling among VB states, which is absent in vacuum. Therefore ψ°i ’s are no longer the eigenstates of the solute-solvent system; rather, the latter are given by the mixtures of the former
ψsa ) ∑caiψi°
(2.3)
i
where a is the solute electronic state quantum number in solution. The state coefficients cai’s determined by solving a due to the Scho¨dinger equation with H ˆ el (eq 2.1) change with b
1394 J. Phys. Chem., Vol. 100, No. 4, 1996
Bursulaya et al.
Figure 1. Schematic diagrams for (a) diatomic and (b) triatomic solutes, employed in the MD simulation studies. The LJ potential associated with each constituent atom is indicated as a sphere of radius σ/2. The mass point M in part (b) is completely embedded in the LJ spheres of the two E sites. It is slightly off-center by 0.505 Å from the line connecting the two E’s, so that the solvent molecules can approach the M site closer than the latter. The molecular coordinate system is chosen such that the three sites are on the yz-plane with the E-to-E direction parallel to the y-axis. Thus the x-axis is perpendicular to the molecular plane.
Figure 3. Probability distribution functions for the QP solute dipole moment in water. (a) Diatomic solutes: QP0 (-‚-); QP1 (-‚‚‚-); QP2 (s); QPS (‚‚‚). Except for the nonpolar QP0 solute, the probability distribution functions are quite asymmetric. (b) Triatomic solutes: QP3X (‚‚‚); QP3Y (s). Even though the average dipole moment and average polarizability for the QP3X and QP3Y models are identical except for the direction of the polarizability anisotropy, their dipole distributions are markedly different. The two small secondary peaks for QP3X, absent in QP3Y, are closely related to weak hydrogen bonds with water.
solution. The energy eigenvalues Esa for the state ψsa 2 ° Esa ) 〈ψsa|H ˆ el|ψsa〉 ) ∑cai Ei + ∑caicajb µ ij° ‚ b + VLJ + Vsolv i
Figure 2. Solute dipole moment variation along an equilibrium MD trajectory for the QP0 model. Due to the instantaneous readjustment of the solute electronic charge distribution to the changing solvent environment, its dipole moment fluctuates significantly around its average value 〈µ〉 ) 0 (cf. Table 2A).
solvent-dependent electronic coupling term pˆ ‚ b. As a result, the solute electronic charge distribution and dipole moment vary with b . We also note that the instantaneous solute dipole b µ a, equilibrated to the solvent
b µ a ) 〈ψsa|pˆ |ψsa〉 ) ∑caicajb µ ij°
(2.4)
i,j
depends on the gas-phase transition dipole moments b µ °ij (i * j) as well as the diagonal components b µ °ii. Thus even a nonpolar solute with vanishing b µ °ii’s can be momentarily polarized and become stabilized in solution (see Figures 2 and 3). For completeness, we briefly review here the solvent coordinate formulation of the electronic free energy surface in
i,j
(2.5)
define the Born-Oppenheimer (BO) electronic potential energy surfaces upon which the nuclear motions for the solute and solvent occur. These BO levels, in general, depend on all the nuclear coordinates of the system. Since we are interested only in the solvent influence on the solute electronic structure and resulting modulation of the BO surfaces, we simplify our description by projecting out those nuclear degrees of freedom which do not couple to the solute. Recognizing that the solute electronic structure depends only on b in the current dipole approximation, we integrate out all nuclear degrees of freedom except for b via an equilibrium ensemble average. This b) projection procedure defines a BO free energy surface Gsa( as
b) ≡ -kBT ln[∫dQ exp(-Esa/kBT)δ( b(Q) - b )] (2.6) Gsa( (Q) is where kB is Boltzmann’s constant, T is the temperature, b the solvent electric field variable that depends on the system nuclear coordinates, and ∫dQ denotes the configuration space integration. We note here that the ensemble average is taken in the presence of a given BO electronic state ψsa. Thus the
Polarizable Solute Solvation in Water
J. Phys. Chem., Vol. 100, No. 4, 1996 1395
electronic states are not thermally averaged. This projection reduces the dimensionality of each BO electronic surface to three, i.e., the three components of b , evaluated at the solute dipole location.21 These three components gauging collective b) in the point solvent motions are sufficient to describe Gsa( dipole approximation and are often referred to as solvent µij’s are all parallel to coordinates.22 If we can assume that b each other, Gsa can be further reduced to a one-dimensional curve in terms of the solvent electric field component |, which is parallel to the solute dipole.23 A similar projection of the kinetic energy terms will introduce a velocity variable ˘ |. The force constant24,25 and mass26,27 associated with | can be determined by utilizing the equipartition principle
ksa )
k BT
; msa )
〈δ|2〉a
k BT
)
msa
)
〈δ˘ |2〉a
∑
2〈ψi°|pˆ |ψj°〉〈ψj°|pˆ |ψi°〉
j(*i)
(2.9)
Ej° - Ei°
∑ b(*a)
2〈ψsa|pˆ |ψsb〉〈ψsb|pˆ |ψsa〉 Esb
-
Esa
)
2〈ψsa|pˆ |ψsb〉〈ψsb|pˆ |ψsa〉
∑ b(*a)
µ (D)
0.0 (0.3 (0.51 (0.9
0.0 3.6 6.1 10.8
solute
(q1 (e)
(q2 (e)
(qex (e)
E02 - E01 (cm-1)
R0a (Å3)
QP0 QP1 QP2 QPS
0.0 (0.45 (0.3 (0.16
0.0 (0.9 (0.9 (0.16
(0.41 (0.15 (0.3 (0.55
30 000 17 200 20 000 17 000
8.1 1.58 6.52 25.8
C. Quantum Triatomic Solutes QP3X QP3Y QP3IW QP3IS
R0,xx (Å ) R0,yy (Å3) R0,zz (Å3) R j 0 (Å3) Ei0 - E01 b (cm-1) 3
11.1 0.0 3.7 11.1
0.0 11.1 3.7 11.1
0.0 0.0 3.7 11.1
3.7 3.7 3.7 11.1
30 000 30 000 30 000 30 000
a
The vacuum polarizability along the solute molecular axis. b The energy gap between the vacuum adiabatic states ψ01 and ψi0 (i ) 2, 3, 4).
TABLE 2: Equilibrium Solute Electronic Structure in Water A. Diatomic Solutes solute
〈q〉a (e)
〈µ〉b (D)
R0c (Å3)
〈Rs〉d (Å3)
〈R j s〉e (Å3)
QP0 QP1 QP2 QPS
0.00 ( 0.025 0.51 ( 0.024 0.51 ( 0.106 0.47 ( 0.072
0.0 ( 0.30 6.1 ( 0.29 6.1 ( 1.27 5.6 ( 0.86
8.1 1.58 6.52 25.8
8.1 5.4 15.5 13.7
2.7 1.8 5.17 4.57
B. Triatomic Solutes solute
f
〈µx〉 (D)
〈µy〉f (D)
〈µz〉f (D)
〈R j s〉g (Å3)
QP3X QP3Y QP3IW QP3IS
-0.1 ( 1.70 0.0 0.01 ( 0.20 0.05 ( 1.70
0.0 -0.01 ( 0.47 0.00 ( 0.14 -0.01 ( 0.50
0.0 0.0 -0.01 ( 0.19 -0.01 ( 1.09
3.7 3.7 11.5
a
we define a similar tensor rsa in solution associated with ψsa as
rsa )
(q (e)
B. Quantum Diatomic Solutes
(2.8)
〈δ|2〉a
represents the promptness of the collective solvent response via inertial motion along |. In view of eqs 2.5-2.7, ksa, msa, and ωsa, in general, depend upon the solute electronic state.2,19 This is a direct consequence of the solute electronic structure variation with the solvent. For later use, we introduce polarizability tensors in vacuum and in solution. In view of the vacuum solute polarizability tensor ri0 for the state ψ°i
ri0 )
solute NP CN1 CN2 CN3
solute
x x ksa
A. Classical Diatomic Solutes
(2.7)
〈δ˘ |2〉a
where 〈 〉a means equilibrium ensemble average in the presence of a solute with ψsa and δ denotes the deviation of a variable from its equilibrium average value. The force constant ksa gauges the ease with which the solvent configuration can fluctuate on the electronic BO surface, while msa measures the degree of difficulty, or inertia, associated with the solvent velocity changes. Thus the solvent frequency ωsa given by
ωsa
TABLE 1: Parameters
Gsb
-
Gsa (2.10)
r0,s are the symmetric second-rank tensors, whose trace terms, in general, do not vanish. We thus reduce them into two irreducible tensors, i.e., isotropic polarizability r0,s and a traceless anisotropy term ∆r0,s
1 j 0,sI R j 0,s ≡ Tr r0,s; ∆r0,s ≡ r0,s - R 3
(2.11)
where I is an identity matrix. Since both the solvent-induced electronic coupling in eq 2.2 and r’s in eqs 2.9 and 2.10 involve pˆ , we see that the solute electronic structure variation with the solvent configuration will be closely related to its polarizability.28 Also due to the solvent modulations of the transition dipole moments and the relative free energies, the solute polarizability and its anisotropy vary with the solvent configura-
Average solute charge separation and its fluctuations in water. Average solute dipole moment and its fluctuations in water. c Vacuum polarizability along the solute molecular axis. d Average polarizability in water along the solute molecular axis. e Average isotropic polarizability in water. f Average solute dipole moment and its fluctuations along the molecular x, y, and z axes in water. g Average isotropic polarizability in water. b
tions; thus their solution-phase average values can be quite different from their vacuum values (see Table 2). 3. Models and Simulation Methods As in our previous work,2 the simulation cell consists of a single rigid solute molecule immersed in 256 rigid solvent molecules. The solute and solvent internal geometries are constrained via the SHAKE algorithm.29 The SPC/E water model is used to compute solvent interactions.30 We have studied both linear diatomic and nonlinear triatomic solutes, interacting with the solvent through Lennard-Jones (LJ) and Coulomb potentials. For all the diatomic solute models (Figure 1a), the LJ parameters, σ ) 4.0 Å and /kB ) 200 K, and mass, m ) 80 amu, are identical for each constituent atom, separated by d ) 2.5 Å. The partial charges on the solute molecule are varied within the diatomic models discussed below; the sites with a partial positive and a negative charge will be denoted as + and -, respectively.
1396 J. Phys. Chem., Vol. 100, No. 4, 1996 We have considered two distinct types of diatomic solute molecules, viz., a classical nonpolarizable (CN) model with fixed charges on the atomic centers and a quantum polarizable (QP) model with atomic charges that fluctuate with the surrounding solvent configuration.2 Four different CN models are explored; CN1 with q ) (0.3 e and µ ) 3.6 D, CN2 with q ) (0.51 e and µ ) 6.1 D, CN3 with q ) (0.9 e and µ ) 10.8 D, and neutral pair (NP) with vanishing partial charges and dipole moment. We have studied four different polarizable solute models (QP0, QP1, QP2, and QPS), whose electronic structure variations are described via two valence bond states, ψ°1 and ψ°2 (E°1 < E°2).2 The dipole moments of these two states as well as their transition dipole moment are represented by partial point charge pairs (q1,2,ex centered on the atomic sites. Thus, the model quantum diatomic solutes are anisotropically polarizable along the molecular axis (see triatomic solutes below). We have adjusted the parameters for the QP1 and QP2 models so that their average dipole moment in water is ∼6 D; this aVerage dipole moment coincides with the fixed dipole moment of CN2. To analyze and understand the polarizability anisotropy effects and their directional dependence, we have also studied the behavior for nonlinear triatomic solutes (Figure 1b), comprised of two identical LJ atoms E (each with m ) 60 amu, σ ) 4.0 Å, and /kB ) 200 K) and a third atom M of mass 40 amu with vanishing LJ parameters. They are in a bent geometry with C2V symmetry with a bond angle of 136° about atom M and an E-E separation of d ) 2.5 Å; the latter value is the same as the intersite separation for the diatomic solutes. Thus the site M is merely a mass point, embedded in the LJ spheres associated with the two E sites. As in the diatomic case, two distinct types of solute modelssi.e., classical nonpolarizable (CN3) and quantum polarizable (QP3)sare examined (the subscript 3 is to distinguish from the diatomic solutes). To avoid possible confusion due to the plethora of models, we confine ourselves to nonpolar triatomic models, so that none of the CN3 and QP3 solutes have a permanent dipole moment. Thus the nonpolarizable CN30 model is completely characterized by two LJ sites (E) and one mass point (M); these three sites will be referred to as real sites. As for the QP3 solutes, we employ four VB states, ψ°1,2,3,4, to model both isotropic and anisotropic polarizability; ψ°1 is the lowest-energy state among the four VB states, while the remaining three are assumed to be degenerate in energy.31 In order to describe the transition dipole moments among the four VB states, we introduce three more pairs of sites in addition to the three real atomic sites. The positions of these additional six sites are ((0.3, 0, 0), (0, (0.3, 0), and (0, 0, (0.3-0.337) (units: Å) in the molecular coordinate system with the origin located at M. Thus they are situated at the corners of an octahedron with its center at (0, 0, -0.337), which respects the C2V molecular symmetry, while breaking the axial symmetry associated with the LJ potentials of E. These sites are not real in that they are massless and will be referred to as fictitious sites. We have considered two isotropically (QP3IW, QP3IS) and two anisotropically (QP3X, QP3Y) polarizable models (Table 1C) to analyze the effects of polarizability anisotropy. For the former isotropic models, the transition dipole moments between the lower ψ°1 and higher ψ°2,3,4 states are represented by partial point charge pairs on the fictitious sites along the x, y, and z molecular axes, respectively. The transition moments between the three higher states are assumed to vanish. The QP3IS model is three times more polarizable than QP3IW. By contrast, the QP3X and QP3Y models are polarizable only along the molecular x and y axes (Figure 1b), respectively. Thus all the r0 components vanish except for R0,xx or R0,yy. We have chosen the parameters so that their isotropic
Bursulaya et al. TABLE 3: Solute Rotational Dynamics solute
〈µ〉a (D)
〈Rs〉 (Å3)
ζb
ζ(t ) 0)c
ζDFb
ζDF(t ) 0)c
NP CN1 CN2 CN3 QP0 QP1 QP2
0.0 3.6 6.1 10.8 0.0 ( 0.30 6.1 ( 0.29 6.1 ( 1.27
0.0 0.0 0.0 0.0 8.1 5.4 15.5
5.6 7.6 7.7 25.9 7.4 9.4 9.9
41.7 54.1 83.1 179.9 46.8 86.1 91.6
0 2.0 2.1 20.3 1.8 3.8 4.3
0 6.4 41.4 138.2 5.1 44.4 49.9
a Average solute dipole moment and its fluctuations in water. b Units: ps-1. c Units: ps-2.
polarizability (eq 2.11) is the same as that for QP3IW, while their respective nonvanishing r0 componentssi.e., R0,xx for QP3X and R0,yy for QP3Ysare the same as the corresponding elements for QP3IS. The simulations were conducted in the canonical ensemble at 298 K using the extended system method of Nose´.32 Periodic, truncated octahedral boundary conditions33 were employed. The central containing cube of the truncated octahedron is 25 Å in length, yielding a system density of 1.0 g/cm3. The long-range electrostatic interactions were computed with the Ewald method,34 resulting in essentially no truncation of these interactions. The trajectories were integrated with the Verlet algorithm35 and a time step of 2.0 fs. At each step of the QP model simulations, the lower adiabatic energy surface was followed by diagonalizing the Hamiltonian matrix in solution2,25 to yield the charges on the polarizable solute.36 The simulations were carried out with 20 ps of equilibration, followed by a 300 ps trajectory from which equilibrium averages were computed. 4. Results and Discussion The MD simulation results are summarized in Tables 2 and 3 and displayed in Figures 2-11. Some of the results are repeated for convenience and clarity. We first discuss the solute electronic structure in vacuum and in water. This is followed by a detailed analysis of solute-solvent structure and solute rotational dynamics. A. Solute Electronic Structure. The data summarizing the diatomic QP solute electronic structure in vacuum are presented in Table 1B, and the resulting ensemble averages in water are presented in Table 2A. We first notice that the average charge separation (〈q〉 and dipole moment 〈µ〉 ()〈q〉d) associated with the lower adiabatic state of the QP solutes in water are significantly enhanced compared to those in vacuum, except for the QP0 model. (Recall that the charge separation and dipole moment for the lower in-vacuo adiabatic state are given by (q1 and q1d.) The solvent environment responds to a nonvanishing solute dipole by becoming orientationally polarized. This in turn electronically polarizes the solute charge distribution further and thus increases its dipole moment compared to that in vacuum. (Here and hereafter, the electronic state index a will be dropped and the equilibrium ensemble average will be denoted simply as 〈 〉 because only the lowest adiabatic state is considered for the QP models.) We stress, however, that the QP solute charge distributions and dipole moments fluctuate in solution, including the QP0 model. Their fluctuation amplitudes, measured by the standard deviation, tend to increase with the solution-phase average polarizability 〈Rs〉. (Despite the same 〈Rs〉 value, the dipole moment fluctuations for the triatomic QP3X and QP3Y solutes are quite different (Table 2B). See below for details.) This is due to the instantaneous solute charge readjustment to the fluctuating solvent environment. The better the solute adjusts its electronic structure in solution, the more its dipole moment fluctuates with the solvent. In Figure 2, the QP0 dipole moment variation during the first 10 ps of its
Polarizable Solute Solvation in Water
J. Phys. Chem., Vol. 100, No. 4, 1996 1397
Figure 5. Hydrogen-solute radial distribution functions (a) gEH(r) and (b) gMH(r) for various nonpolar triatomic molecules. The QP3IS solute (- -) is highly and isotropically polarizable, while QP3X (‚‚‚) has strong anisotropic polarizability along the molecular x-axis. Both models exhibit weak hydrogen bond structures at the E and M sites, despite their vanishing average dipole moment in water. By contrast, the nonpolarizable CN30 solute (s) does not show any structure around 2 j r j 3 Å. The corresponding hydrogen g(r) functions for QP3IW with weak, isotropic polarizability and for QP3Y with strong, anisotropic polarizability along the molecular y-axis are almost indistinguishable from those for CN30 and thus are not shown.
Figure 4. Solvent radial distribution functions for the QP2 (s) and CN2 (‚‚‚) solutes: (a) g+O(r); (b) g-O(r); (c) g+H(r); (d) g-H(r). The NP solute distribution functions gnO(r) and gnH(r) are also plotted (- -) for comparison. The standard errors for g(r) are quite small; they fall in the range 0.01-0.15, depending on r.39
equilibrium trajectory is displayed. This figure graphically illustrates that, even though its average value vanishes, the instantaneous QP0 dipole moment can be considerably different from zero. It should be noticed there that the solute dipolesboth its magnitude and polarityschanges very rapidly; its typical time scale is much shorter than 1 ps. This is due to ultrafast solvation dynamics, reflected through the solute electronic structure variation (see Figure 11 below). Another interesting feature we observe in Table 2 is that the polarizability change upon solvation depends on the solute models we have investigated. To be more specific, 〈Rs〉 values for the QP1 and QP2 solutes are considerably higher than the corresponding R0 values, while the QPS model becomes less polarizable in water than in vacuum (i.e., 〈Rs〉 < R0). By contrast, 〈Rs〉 = R0 for QP0 (and also QP3 models). This arises from the fact that the solvent stabilization of each VB state strongly depends on its electronic character. In other words, if a particular VB state, say ψ°2, has a larger dipole moment than ψ°1, then the former will be better stabilized in water than the latter. This is exactly the case with the QP1 and QP2 solutes. Since ψ°2 is higher in energy than ψ°1 for these solutes, the energy gap between the two will be significantly reduced upon solvation. Since the polarizability is inversely proportional to the energy gap (eqs 2.9 and 2.10), it will be enhanced in water compared to vacuum. By contrast, the dipole moments for the ψ°1 and ψ°2 states are the same for QP0, i.e., µ°11 ) µ°22 ) 0. As
1398 J. Phys. Chem., Vol. 100, No. 4, 1996
Figure 6. (a) Crot 1 and (b) Cω(t) TCF’s for various classical nonpolarizable solutes: NP (- -); CN1 (-‚‚‚-); CN2 (‚‚‚); CN3 (-‚-). TCF’s are computed by averaging over two 300 ps trajectories. With increasing solute dipole moment, Crot 1 in part (a) becomes slower, while the corresponding Cω(t) in part (b) becomes faster. The statistical error estimated according to ref 51 is ∼0.06 for Crot 1 of CN1 around t ≈ 5 ps. Thus the crossing of the CN1 and CN2 rotational TCF’s and plateau-like behavior for CN1 do not seem to be real effects.
such, the solvent will not favor one state over the other and they will be equally stabilized. Thus the solute electronic structure will remain more or less unchanged upon solvation and so does the relative VB energy difference. Hence the polarizability will be little affected by solvation.37 The situation for the QPS solute is very similar to that for QP0 with one important exception; because of its finite dipole moment, the solvent electric field in equilibrium with QPS does not vanish. b, as This induces electronic coupling between ψ°1 and ψ°2 via pˆ ‚ pointed out in section 2 above. This results in two effects. First, the adiabatic gap in solution becomes larger than in vacuum. Second, due to mixing, the magnitude of the transition dipole moment between the two solution-phase adiabatic states becomes decreased compared to that in vacuum. From eqs 2.9 and 2.10, we see that both of these effects contribute to the polarizability reduction in solution. This rather differing polarizability behavior is a result of hyperpolarizability, i.e., polarizability variation with the applied external electric field strength. Since the solute dipole moment µ is a function of the external electric field, hyperpolarizability can also be interpreted as the µ-dependence of the solute polarizability. From Tables 1B and 2A, we see that the hyperpolarizability trend for the QPS solute is the opposite of those for QP1 and QP2. Both µ and polarizability increase upon solvation for the latter two solutes, so that ∂Rs/∂µ > 0, while their changes are in the opposite directions (∂Rs/∂µ < 0) for the former. These polarizability features can be further illustrated by calculating from the simulations distributions of the instanta-
Bursulaya et al.
Figure 7. (a) Crot 1 and (b) Cω(t) TCF’s for quantum polarizable solutes: QP0 (-‚-); QP1 (-‚‚‚-); QP2 (s). For comparison, the TCF’s for NP (- -) and CN2 (‚‚‚) are also exhibited. As the solute becomes more polarizable, Crot 1 tends to decay more slowly, while the Cω(t) relaxation seems to become slightly accelerated.
neous dipole moments along the molecular axis for each of these models. The resulting probability distribution functions for µ are displayed in Figure 3a for the QP0, QP1, QP2, and QPS models. The µ distribution function for QP0 is symmetric about µ ) 0. The direction of the instantaneous dipole induced by the solvent fluctuates in both directions along the molecular axis, since neither VB state has a net charge separation (cf. Figure 2).37 The remaining distributions all exhibit a net dipole moment along the molecular axis coincident with the gas phase direction. These distributions have similar 〈µ〉 values, 6.1 D for QP1 and QP2 and 5.6 D for QPS, even though their distributions are markedly different. Both QP1 and QP2 show µ distributions skewed toward larger magnitudes with respect to 〈µ〉, whereas the distribution for QPS is skewed toward lower magnitudes. This trend is closely-related to hyperpolarizability, mentioned above. For the QP1 and QP2 models (in particular, the latter), hyperpolarizability favors the higher dipolar solute states rather than the lower ones. As noted above, this arises from easier mixing of the two VB states ψ°1 and ψ°2 due to the decreasing energy gap, as the solvent electric field strength increases. By contrast, the lower dipolar states are preferred for QPS. The solvent polarization actually increases the energy gap, thereby making further electronic polarization more difficult, skewing the distribution of µ toward smaller values than the average. We also add a cautionary remark here. Although their skewness can be understood in terms of hyperpolarizability, it is by no means hyperpolarizability alone that completely specifies the µ distributions. As a revealing example, we consider the nonpolar, polarizable triatomic QP3X and QP3Y solutes. These two models are identical except for the polar-
Polarizable Solute Solvation in Water
Figure 8. Time-dependent frictions for the polarizable and nonpolarizable solute models: (a) ζ(t); (b) ζDF(t) (units: ps-2). QP0 (-‚-); QP1 (-‚‚‚-); QP2 (s); NP (- -); CN2 (‚‚‚). Regardless of the polarizability, ζ(t) for the dipolar solutes relaxes faster than for the nonpolar solutes. Within the models with the same average dipole moment, the initial value of ζ(t) increases and its initial temporal decay becomes faster with solute polarizability (anisotropy).
izability direction (Table 1C). Thus their hyperpolarizabilities are the same (again except for the directions). In Figure 3b, their probability distribution functions for µ are exhibited. Even though both models yield symmetric distributions about µ ) 0, their overall shapes are drastically different! QP3X has two small secondary peaks in addition to a main peak centered at µ ) 0, while the QP3Y distribution is characterized by a single peak that is quite narrow compared to the former. As a result, the dipole moment fluctuation for the former is more than three times larger than that for the latter, despite the same 〈Rs〉 (Table 2B). This clearly shows that more than just hyperpolarizability is involved in determining the µ probability distribution. In section 4.B below, it is argued that the interplay and competition between the short-range interactions and polarizability actually play a crucial role; the two secondary peaks for QP3X are related to weak hydrogen bonding with solvent water. B Solute-Solvent Structure. To study the polarizability effects on the solute-solvent structure, we have investigated the radial distribution function g(r) of water molecules around the solute. The g(r) functions of the water hydrogen and oxygen atoms with respect to the positive and negative sites of the CN2 and QP2 models are displayed in Figure 4. For reference, similar g(r) functions for the nonpolar and nonpolarizable NP solute are also plotted; the two neutral atomic sites of the NP solute are equivalent and will be denoted as n. We notice there that the permanent solute dipole moment plays a primary role in the modulation of the solvent structure. Since this is wellknown and has already been analyzed in connection with ionic systems in detail,38 it will not be considered here. Rather we
J. Phys. Chem., Vol. 100, No. 4, 1996 1399
Figure 9. Rotational dynamics for triatomic solutes: Crot 1 for the molecular (a) x and (b) y axes (cf. Figure 1b). CN30 (- -); Qp3IS (-‚-); QP3X (s); QP3Y (‚‚‚). As in the diatomic case, TCF’s are computed by averaging over two 300 ps trajectories. Despite its high isotropic polarizability, the rotational dynamics for QP3IS are faster than those for the less polarizable QP3X and QP3Y solutes.
Figure 10. Comparison between the electric TCF F|(t) for QP2 (-‚-) and CN2 (-‚‚‚-) and the normalized dynamic dielectric friction ζDF(t)/ ζDF(t ) 0) for QP2 (s) and CN2 (‚‚‚). Due to the neglect of the random projected dynamics8b and also due to the nonlinearity arising from solute polarizability,2 there is significant discrepancy between F|(t) and ζDF(t)/ ζDF(t ) 0). Though not shown explicitly, F⊥(t) for both CN2 and QP2 are very close to F|(t) for CN2. Thus the perpendicular electric field dynamics are little affected by the solute polarizability along its molecular axis.
will focus on the polarizability effects that are present over and above the permanent dipole moment effects. We first examine the oxygen distributions around the solute sites. The g+O(r) functions of both the CN2 and QP2 solutes are quite similar to gnO(r) of the NP, exhibiting a strong, broad first peak centered at 3.9 Å and a second, weaker peak at 6.6 Å. Compared to the corresponding NP one, the trend in CN2
1400 J. Phys. Chem., Vol. 100, No. 4, 1996
Figure 11. Dipole TCF (a) Clµ(t) and (b) ln Clµ(t) (l ) 1, 2) for QP2. The corresponding orientational TCF Clrot(t) values are plotted (s) for comparison. Due to a 20 fs time interval used to save the configurations generated by the MD simulations, the calculated Clµ(t) values at very short time are not smooth. Nonetheless, they clearly reflect the abrupt initial drop in F|(t) in Figure 10.
and QP2 is to shift the first peak to shorter distances, accompanied by a slight broadening, more evident in the QP2 function than the CN2 one, while their secondary peaks move out. This is a result of the electrostatic interactions of the positively charged solute site and the negatively charged oxygen atoms in water. Since the LJ parameter σ is fairly large compared to the water oxygen and the charge magnitude on the solute is not great, the structural change in the first peak is somewhat modest compared to that for a small cationic solute.38 g-O(r) functions show a marked decrease in the magnitude of the first peak along with a large shift toward smaller separations and a flattening of structure at larger distances compared to gnO(r). This trend is also more pronounced in the polarizable model and is attributed to weak hydrogen bonding of the water to the negatively charged site (see below). Similar to the oxygen distributions, the structure associated with the hydrogen distributions becomes altered for the polar solutes compared to the nonpolar one. The first peak in the g+H(r) function increases significantly while narrowing and shifting to a slightly larger distance. The enhanced structure in g+H(r) is again more pronounced in the polarizable QP2 model compared to CN2. The most pronounced structural effects due to introducing a partial charge pair on the solute are found in the g-H(r) functions. The CN2 and QP2 distribution functions both exhibit a well-resolved peak at 2.4 Å, which is absent in the nonpolar distribution; such a peak is characteristic of solute-water hydrogen bonding. The main peak centered around ∼4.5 Å in gnH(r) is reduced as a result of the hydrogen bond peak. Also there is an overall decrease in
Bursulaya et al. structure at separations greater than ∼6 Å. Again, we see that the change in hydrogen structure between the nonpolar and polar solute is enhanced for the QP2 model compared to CN2. The hydrogen-bonding peak of g-H(r) becomes higher for QP2 than for CN2 and shifts to smaller r by ∼0.05 Å. The main peak at 4.5 Å becomes a single broad peak for QP2, whereas the CN2 and NP solutes appear to have two peaks in the 4-5 Å region. The more pronounced polarizability effects observed for the hydrogen sites compared to the oxygen sites is due to the asymmetric charge distribution of the water molecule with respect to the center of its LJ sphere. The positively charged hydrogen atoms can approach the negative solute site more closely, to ∼2.4 Å, than the oxygen can approach the cationic solute site, ∼3.8 Å. Similarly, structural effects due to polarizability are more pronounced in the hydrogen pair correlation functions when comparing QP2 and CN2, although such effects are observed also in the oxygen-solute distribution functions. Since the CN2 and QP2 models both have the same 〈µ〉, the enhanced structure observed for the polarizable model is a result of structural changes associated with dipole fluctuations greater than 〈µ〉 compared to those accompanied by fluctuations to smaller values. This can be qualitatively explained by the following. Suppose a certain solvent molecule polarizes the solute and increases its dipole moment. The rest of the solvent molecules will respond to this enhanced solute dipole. Since the solvation free energy of a dipolar solute depends, to a good approximation, on its dipole moment squared, and the distribution functions for the solute are exponentially weighted by this free energy, the polarizable solute will exhibit structural effects which tend to mimic a nonpolarizable solute with a higher dipole moment. We also note that this is essentially a nonpairwise additive, many-body effect. All the distribution functions shown here display this trend. In addition, as seen in Figure 3a, the QP2 dipole moment is asymmetrically distributed, favoring µ greater than the average value. This will also tend to yield more structure in the QP2 model compared to CN2. Considering rather small statistical errors for g(r),39 the polarizability-induced structural enhancement observed here seems to be a real effect. Finally, we note that even the nonpolar polarizable solutes QP3IS and QP3X can exhibit weak hydrogen bonding to the electrically neutral LJ atomic E sites as well as the embedded site M. The hydrogen-solute distribution functions are displayed in Figure 5 for these nonpolar triatomic solute models as well as for CN30. (The corresponding g(r) functions of the weakly polarizable QP3IW and highly and anisotropically polarizable QP3Y models are almost indistinguishable from those of CN30 and thus are not shown explicitly there.) The hydrogen-bonding peak, located at ∼2.2 Å in gMH(r) and at ∼2.5 Å in gEH(r), is well-resolved in both the strong polarizable isotropic QP3IS and anisotropic QP3X solutes displayed. This is due to the transition dipole connecting the nonpolar VB states, absent in the CN3 solute. However, we also point out that the QP3Y model with the same Rs as QP3X does not exhibit at all the weak hydrogen-bonding peak present in the latter. This is due to the polarization direction coinciding with the E-E direction along the LJ sphere centers for QP3Y (Figure 1b). The LJ potential on the solute prevents water molecules from approaching the induced charge centers close enough for hydrogen bonding to occur. Due to different hydrogen-bonding character, the µ probability distributions for QP3X and QP3Y in Figure 3b exhibit very different structures. The two secondary peaks around µ ∼ (3 D present in QP3X correspond to dipolar states with large charge separation, which become stabilized through hydrogen bonding to water molecules.
Polarizable Solute Solvation in Water
J. Phys. Chem., Vol. 100, No. 4, 1996 1401
C. Solute Rotational Dynamics. As mentioned in section 1, solute rotational dynamicssin particular, dielectric frictionscan be influenced by its electronic structure variation with the fluctuating solvent environment. To address this, we briefly review the stochastic treatment of the solute rotational dynamics, usually analyzed via the orientational time correlation functions (TCF’s) Crot l (t)
Clrot(t) ) 〈Pl(uˆ ‚uˆ (t))〉 ) exp[-l(l + 1)DRt]
(4.1)
where uˆ is a unit vector representing the solute orientation, Pl’s are the Legendre polynomials, and DR is the solute rotational diffusion coefficient. The second equality in eq 4.1 follows in the diffusion limit attained at long time when the solute angular velocity ω b becomes equilibrated to its surrounding. In the generalized Langevin equation description, the equation of motion for ω b is given by5,7,8,11
dω b dt
b (τ) + B R(t) ) -∫0 dτ ζ(t - τ) ω t
(4.2)
where B R(t) is a random torque divided by the solute moment of inertia I and ζ(t) is the time-dependent friction given by
ζ(t) )
I 〈R B‚R B(t)〉 nrotkBT
(4.3)
Here nrot is the number of solute rotational degrees of freedom; thus, nrot ) 2 for the diatomic solutes, while nrot ) 3 for the nonlinear triatomic solutes. It is easy to show from eq 4.2 that ζ(t) is related to the solute angular velocity TCF Cω(t) by
1 ; Cω(t) ≡ 〈ω b ‚ω b (t)〉/〈ω b ‚ω b〉 C ˜ ω[z] ) z + ζ˜ [z]
(4.4)
where ˜f[z] is the Laplace transform of a time-dependent function f(t)
˜f [z] ) ∫0 dt e-ztf(t) ∞
(4.5)
Thus one can determine the total rotational friction ζ by integrating Cω(t)
ζ ) [∫0 dt Cω(t)]-1; ζ ≡ ∫0 dt ζ(t) ∞
∞
(4.6)
which is related to DR by
DR ) kBT/Iζ
does not obey the dielectric continuum theory predictions for point dipolar solutes; viz., ζDF increases quadratically with the solute dipole moment. Possible explanations for the discrepancy include the solvent molecularity,8b specific solvation, and the finite character of the solute charge distributions,10 as well as the breakdown of the perturbative aspects of various theoretical approaches.5,6 In Figure 7, Crot 1 (t) and Cω(t) for the QP solutes are compared with those for the nonpolarizable ones. The time correlation functions were calculated by averaging over two 300 ps trajectories. We notice that the solute polarizability has a pronounced effect on its rotational motions. (As we will see below, it is actually polarizability anisotropy that modulates solute rotational dynamics.) For instance, the Crot 1 (t) behavior for the nonpolar QP0 solute is closer to that for the dipolar CN2 model than to that for the nonpolar CN0. We observe that Crot 1 (t) becomes slower for the polarizable solutes than for the nonpolarizable models with the same average dipole moment, while Cω(t) for the former is faster than for the latter. Furthermore, as the solute polarizability increases, Crot 1 (t) tends to become slower. The ζ values for QP1 and QP2 determined via eq 4.6 are, respectively, 9.4 and 9.9 ps-1, while it is 7.4 ps-1 for QP0. Compared with the NP (ζ ) 5.6 ps-1) and CN2 (ζ ) 7.7 ps-1) solutes, we find that the polarizable solutes are characterized by higher frictions. We note that these trends with increasing solute polarizability are similar to those with increasing solute dipole moment, noted in Figure 6. To understand the increasing friction trend with polarizability, we have investigated the dielectric friction ζDF defined as5-12
ζDF(t) ≡ ζ(t) - ζNP(t); ζDF ) ∫0 dt ζDF(t) ∞
where ζNP(t) is the rotational friction for the NP solute. In view of Figure 4, the difference in both specific and nonspecific solvation between the polar and nonpolar solute is included in the definition of ζDF. The time-dependent friction was determined directly from the simulations by inverting eq 4.4.40 The results for ζ(t) and ζDF(t) are shown in Figure 8 and summarized in Table 3. We notice several interesting features there. First, the initial value of ζ(t) is higher for the dipolar solutes than for the nonpolar NP model. This is mainly due to the solutesolvent electrostatic interactions that contribute to the torque on the solute; this feature is absent for NP. To see this a little better, we approximate the initial value ζDF(t ) 0) from eq 4.3 as
(4.7) Crot 1 (t)
and Cω(t) are In Figure 6, the MD results for presented for various nonpolarizable solutes. We notice that, among the classical solutes, rotational motions become slower with increasing solute dipole moment, while the corresponding ω relaxation becomes faster. These trends arise from the increasing solute-solvent electrostatic interactions with the growing solute dipole moment, which tend to hinder the solute rotational motion and thus facilitate its rotational energy dissipation. From eq 4.6, the rotational friction ζ values are found to be 5.6, 7.6, 7.7, and 25.9 ps-1 for the NP, CN1, CN2, and CN3 solutes, respectively (Table 3); thus ζ grows with µ. This is a so-called dielectric friction effect. Since it has been extensively studied both theoretically5-11 and experimentally,12 we will not consider its details associated with the nonpolarizable classical solutes here. We will instead concentrate on the solute polarizability effects on ζ. We nonetheless point out here that, within the CN models we have investigated, dielectric friction, ζDF, in Table 3 (see eq 4.8 below for its definition)
(4.8)
ζDF(t ) 0) =
〈τ b‚τ b〉
; b τ )b µ×b
nrotIkBT
(4.9)
where we have assumed that torque due to the short-range interactions cancels exactly between the NP and polar solutes, so that only the long-range interactions make significant contributions. (Note that, at t ) 0, the random torque is the same as the deterministic torque.) The right-hand side of eq 4.9 can be evaluated directly from the MD simulations, independent of the ζDF(t) calculations. The results for ζDF(t ) 0) are 49.1, 47.8, and 52.2 ps-2 for the CN2, QP1, and QP2 solutes, respectively, while it is quite small (
