J. Phys. Chem. C 2007, 111, 2565-2572
2565
Dynamics of Water in Solvation Shells and Intersolute Regions of C60: A Molecular Dynamics Simulation Study Niharendu Choudhury* Theoretical Chemistry Section, Chemistry Group, Bhabha Atomic Research Center, Mumbai 400 085, India ReceiVed: October 20, 2006; In Final Form: December 12, 2006
We present dynamical characteristics of the solvent in various solvation shells of a C60 molecule as well as intersolute region between two fullerene molecules. Results are obtained from extensive atomistic molecular dynamics simulations of aqueous solutions of C60 in the isothermal-isobaric ensemble. Information about translational mobilities of water in and around fullerene molecules are obtained from velocity auto correlation function and mean square displacements as calculated from simulation data. Reorientational behavior also is investigated by analyzing first and second rank dipole-dipole correlation functions as calculated from simulation trajectories. Both the translational and rotational mobilities of water in the first solvation shell around a C60 molecule are found to be much slower than the same in the second and third solvation shells. Slower decay of the occupation time correlation function for the water in the first hydration shell indicates slower exchange of these water molecules with those from the bulk. Bimodal nature of the relaxations of both the reorientational and occupation time correlation functions is observed. Dynamics of water in the intersolute region between two C60 molecules also is investigated. Translational and orientational dynamics of water in the intersolute region are found to be even slower as compared to the same in the first solvation shell. Occupation time correlation function for intersolute water, however, decays faster than that of the first solvation shell water. In general, dynamics of water in the vicinity of a nonpolar molecule like C60 is qualitatively the same as that found near a biomacromolecule or micelle.
I. Introduction Water at the interfacial region in general and in solvation shells and intersolute regions of large molecules or molecular assemblies in particular have been the focus of many recent investigations.1-9 The ability of the interfacial water to influence structure, dynamics, function, and aggregation behavior of macromolecules or assemblies in aqueous solution has made their investigations very interesting. As for example, water in the vicinity of a protein molecule or in the confined region between two domains of a protein molecule plays a very important role in determining their native structure and function. Many biomolecules are even biologically inactive without a threshold level of hydration. Understanding the structure and dynamics of the interfacial water molecules is therefore of paramount importance in elucidating the mechanism of many processes and phenomena such as protein folding, molecular recognition, enzymatic reaction, self-assembly of nanomaterials, etc. Water in the immediate surrounding of a solute is generally called hydration water. If the solute is a biomolecule, the surrounding water is sometimes referred7 to as biological water. It has been shown in many recent studies4-7,9-11 that the structure and dynamics of the so-called biological water are distinctly different from the same of the bulk water. The hydration water stays for some time in the hydration shell of the solute and frequently gets exchanged with the same from the bulk. Thus, these water molecules are not as free as bulk water and at the same time not as immobilized as the interior water (e.g., crystallographic water of a protein molecule). The * Corresponding author. E-mail:
[email protected] or niharc2002@ yahoo.com.
hydration water being closer to the solute surface are well correlated with the solute molecule and because of this correlation nature of the dynamics of the water becomes significantly altered, particularly in the presence of an extended solute surface. Recently many experimental investigations6,12-18 have been devoted to understand the dynamics of water in solvation shells of various proteins, DNAs, and micelles. Many recent experimental techniques such as high resolution inelastic X-ray and neutron scattering (INS),14 infrared spectroscopy,15 and nuclear magnetic relaxation dispersion (NMRD)12,13 are capable of exploring the dynamics of the hydration water of a biomolecule. In particular, translational and rotational diffusions and residence time of hydration water have been investigated in many experimental studies using INS and NMRD. The NMRD technique has been shown to provide a more reliable and unified picture as compared to INS spectroscopy. Molecular dynamics (MD) simulations19 on the other hand can be used to explore the dynamics of hydration water around a large molecule. This method has the advantage of atomistically detailed description of all the molecular motions, and thus one can extract individual as well as collective dynamics from a MD simulation trajectory. Recently many computational investigations on dynamics of hydration water have been reported. Useful informations on the dynamics of the hydration water of proteins have been obtained from many simulation studies.11,20,21 It is now well recognized that the hydration water of a protein exhibits slower translational as well as rotational relaxations as compared to the bulk water. Recently a large number of simulation studies on the dynamics of water in the vicinity of a micelle has also been reported.22-24 Apart from the translational and rotational mobilities, hydrogen bond dynamics near
10.1021/jp066883j CCC: $37.00 © 2007 American Chemical Society Published on Web 01/25/2007
2566 J. Phys. Chem. C, Vol. 111, No. 6, 2007 hydrophilic solute surfaces have been investigated as well. The lifetime of a hydrogen bond near the hydrophilic solute surface involving a polar head group and water has been found23-25 to be much more as compared to the same between two water molecules in the bulk. In essence, water near an extended hydrophilic surface exhibits slower reorientational, translational, and hydrogen bond dynamics as compared to the same in the bulk. Another intriguing aspect of the hydration water is the existence of bimodal dynamics involving a slow and a faster relaxation processes. In some cases, the bimodal relaxation of the hydration water follows a stretched biexponential decay reminiscent of the same observed26 in supercooled liquids. In sharp contrast to the biomacromolecules or their assemblies with extended hydrophilic surfaces, there are many molecules or molecular assemblies that also have extended surfaces but with hydrophobic characters. Fullerenes, carbon nanotubes, and their assemblies fall into this category. These molecules have many interesting physical and chemical properties leading to potential applications27-30 in many areas as diverse as biology, chemistry, material science, and chemical engineering. The knowledge of dissolution and aggregation behaviors of these materials in aqueous solution will help material scientists and engineers in preparing nanomaterials with tailor-made properties of one’s interest. Recently, a concentrated aqueous solution of C60 has been prepared31-34 as a stable dispersion of nanoaggregates of C60 molecules. The mechanism by which these aggregates form is still not completely known. Many recent studies35-43 have been devoted to elucidate the mechanism of such hydrophobic association or aggregation processes in aqueous solution. The possibility of a dewetting induced collapse mechanism for the aggregation of nanoscopic or larger hydrophobic materials has been investigated. From the detailed computational investigations of behavior of water in and around large hydrophobic solutes, it has been demonstrated that a phase transition from the dewetted to wet state bridged by an intermediate empty-filled state can be observed just by tunning solute-water van der Waals’ interaction.42 The dynamics of water in the intersolute region between two model hydrophobic solutes of nanoscopic dimensions has also been investigated.40 Like hydration water of a biomacromolecule, water in and around nanoscopic or larger hydrophobic solutes have been found to play a major role in determining their aggregation and dissolution behavior in water. Therefore, it will be interesting to investigate the structure and dynamics of the water molecules at the C60 surface and in the confined region between two C60 molecules that may shed light on the stability of these materials and their aggregates in solutions. Many recent computational studies44-47 have reported a structure of water around a fullerene molecule. To reduce the computational requirements, coarse graining of the fullerene-fullerene48 and fullerene-water49-51 interactions have also been proposed. Although spatial and orientational structures of water around a fullerene molecule and in the intersolute region between two fullerene molecules have been studied in many recent investigations, 44-47,49 the dynamics of water in the solvation shell or in the intersolute region between fullerene molecules has not been investigated so far. In the present investigation, we intend to study the dynamics of water in various solvation shells of C60 using atomistic molecular dynamics simulation technique. We compare the dynamics of water at different solvation shells with each other. We are interested also in the dynamics of water in the intersolute region between two fullerene molecules. II. Models and Simulation Procedure Molecular dynamics simulation in the isothermal-isobaric
Choudhury TABLE 1: Different Spherical Shells of Water Considered around a C60 Molecule shell no
r1 (Å)
r2 (Å)
remarks
1 2 3
0 (or 6.0) 8.5 10.5
8.5 10.5 12.5
first solvation shell second solvation shell third solvation shell
(NPT) ensemble has been employed to study the aqueous solution of C60. In the present simulation study, we have considered a single buckyball as well as two buckyballs in water. A fully atomistic description of the buckyball has been used. Buckyballs are placed in a cubic simulation box containing ∼1000 water molecules. The overlapping water molecules from the immediate vicinity of the solute were removed followed by a steepest decent minimization. The water molecule was represented by the standard SPC/E52 model. Carbon atoms of the solute were represented by Lennard-Jones (LJ) interaction. The interaction between the oxygen atom of the water and the carbon atom of the C60 was also represented by LJ interaction with the energy and size parameters obtained by combining the individual parameters using Lorentz-Berthelot mixing rule. In the present study, we used LJ parameters for the carbon atoms of the buckyball σCC ) 3.47 Å, and CC ) 0.27647 kJ/mol as used in the original parametrization of Girifalco.48 The cross parameter for interaction CO as obtained from the present set of parameters is 0.10 kcal/mol, which is close to the value (CO ) 0.094 kcal/mol) adopted by Werder et al.53 to reproduce the macroscopic contact angle of a water droplet on a graphite surface. Details of simulation methods are described elsewhere.49 The solutes were rigid and kept fixed during simulation runs. All the systems were simulated at a target pressure of 1 atm and a target temperature of 298 K. Each of the calculated quantities was obtained from the postprocessing of a 1 ns trajectory saved at every 20 femtoseconds. For the investigation of the dynamics of water inbetween two C60 molecules, we have considered water molecules in the interstitial region between the two solutes by considering49 a cylindrical region enclosing both the solutes with the line joining their centers lying along the cylindrical axis. The dynamical properties of a fluid are generally obtained by calculating a time correlation function C(t) of the form N
C(t) ) 〈
∑ i)1
N
fi(t0)‚fi(t0)〉 ∑ i)1
fi(t + t0)‚fi(t0)〉/〈
(1)
where fi(t) is a vector function of positions or velocities of a molecule i at time t. Angular brackets denotes averaging over time origins t0. In the present analysis of time correlation functions for water in solvation shells or confined spaces between solute molecules, we adopt the convention that the entire dynamical history of a molecule is classified according to its positions at initial time. Thus, in the above equation the sum over N includes only contributions from those molecules that are within the specified region at an initial time, irrespective of there actual position at time t later or at any intermediate time between 0 and t. Therefore, eq 1 can be modified40 for the present case as N
〈 C(t) )
[fi(t + t0)‚fi(t0)]θi(t0)〉 ∑ i)1 N
〈
|fi(t0)|2θi(t0)〉 ∑ i)1
(2)
Molecular Dynamics Simulation Study of C60
J. Phys. Chem. C, Vol. 111, No. 6, 2007 2567 TABLE 2: Fitting Parameters for the Occupation Time Distribution Function R(t) as Calculated by Using Eq 3 and Fitted to a Curve According to Eq 4 shell no
A
τs
βs
B
τl
βl
1 2 3 intersolute
0.01 0.53 0.56 0.18
6.07 4.33 3.29 0.37
0.99 0.71 0.59 0.42
0.99 0.47 0.44 0.82
31.33 32.15 40.03 19.64
0.69 0.30 0.27 0.71
of interest at time t0 + t (or at time t0), otherwise fi(t0 + t) (or fi(t0)) is set equal to zero. Thus we can write a normalized occupation time distribution function R(t) as N
〈 R(t) )
θi(t0)θi(t + t0)〉 ∑ i)1 N
〈
Figure 1. (a, top) Solute-water radial distribution functions for the oxygen (gSO(r)) atoms of water. Three arrows with labels first, second, and third indicate positions of the outer boundaries of respective solvation shells. (b, bottom) Occupation time distribution functions of water in the solvation shells of a fullerene molecule. The dashed line is for the first solvation shell, the solid line is for second solvation shell, and the dotted line is for the third solvation shell. In the inset, we show the behavior of the same set of functions magnified at initial times.
where the function θi(t0) is 1, if the i-th molecule is within the region of interest at an initial time t0 irrespective of its position at any other time and is 0 otherwise. It is important to note that the dynamical behavior at short and intermediate time scales will be well characterized by this scheme, as one takes into consideration the crossing of molecules through the shell boundary.54 In this scheme, a molecule which leaves the region of interest can come back and contribute to the dynamics at any later time. At a long time however, dynamics may not reflect the effect of the confinement. III. Results and Discussions In this section, we present our results on various dynamical quantities of water in the solvation shells and in the intersolute regions of C60 molecules in aqueous solutions. A. Dynamics at Solvation Shells. The results for various dynamical properties of water in three different solvation shells around a fullerene molecule are presented. A solvation shell is defined as a spherical shell around the C60 molecule with inner and outer radii r1 and r2 respectively (see Table 1). The first solvation shell is considered to be the layer of water molecules nearest to the outside surface of the fullerene molecule and it essentially incorporates all the water molecules in the first major peak in solute-water radial distribution function, gSO (r) as shown in Figure 1(a). While second solvation shell consists of water molecules in the second gSO peak, third solvation shell contains water molecules beyond the second peak of gSO (r) with almost homogeneous distribution (see the locations of outer shell boundaries as shown by arrows in Figure 1(a)). 1. Occupation Time Distribution Function. The occupation time distribution function is calculated here from the time correlation function of the form represented by eq 1 with fi(t0 + t) (or fi(t0)) equal to 1 if a water molecule i is in the region
(3)
θi(t0)θi(t0)〉 ∑ i)1
where θi(t) is 1 if the oxygen atom of a water molecule i is in the region of interest at time t and is zero otherwise. The angular bracket in the above expression denotes averaging over the time origins. We have shown in Figure 1(b) R(t) as a function of time for water in the three solvation shells around the fullerene molecule. The function R(t) provides information about the local dynamics of the hydration water molecules. As seen in Figure 1(b), the occupation time distribution functions of the water in the second and third solvation shells decay initially very rapidly followed by a slower decay while the same for the first solvation shell, which is closer to the solute surface, decays slowly through out the entire time. The nature of these decay profiles indicate a bimodal character. Thus,at least two time scales are required to describe the relaxation processes over the entire period of time. Biexponential and stretched biexponential functions are generally used5 to describe hydration water dynamics of proteins and micelles. For the present case of hydration water around a nonpolar fullerene molecule, both the fast and the slow relaxation processes of the occupation time distribution function have been found to be well described by the stretched biexponential functions instead of the usual biexponential function. The relaxation processes of R(t) are thus represented by
R(t) ) A exp[-(t/τs)βs] + B exp[-(t/τl)βl]
(4)
where τs and τl respectively represent time constants for short and long time relaxations, while βs and βl are the corresponding stretched exponents. As shown in Table 2, the time constant for the fast decay process increases as the correlation due to the presence of the solute increases. However, the slower decay process shows the reverse trend. Relative proportions of the two processes are given by the preexponents A and B. It is important to observe (see Table 2) that for the first solvation shell the long time component dominates the decay with around 99% contributions, while proportion of the short time component is negligibly small. For the other two solvation shells, proportion of the two decay processes are almost equal. Short occupation time may be a result of the quick exchange of the water molecules lying near the shell boundaries. One of the primary features that characterizes the average occupation time in a particular solvation shell is the occupation half-life, the time at which half of the water molecules initially present in a solvation shell remain in the same shell. It is evident from these decay profiles that the half-life of the water in the first solvation shell is around 18.2 ps whereas that of the water in the second and third solvation shells are around 3.6 and 2.7 ps, respectively.
2568 J. Phys. Chem. C, Vol. 111, No. 6, 2007
Choudhury TABLE 3: Diffusion Constants Calculated from MSD as well as VCF and Fitting Parameters (Eq 6) for MSD of Water at Different Solvation Shells D/10-5 (cm2/s)
Figure 2. (a, top) Velocity-velocity autocorrelation functions for the water oxygens belonging to the three solvation shells of a C60 molecule. Inset shows the same plots magnified near the initial time over a smaller time window. Legends are the same as in Figure 1(b). (b, bottom) Mean square displacements of water oxygens for the first, second, and third solvation shells of a C60 molecule in its aqueous solution. Legends are the same as in Figure 1(b). Inset shows long time behavior of the same set of MSDs.
The longer half-life of the first solvation shell water indicates that the local environment of the water molecules is rigid as compared to the same in other solvation shells. The stretched exponential decay also known as Kohlrausch-Williams-Watts (KKW) law is often employed26 in supercooled liquids and recently it has been used5,20,24 to explain the residence times of the hydration water of proteins and micelles. It is important to note that KKW decay for the occupation distribution time correlation function is observed for the hydration water of a nonpolar nanoscopic molecule like C60. 2. Translational Mobilities and Diffusion. Translation mobilities of the water at different hydration shells as obtained from velocity autocorrelation function (VCF) and mean square displacements (MSD) are presented below. (i) VCF. We use the velocity autocorrelation function as a sensitive measure of the effect of the local environment of the water molecules in the solvation shells of the C60 molecules. VCF has been calculated for all of the three solvation shells and is shown in Figure 2(a). For the second and third solvation shells, VCFs are quite similar, while for the first solvation shell oscillations in the VCF are more pronounced because of the stronger caging effect arising from solute-water correlations. An initial small bump or shoulder observed for all the three states considered here is due to rebound of the water molecule from the shell of their neighbor.54 Diffusion constant, a measure of the mobility of the molecules, can be calculated from the integrals of the velocity autocorrelation function. Numerical integration of the VCF for all the three states have been performed with the upper limit of integration set to 1 ps. Calculated diffusion constants from VCF as well as MSD (discussed later) have been tabulated in Table 3. The diffusion constant of the water in the first solvation shell is much smaller as compared to same in the second and third solvation shells. (ii) MSDs. The MSD (〈∆r2〉) of the oxygen atom of the water molecules at different hydration shells for 60 ps has been shown in Figure 2(b). As all other correlation functions in the present
Parameters (eq 6)
shell no
VCF
MSD
P
a
1 2 3
2.55 3.32 3.55
3.11 3.31 3.45
1.19 1.70 1.83
1.11 1.04 1.03
work, MSD also has been calculated by considering only those water molecules that are in the specified solvation shell at the initial time. It is clear from this figure that slope of the MSD for the first solvation shell water is smaller than the other two cases in this time zone. However at a sufficiently long time when water molecules are far away from the solute surface, presence of the solute has a little effect on the water mobility and thus as expected, all three of the MSD plots shown in the inset of Figure 2(b) are co-incident. We have already seen that the occupation half-life of the first solvation shell water is around 20 ps, and beyond 40-50 ps the water molecule initially present in the first solvation shell is far away from the interface into the bulk and experiences little effect of the solute. The MSD is a good measure of the translational diffusivity of particles in a given environment. From the long time behavior of MSD, one can essentially extract diffusion constant (D) using Einstein’s relation, viz.
〈∆r 〉 1 D ) lim 6 tf∞ ∆t 2
(5)
where ∆t is the time difference corresponding to the measured MSD 〈∆r2〉. It is important to mention that extraction of D from the Einstein’s relation is possible when MSD shows linear time dependence at long time. In many recent studies, MSD of hydration5,24 and confined40 water have been shown to be nonlinear and therefore calculation of D from the Einstein’s relation may not be meaningful. However, here we obtain the diffusion constant from a linear fit of the MSD data calculated from simulation trajectory. The tabulated value of the diffusion constant for the first solvation shell as calculated from MSD differs from the same obtained from VCF (see Table 3). This minor difference between the D values calculated from the two different methods may be a consequence of the two different time scales considered in the two methods. It is important to note that the diffusion constant calculated from the MSD for the first solvation shell water is significantly larger as compared to the same calculated from the VCF. As already mentioned, calculation of the diffusion constant from MSD relies on the long time behavior. In the present case, as water molecules are free to move in and out of the hydration shell, the at long time possibility of crossing the shell boundary and moving to the bulk for a molecule initially present in the solvation shell is very high. Therefore, MSD at long time may have sufficient bulk character. However, calculated values of D from both of the methods show a similar trend in that the diffusion constant increases as one moves away from the solute molecule and eventually it should reach the bulk value far away from the solute. It is also important to note that the MSD of the solvation water is sufficiently nonlinear with time. To measure the departure from the normal behavior of the MSD, we have used
〈∆r2〉 ) PtR
(6)
Molecular Dynamics Simulation Study of C60
J. Phys. Chem. C, Vol. 111, No. 6, 2007 2569 TABLE 4: Fitting Parameters Obtained from the Fitting of the Rotational Correlation Function Γ1(t) to the Stretched Biexponential Decay Curve Given by Eq 8 shell no
A
τs
βs
B
τl
βl
1 2 3 intersolute
0.49 0.54 0.46 0.43
3.11 2.61 2.40 2.69
0.60 0.64 0.63 0.61
0.51 0.46 0.44 0.57
6.18 4.46 4.19 7.13
1.0 1.0 1.0 1.0
TABLE 5: Same as in Table 4 but for the Rotational Correlation Function Γ2(t)
Figure 3. (a, top) First-order orientational correlation function Γ1(t) (See eq 7 with l ) 1) of the dipoles of the water molecules in the three solvation shells around a C60 molecule in its aqueous solution. Inset: same plots showing fast decay at initial times. Legends are the same as in Figure 1(b). (b, bottom) Second-order orientational correlation function Γ2(t) (see eq 7 with l ) 2) of the water dipoles corresponding to the three solvation shells around a C60 molecule in its aqueous solution. Inset: same plots showing fast decay at initial times. Legends are the same as in Figure 1(b).
to fit our MSD data. The exponent R is a measure of the deviation from normal behavior. The value of R is 1 for a true Brownian diffusion observed in normal liquid. Unlike the sublinear behavior (R < 1) of MSD observed in many earlier studies,20,21,24,40 here we find the exponent R to be greater than 1. The values of P and R are given in Table 3. As expected, deviation from nonlinearity is maximum for the first hydration shell and decreases gradually toward unity as one goes through second and third solvation shells. On the other hand, proportionality constant P, which is a measure of diffusion, increases as the distance from the solute surface increases. Both the VCF and the MSD analysis show that translation mobility and diffusibility of water in the vicinity of the solute are hindered because of solute-solvent correlation. 3. Orientational Dynamics. Orientational dynamics of water in the solvation shells has been analyzed through the reorientational dynamics of the dipole moment vector µ. The time evolution of µ is generally determined by calculating the time correlation functions Γl defined as the autocorrelation function of the Legendre polynomials Pl(cos θ(t)), viz.
Γl(t) ) 〈Pl(µˆ (t0 + t)µˆ (t0))〉
(7)
where µˆ (t) is the unit vector along the dipole moment vector at time t. Here angular bracket represents both averaging over molecules and time origins. Like all other time correlation functions of this work, in the present case also we have considered only those water molecules that stay in the solvation shell at an initial time t0 through the introduction of the variable θi (cf. eq 2). The first (l ) 1) and second (l ) 2) order Legendre Polynomial are generally evaluated from MD simulation trajectory, as they can be obtained from experiments as well. Firstorder correlation function Γ1 can be obtained from infrared spectroscopy, whereas the second-order function Γ2 is obtained from NMR experiments. We present both first-order and second-
shell no
A
τs
βs
B
τl
βl
1 2 3 intersolute
0.44 0.33 0.35 0.61
0.24 0.08 0.09 0.54
0.34 0.40 0.42 0.33
0.56 0.67 0.65 0.39
2.32 1.64 1.66 2.67
0.86 0.86 0.90 0.85
order orientational time correlation functions Γ1(t) and Γ2(t) in Figure 3(a),(b) respectively. Comparison of the two functions shows that Γ2(t) in general decays much faster than Γ1(t) for all the three solvation shells. Both of the functions follow the same general trend that a longer time is required for the relaxation of the orientational correlation of water, which is in the close proximity of the solute as compared to those away from it. For both Γ1(t) and Γ2(t) of all the three states considered here, initially there is a very fast decay with a minimum or bump within 0.1 ps (see Figure 3(a) and (b)) followed by a slower decay with a longer time constant. The quick initial decay is due to librational oscillations and inplane rotations. Careful observation shows that initially time scales of all the three relaxation processes are almost the same because the effect of the difference in environments for the three states has not been manifested in this time domain. The nature of both Γ1(t) and Γ2(t) indicates that two different time scales exist and thus a biexponential decay may be required to fit the simulation data. However, it has been observed that Γ1(t) or Γ2(t) are not so well described by a biexponential decay and instead stretched biexponential decays describe both the functions quite well. Thus, the form of the equation with which both Γ1(t) and Γ2(t) have been fitted is as
Γl(t) ) A exp[-(t/τs)βs] + B exp[-(t/τl)βl]
(8)
where meanings of the symbols on the right-hand side of the above equation are similar to those of eq 4. The values of the time constants (τ), stretched exponents (β), and the relative proportions of the two components (A and B) for Γ1(t) and Γ2(t) are tabulated in Table 4 and 5, respectively. As expected, decays of both the longer and shorter processes in Γ1(t) as well as Γ2(t) becomes faster as one moves away from the solute surface. It is noted that the stretched biexponential decay of the rotational correlation functions has already been observed in the case of the hydration water of proteins5,20,21 and micelles24 and for the confined water between two nanoscopic hydrophobic plates.40 In general, we may conclude that the close proximity of the solute hinders the rotational or orientational motion of the water molecules. B. Dynamics in the Intersolute Region between Two C60 Molecules. In this subsection, we present our results on the dynamical behavior of water in the intersolute space between two C60 molecules and compare them with the corresponding results for water in the first solvation shells. In the present simulation study, two C60 molecules are kept in the middle of the simulation box containing a large number of water molecules. During simulation, two C60 molecules are kept fixed
2570 J. Phys. Chem. C, Vol. 111, No. 6, 2007
Choudhury
Figure 4. Occupation time distribution function of water in the intersolute region between two C60 molecules is shown by solid line. The same in the first solvation shell of a C60 molecule (dashed line) is also shown for comparison.
with an intersolute distance of around 12.66 Å corresponding to the first solvent separated state in the potential of the mean force45 of the C60 in water. The confined space between the two solutes is defined49 as a hypothetical cylindrical space enclosing the two C60 molecules with the length of the cylinder being the distance between the centers of two solutes and the radius being the same as that of a C60 molecule (for illustration see Figure 6(a) of ref 49). We consider the dynamics of those water molecules that occupies this cylindrical space at an initial time irrespective of their position at a later time t. Therefore, all the time correlation functions have been calculated using the variants of eq 2. The occupation time distribution function of the water molecules in the intersolute space along with the same of the first solvation water are shown in Figure 4. In each case, the nature of the decay is similar with a slower component at long time, but the initial decay is faster for the intersolute water. With the intersolute space being small, a large fraction of water occupying boundary region can easily move out of the confined space and thus contribute to the faster decay of the occupation distribution function. As for the solvation shell water, the bimodal nature of the relaxation is observed for the intersolute water also, and the relaxation process is well described by a stretched biexponential function. Time constants for both the faster and the slower components for the intersolute water are smaller than the respective quantities for the first solvation shell water (see Table 2). Thus, on an average the occupation time of water in the intersolute region is smaller as compared to the same of the first solvation shell water. The half-life of water in the intersolute region is around 7.2 ps as compared to 18.2 ps in the first solvation shell. The translational mobilities of water in the intersolute region also have been investigated by calculating MSD and VCF. The VCF and MSD for the intersolute water along with the same for the first solvation shell water are shown in Figure 5(a),(b), respectively. From both the VCF and the MSD plots, it is evident that the translational mobility of the intersolute water is less as compared to the same of the first solvation shell water. By definition of the time correlation function considered here, we impose restriction only on the initial position, and thus for a small intersolute region where average occupation time is of the order of a few pico seconds one can expect the molecule under consideration to be in the bulk within a very short period of time. The rotational or orientational dynamics of the water molecules in the intersolute region also has been investigated by calculating both first and second rank dipole-dipole correlation functions Γ1(t) and Γ2(t). These are shown and compared with the corresponding functions for the first solvation shell water and bulk water in Figure 6. Like the first solvation shell water
Figure 5. (a, top) Velocity-velocity autocorrelation function of water oxygens in the intersolute region between two C60 molecules is shown by the solid line. For comparison, the same for water in the first solvation shell around a C60 molecule also is shown (dashed line). (b, bottom) Mean square displacements of water oxygens in the intersolute region between two C60 molecules. The MSDs for water in the first solvation shell of a fullerene and in the bulk are also shown for comparison. The solid line, dashed line, and open circles represent MSD of water in the intersolute region, first solvation shell, and bulk, respectively.
Figure 6. (a, top) First-order orientational correlation function Γ1(t) (see eq 7 with l ) 1) of the dipoles of the water molecules in the intersolute region (solid line), in the first solvation shells (dashed line), and in the bulk (open circles). (b, bottom) Same as in Figure 6(a) but for the second-order orientational correlation function Γ2(t) (see eq 7 with l ) 2).
and bulk water,40 both Γ1(t) and Γ2(t) have two relaxation time scales. Relaxations of the intersolute water are seen to be slower than both bulk water and first solvation shell water. Orientational relaxation of the intersolute water also follow KKW decay with slower and faster time scales. Time constants and related quantities for the two relaxation processes as obtained from the
Molecular Dynamics Simulation Study of C60 fitting of Γ1(t) and Γ2(t) are shown in Tables 4 and 5, respectively. IV. Conclusions We have investigated the effect of the correlation of a nanoscopic nonpolar solute on the dynamical behavior of the water in the vicinity. Two types of confined environments have been investigated. In the first case, water in the various solvation shells around a C60 molecule have been considered, while in the second case we have studied confined water in the region between two C60 molecules in their first solvent-separated state. In general, translational and reorientational mobilities of water have been found to be slower for the water in the first solvation shell than in the other solvation shells. Both reorientational correlation function and occupation distribution time correlation function show KKW type relaxation. In the intersolute region, translational as well as reorientational dynamics are even more slower. Mean square displacements of water in both of the environments show anomalous nonlinear behavior. In essence, stronger is the correlation of the solute molecule with water and slower is the dynamics of water in general. It is important to note that slower translational mobilities and slower relaxations of the reorientational correlation functions also have been observed in the case of the hydration water of proteins and micelles. It is intriguing to observe such similarities in the dynamical behavior of the hydration water of two distinctly different classes of solutes. It will therefore be interesting to investigate the origin of slowness of the dynamics in these two apparently disparate systems. Excluded volume effect may play an important role in governing the dynamics of the hydration water. Investigations in this direction are in progress. Similarities in dynamics of solvation water of these two different classes of molecule also indicate that the effect of fullerene-water dispersion interaction may not have a considerable effect on the dynamics of neighboring water. However, in light of the recent observation53 that macroscopic contact angle of a water droplet on a model graphite surface depends on solute-water van der Waals interaction it will be interesting to investigate the effect of the solute-water dispersion interaction on the dynamics of solvation water. Acknowledgment. It is a pleasure to thank Dr. S. K. Ghosh and Dr. T. Mukherjee for their kind interest and encouragement. Thanks to the Computer Division, Bhabha Atomic Research Centre, Mumbai, for providing Anupam supercomputing facilities and support. References and Notes (1) Kauzmann, W. AdV. Protein Chem. 1959, 14, 1. (2) Protein-SolVent Interactions; Gregory, R. B., Ed.; Dekker: New York, 1995. (3) Kuntz, I. D.; Kauzmann, W. AdV. Protein Chem. 1974, 28, 239. (4) (a) Makarov, V. A.; Pettitt, B. M.; Feig, M. Acc. Chem. Res. 2002, 35, 376-384. (b) Lounnas, V.; Pettitt, B. M.; Philips, G. N., Jr. Biophys. J. 1994, 66, 601. (5) Bizzarri, A. R.; Cannistraro, S. J. Phys. Chem. B 2002, 106, 6617. (6) Halle, B. In Hydration Processes in Biology; Bellissent-Funel, M.C., Ed.; IOS Press: Amsterdam, 1999. (7) Bagchi, B. Chem. ReV. 2005, 105, 3197. (8) Nandi, N.; Bhattacharyya, K.; Bagchi, B. Chem. ReV. 2000, 100, 2013. (9) Pal, S. K.; Peon, J.; Zewail, A. H. J. Phys. Chem. B. 2002, 106, 12376. (10) Marchi, M.; Sterpone, F.; Ceccarelli, M. J. Am. Chem. Soc. 2002, 124, 6787. (11) Hua, L.; Huang, X.; Zhou, R.; Berne, B. J. J. Phys. Chem. B. 2006, 110, 3704.
J. Phys. Chem. C, Vol. 111, No. 6, 2007 2571 (12) (a) Polnaszek, C. F.; Bryant, R. G. J. Am. Chem. Soc. 1984, 106, 428. (b) Denisov, V. P.; Halle, B. Faraday. Discuss. 1996, 103, 227. (13) (a) Denisov, V. P.; Peters, J.; Horlein, H. D.; Halle, B. Nat. Struct. Biol. 1996, 3, 505. (b) Denisov, V. P.; Jomsson, B. H.; Halle, B. Nat. Struct. Biol. 1999, 6, 253. (14) (a) Bellissent-Funel, M. C.; Teixeria, J.; Brandley, K. F.; Chen, S. H.; Crespi, H. L. Physica B 1992, 740, 180-181. (b) Settles, M.; Doster, W. Faraday Discuss. 1996, 103, 269. (c) Zanotti J.-M.; Chen, S. H. Faraday Discuss. 1996, 103, 281. (d) Dellerue, S.; Bellissent-Funel, M. C. Chem. Phys. 2000, 258, 315. (15) (a) Zubavicus, Y.; Grunze, M. Science 2004, 304, 974. (b) Wernet, P.; Nordlund, D.; Bergmann, U.; Cavalleri, M.; Odelius, M.; Ogasawara, H.; Naslund, L. A.; Hirsch, T. K.; Ojamae, L.; Glatzel, P.; Pettersson L. G. M.; Nilsson, A. Science 2004, 304, 795. (c) Schultz, C. P. Nat. Struct. Biol. 2000, 7, 7. (d) Troullier, A.; Reinstadle, D.; Dupont, Y.; Naumann D.; Forge, V. Nat. Struct. Biol. 2000, 7, 78. (16) Halle, B.; Denisov, V. P.; Venu, K. In Biological Magnetic Resonance; Ramakrishna, N., Berliner, L. J., Ed.; Kluwer Academic/ Plenum: New York, 1999; Vol. 17, p 419. (17) (a) Otting, G. In Biological Magnetic Resonance; Ramakrishna, N., Berliner, L. J., Ed.; Kluwer Academic/Plenum: New York, 1999; Vol. 17, p 485. (b) Otting, G.; Liepinsh, E. Acc. Chem. Res. 1995, 28, 171. (18) (a) Sarkar, N.; Datta, A.; Das, S.; Bhattacharyya, K. J. Phys. Chem. 1996, 100, 15483. (b) Datta, A.; Mandal, D.; Pal, S. K.; Bhattacharyya, K. J. Mol. Liq. 1998, 77, 121. (19) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University: New York, 1987. (20) (a) Rocchi, C.; Bizzarri, A. R.; Cannistraro, S. Phys. ReV. E 1998, 57, 3315. (b) Bizzarri A. R.; Paciaroni, A.; Cannistraro, S. Phys. ReV. E 2000, 62, 3991. (21) Ghosh Dastidar, S.; Mukhopadhyay, C. Phys. ReV. E 2003, 68, 021921. (22) (a) Balasubramanian, S.; Bagchi, B. J. Phys. Chem. B 2001, 105, 12529. (b) Balasubramanian, S.; Bagchi, B. J. Phys. Chem. B 2002, 106, 3668. (23) Balasubramanian, S.; Pal, S.; Bagchi, B. Phys. ReV. Lett. 2002, 89, 115505. (24) (a) Ghosh Dastidar, S.; Mukhopadhyay, C. Phys. ReV. E. 2004, 70, 061901. (b) Chatterjee, C.; Majumder, B.; Mukhopadhyay, C. J. Phys. Chem. B 2004, 108, 7430. (25) Pal, S.; Balasubramanian, S.; Bagchi, B. Phys. ReV. E 2003, 67, 061502 and references therein. (26) Phillips, J. C. Rep. Prog. Phys. 1996, 59, 1133. (27) Taylor, R. The Chemistry of Fullerenes; World Scientific: Singapore, 1995. (28) Gorman, J. Sci. News (Washington, D. C.) 2002, 162, 26. (29) Cagle, D. W.; Kennel, S. J.; Mirzadeh, S.; Alford, J. M.; Wilson, L. J. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 5182. (30) Wharton, T.; Alford, J. M.; Husebo, L. O.; Wilson, L. J. In Recent AdVances in the Chemistry and Physics of Fullerenes and Related Materials; Kadish, K., Ruoff, R. S., Eds.; Electrochemical Society: Pennington, NJ, 2000; Vol. 9, p. 258. (31) Fortner, J. D. et al. EnViron. Sci. Technol. 2005, 39, 4307. (32) Deguchi, S.; Alargova, R. G.; Tsujii, K. Langmuir 2001, 17, 6013. (33) Alargova, R. G.; Deguchi, S.; Tsujii, K. J. Am. Chem. Soc. 2001, 123, 10460. (34) Scrivens, W. A.; Tour, J. M. J. Am. Chem. Soc. 1994, 116, 4517. (35) Pratt, L. R.; Pohorille, A. Chem. ReV. 2002, 102, 2671. (36) (a) Lum, K.; Chandler, D.; Weeks, J. D. J. Phys. Chem. B 1999, 103, 4570. (b) Chandler, D. Nature 2005, 437, 640. (37) Huang, X.; Margulis, C. J.; Berne, B. J. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 11953. (38) Ashbaugh, H. S.; Paulaitis, M. E. J. Am. Chem. Soc. 2001, 123, 10721-10728. (39) Choudhury, N.; Pettitt, B. M. J. Am. Chem. Soc. 2005, 127, 3556. (40) Choudhury, N.; Pettitt, B. M. J. Phys. Chem. B 2005, 109, 6422. (41) Choudhury, N.; Pettitt, B. M. Mol. Simul. 2005, 31, 457. (42) Choudhury, N.; Pettitt, B. M. University of Houston, Houston, TX. Unpublished work, 2006. (43) Choudhury, N.; Pettitt, B. M. J. Phys. Chem. B 2006, 110, 8459. (44) Walther, J. H.; Jaffe, R. L.; Kotsalis, E. M.; Werder, T.; Halicioglu, T.; Koumoutsakos, P. Carbon 2004, 42, 1185. (45) Li, L.; Bedrov, D.; Smith, G. D. Phys. ReV. E. 2005, 71, 6445.
2572 J. Phys. Chem. C, Vol. 111, No. 6, 2007 (46) Li, L.; Bedrov, D.; Smith, G. D. J. Chem. Phys. 2005, 123, 204504. (47) Hotta, T.; Kimura, A.; Sasai, M. J. Phys. Chem. B. 2005, 109, 18600. (48) Girifalco, L. A. J. Phys. Chem. 1992, 96, 858. (49) Choudhury, N. J. Chem. Phys. 2006, 125, 034502. (50) Vaitheeswaran, S.; Yin, H.; Rasaiah, J. C.; Hummer, G. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 17002.
Choudhury (51) Izvekov, S.; Violi, A.; Voth, G. A. J. Phys. Chem. B 2005, 109, 17019. (52) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. B 1987, 91, 6269. (53) Werder, T.; Walther, J. H.; Jaffe, R. L.; Halicioglu, T.; Koumoutsakos, P. J. Phys. Chem. B 2003, 107, 1345. (54) Zichi, D. A.; Rossky, P. J. J. Chem. Phys. 1986, 84, 2814.
