J. Phys. Chem. 1994,98, 62256230
6225
Sulfate Anion in Water: Model Structural, Thermodynamic, and Dynamic Properties William R. Cannon,'$tv#B. Montgomery Pettitt,t*#and J. Andrew McCammont*# Department of Biochemical and Biophysical Sciences and Department of Chemistry, University of Houston, Houston, Texas 77204 Received: February 3, 1994"
+
Binding energies, intermolecular distances, and partial charges from ab initio studies of Sod2- H2O were used to develop a molecular mechanics model for Sod2- in water. Structural, dynamic, and thermodynamic properties for the resulting model used in condensed-phase simulations agree well with experimental estimates. Thirteen waters are found to be present in the first solvation shell, and the residence time of these waters has been calculated to be 23 ps. The model anion has a free energy of solvation relative to xenon of -275 kcal mol-', which compares well with the experimental estimate of -266 kcal mol-'.
Introduction Of all the properties of ions in condensed phases, perhaps the most fundamental are the solvation properties. The detailed understanding of electrolyte solutions requires knowledge of the energetics and dynamics of ion solvation. Experimental data on equilibrium and dynamic properties of electrolyte solutions have been collected for many systems; cf. reviews by Marcus'J and Ohtaki and RadnaL3 Theoretical studies complement experimental work on electrolyte systems by providing microscopic models that aid in the analysis of the experiments4 and by permitting predictions of the properties of systems that have yet to be studied in the laboratory. Specifically,theoretical studies can provide detailed information on the energetic cost of creating a cavity in bulk water, subsequently filling that cavity with the ion and realignment of the solvent in the electric field of the ion. The structure and dynamic behavior of the solvent near the ion can be compared to the structure and dynamics of bulk water and water around other solutes. In this paper, we present the development of potential parameters for S042-,a tetrahedral dianion important in biological systems and in the Hofmeister series5 The sulfate model is also important because its simple geometry makes it a good choice for understanding fundamental pol yatomic ion-solventinteractions. Anions are sometimes assumed to be structure-breaking entities in water, in contrast to cations. However, s0d2-iS thought to be an exception.' We present here evidence of the organization of water around SO4*-in a consistent manner. The paper is organized as follows. We first present the developmentof potential functionsbased on quantum mechanical and model gas-phase binding studies, followed by simulations to obtain the free energy of hydration of S042-relative to that of xenon. The resulting model is used in molecular dynamics simulations to provide structural, thermodynamic, and dynamic information on dilute solutions of so42-in water. Finally, the results are discussed in relation to experimental data.
calculated, and basis set superposition errors were determined by the counterpoise method.6 Optimization of the sulfate dianion geometry using HartreeFock calculations and only sulfate atomic centers gives a bond distance in the tetrahedral dianion of 1.4898 A, in agreement with the small molecule crystallographic value of 1.49 A.7 Calculationson bimolecular complexes with H20 were performed with S02-fied in its isolated optimized geometry and water held in a geometry of r = 0.9572 A and OHOH = 104.52O. Charge distributionsweredeterminedfrom the quantum mechanical wave functions by fitting the electrostatic potential to point charges centered on each of the atoms using the CHELPG*T~ facility in Gaussian 92.1° Three low-energy geometriesof the S042--H20complex were studied, as shown in Figure 1: (a) the water coordinated in a bidentate fashion with each hydrogen interacting with a single sulfate oxygen and the water lying in the plane defined by these two sulfateoxygensand the sulfur atom; (b) the water coordinated such that the interacting hydrogen and water oxygen are collinear with a sulfate oxygen and sulfur, and the other water hydrogen is coplanar with the same sulfate oxygen, sulfur, and a third sulfate oxygen; (c) a configuration identical to (b) except that the water is rotated 180' about the four-atom collinear axis. In addition, the interactionenergy at the lowest-energyconfiguration, Figure la, was determined at the HF/6-31 l+G(d) level, and distance-dependent interactionenergies of the same configuration were determined at the HF/6-31+G(d) level. Potential Parameters. Molecular mechanics potential parameters were developed based on the results of the ab initio calculationsand previous model parameters from thermodynamic simulations. Specifically, care was taken to make them as consistent as possible with the OPLS parameter set by utilizing the Lennard-Jones parameters for OPLS thiols, sulfides, and disulfidesll and the OPLS CH30-.12 The intermolecular interactions aredescribed by Coulombicand Lennard-Jonesterms shown in eq 1, where the summation is over all intermolecular atom pairs.
Methods
Ab Initio. Molecular orbital calculations using the 6-31+G(d) basis set were used to determine the ion geometry, optimized intramolecular distances, charge distributions, and interaction energies with water using the supermolecular approach. Correlation energiesup to the Mraller-Plesset MP4SDTQ level were t Department of Biochemical and Biophysical Sciences. t Department of Chemistry. e Abstract published in Advance ACS Abstracts, May 15, 1994.
0022-365419412098-6225%04.50/0
Based on the results from the ab initio calculations shown in Figure 1, little charge is transferred from the soluteto the solvent. Thus essentially the only unknown parameter in the potential function development is the charge distribution within the anion. Two models summarized in Table 1 were developed based on agreement with the ab initio binding energies and optimized distances. The first uses previously established Lennard-Jones 0 1994 American Chemical Society
6226 The Journal of Physical Chemistry, Vol. 98, No. 24, 1994 TABLE 1: Potential Parameters for SO,% model atom charge u,A 2.400 3.55 std 1 S 0 std 2
S 0
-1.100 2.000 -1.000
3.15 3.55 3.15
e,
kcal mol-' 0.250 0.250 0.250 0.200
parameters, and only the charge distribution is varied to produce agreement with the ab initio results. The second model utilizes essentially the charge distribution found from fitting to the quantum mechanical electrostatic potential for the S042-.H20 complex shown in Figure l a and allows for minor adjustment of the Lennard-Jones parameter e. Although decreasing e decreases the depth of the well of the Lennard-Jones potential, it also affects the slope of this potential as well as the slopeof the overallpotential. Thus, depending on the Coulombic interaction, the energy calculated from eq 1 can either decrease or increase when e is decreased. Molecular Dynamics and Thermodynamic Integration. Computationally, the free energy of solvation is obtained most conveniently through either free energy perturbation or thermodynamic integration in which one molecular speciesis changed into another species.l3 In either approach, the Hamiltonian H(X) is changed from that of one system (X = 0) to another (A = l), and the corresponding change in free energy is computed. The fundamental equation of thermodynamic integration at constant temperature and pressure is14
In the application of this equation in computer simulations the ensemble average of the derivativeof the Hamiltonian is collected at discrete values of the integration variable A. A quadrature over these values gives the difference in the Gibbs free energy of the two systems of interest. In order to obtain the relative free energy of solvation of S042to xenon, a thermodynamic cyclel5 is used in which xenon is slowly changed into S042- in solution and in uacuo. Initial equilibration on a solution of sulfate in water was carried out in which both the S and 0 atoms of SO1were uncharged in a periodic boxof 894 SPC/E waters at 1 atmusing theARGOS'6molecular dynamics package. The water positions were taken from a preequilibratedboxof water, and waters that madecontacts closer than 2.8 8,with the solute were removed. The system pressure and temperature were maintained by volume and velocity scaling with a relaxation time of 0.1 ps each." The cutoff distance used was 14 8, for all interactions. A time step of 2.0 fs was used along with the SHAKE algorithm for constraining all bonds. Equilibration was carried out for 5 ps at each of the temperatures 20, 70,120, and 200 K and for 35 ps at 298 K. Convergences of the total system energy and potential energy were used as indicators of equilibration of the system. Both the total energy and the potential energy converged within 20 ps in the simulation at 298 K. Starting from the uncharged sulfate, multiconfiguration thermodyanmic integration14 to sulfate with a -2 charge on the sulfur was performed in 51 steps of the integration variable X in which each step of X included 1000 equilibration steps and 1000 data gathering steps. Using a time stepof 2 fs, the total integration time was 204 ps. Again, bond lengths were constrained by the SHAKE algorithm. Temperature and pressure were restrained to 298 K and 1 atm by velocity and volume scaling with relaxation times of 0.1 ps. Next, the final configurationwas used as a starting structure for 30 ps of MD in which kinetic, potential, and total energies were monitored again. Subsequent thermodynamic integration to model 2 (Table 1) was performed in 21 steps of X with 1000 equilibration and 1000 data gathering steps for a simulation time of 84 ps. A cutoff of 14 A was again used for
Cannon et al. all interactions. Subsequent molecular dynamics of the system was carried out for 100 ps in order to obtain information on the hydration structure and dynamics of the dianion. Integration from the initial uncharged sulfate to xenon was also performed in 21 steps of the integration variable A. The xenon model used was that of Straatsma et for which the hydration free energy is known. Ensemble averages at each step in X were collected using 1000 equilibration steps and 1000 data gathering steps at 2 fs/step for a total simulation time of 84 ps. During the integration the S-0 bond was gradually shortened from the initial 1.49 to 0.70 A. Constraint contributions to the free energy were evaluated a n a l y t i ~ a l l yand , ~ ~contributions due to changes in the moment of inertia as the bonds were shortened were accounted for as the work needed to keep the kinetic energy of the system constant.lg Cutoff distances were again 14 8, for all interactions. This procedure was employed for simulations in solution and in vacuo. The characteristic residence time of water in the first solvation shell of S042- was calculated according to20
in which n(0) is the number of waters in the first solvation shell at time t = 0, n(t) is the number of initial waters still present in the first solvation shell at time t , and T,f is the characteristic residence time for waters around the ion. The number of waters around the ion at a time T that were present at some initial time t was calculated from the autocorrelation functionzl
(4)
Here T is given by T = kA7, where k = 0, 1,2, ... (tmax/AT) and AT is the time step between observations. We used AT = 3 ps in order to exclude contributions from solvent molecules that temporarily leave the first solvation shell and return without ever becoming incorporated into the bulk solution. The results were not found to be sensitive to small changes in this value. The values n ( T ) at various times are then fit to eq 3 by weighted least squares.
Results and Discussion Gas Phase. Interaction energies and optimized intermolecular distances determined by ab initio calculations are given in Table 2. The interaction between the anion and water is expected to bedominated by Coulombic and inductiveeffects. For the lowestenergy configuration studied (Figure la), the effect of including correlation at various levels contributes 0.4 (MP4SDTQ) to 0.8 kcal mol-I (MP3), which amounts to 2-3% of the total interaction energy. The effect of including correlation is even smaller for the non-minimum-energy configuration shown in Figure 1b, for which correlation accountsfor 0.7% of the total interaction energy. The effect of basis set superposition error (BSSE) is more dramatic, accounting for 13% of the total interaction energy at the MP4SDTQ level. The charge distribution found from fitting to the quantum mechanical electrostatic potential for the configurations shown in Figure lb,c is very similar to the charge distribution found for the configuration shown in Figure la. For comparison, determination of gas-phase S042- charges by fluorescence methods22 resulted in qo = -0.93. Binding energies and intermolecular distances determined by the molecular mechanics models are shown in Table 3. Both models perform well in reproducing the ab initio binding energy at the minimum-energy configuration and adequately reproduce energies for other configurations. The two configurations for which the models do not quite reproduce the ab initio results are
The Journal of Physical Chemistry, Vol. 98, No. 24, I994 6227
Sulfate Anion in Water
Effect of Correlation and BSSE on Binding Energies (Corrected for BSSE) for Configurations shown in Figure 1’ method distance basis set HF MP2 MP3 MP4D MP4DO MP4SDTO configuration a -24.87 -25.65 -25.68 -25.45 -25.24 -25.29 6-3l+G(d) 3.44 1.305 3.026 2.840 2.904 2.812 3.326 BSSE -26.03 (-26.45) 6-311+G(d) -18.07 6-3l+G(d)b 4.44 -1 1.65 6-3l+G(d)b 5.44 configuration b -19.34 -19.82 -19.78 -19.61 -19.49 -19.48 6-3l+G(d) 4.21 configuration c 6-31+G(d) -19.34 4.21 a Parenthetical entry assumes additivity of correlation effects from HF/6-31+G(d) to HF/6-31 l+G(d) level. Binding energies are in kcal mol-’. b The basis set superposition error was not accounted for in these calculations.
TABLE 2
~~
~
~~
(-? -0.98
i)
TABLE 4 Thermodynamic Results for the Solvation of
-+
0.55
- o . 9 p . m a
in Water Relative to Xenon. thermodynamic experimental Printegration estimate 3.4 f 0.1 sow xew -6.2 f 0.3 XY,, SOYd -259.9 f 0.8 soy,, szi) 34.1 f 0.8 Szi)4(@ so$;:) AAG -275 f 2 -266.1 *SO, designates a neutral species with no partial charges, S2-0, designatesa specieswith charge centered on sulfur,and S042-deSignateS the final sulfate model 2 with distributed charges. Simulation total includes a -46.8 kcal mol-’ Born cutoff correction.
!304%
O q
-0.97
0.39 -1.00
b -0.99
0.39
-0.99 C
Fipre 1. Configurations of S042--H20complexes studied in ab initio
and model binding studies: (a) water dipole aligned with sulfate electrostatic field; (b) water 0-H bond collinear with sulfateS-0 bond; (c) same as (b) but rotated 180° about collinear axis. TABLE 3 Potential Function Results for SO4%-.H2Ofor Both Potential Models‘ model 1 model 2 configuration rs,o Ebind rs.0 EW a 3.00 -16.33 3.00 -17.69 3.32 -26.04 a 3.34 -26.12 a 4.44 -15.92 4.44 -15.82 -10.29 -10.31 5.44 a 5.44 b 4.24 -20.18 4.24 -19.59 -20.61 -2 1.09 4.09 C 4.09 a Distances at 3.00,4.44, and 5.44 A are fixed while all other distances are the minimum-energy distance for that configuration. the configurations for which BSSE was not accounted for in the ab initio calculations. One source of single ion experimental data that would be valuable for comparison is that from gas-phase binding studies by mass spectrometry. However, ions such as S042-are difficult to study in this manner because of the potential for gas-phase redox reactions. We are not aware of any gas-phase binding studies on S042--H20. Condensed Phase. In any development of potential parameters for use in liquid simulations, care must be taken to make the condensed-phase properties of the molecule of interest as realistic as possible. Properties obtained from simulations can be compared to those obtained from experiment or other computational
methods. Evaluation of the thermodynamic properties of single ions in solution by experimental methods is difficult, however. In the laboratory, single ions cannot be solvated in bulk independently of their counterions for unambiguous determination of thermodynamic properties. A number of methods’ have been used in an attempt to divide up system properties between cation and anion solutes. The free energy of solvation for S042-can be derived from estimates of the enthalpy of s~lvationlJ’-~sbased on several methodsW26 including computation of lattice energies via the Born cycle and knowledge of the entropy of solvation of the proton.27 Estimates of the standard free energy of solvation can vary widely. The simulation results are compared to the free energy of solvation value reported by Marcus’ of -260.5 kcal mol-’ for SO4” and the free energy of solvation value for xenon reported by Ben-Naim28 of 1.3 kcal mol-’. To the former value needs to be added a correction for differences in standard states28 using AG-1, = AGx - RT ln(pi/pD, where p: is the number density of the solute in the liquid phase and pf is the number density of the solute in the gas phase, AGwlv is the free energy of solvation as defined by Ben-Naim, and AGx is the free energy of solvation for the process in which the standard state in the gas phase is the ideal gas a t 1 atm pressure and the standard state in the condensed phase is the hypothetical 1 M ideal aqueous solution. With RTln(pS/pD = 4.3 kcal mol-’, a value of -264.8 kcal mol-’ is appropriate for the solvation free energy of the anion. The free energy of solvation difference between xenon and sulfate is then -266.1 kcal mol-’. The thermodynamic integration results in a total free energy of solvation difference between sulfate and xenon of -275.4kcal mol-’ when a Born cutoff correction29of -46.8 kcal mol-’ is used for the solvation of a 14-A sphere. The individual free energy results for the processes of cavity formation, charging, and distributing the charges are shown in Table 4. Convergence of the solution simulations is shown graphically in Figure 2 where the zero-point energy shown is the value of the accumulated free energy at X = 1 and 1000 data acquisition steps for each value of X. The free energy change is dominated by the affect of charging the neutral species to -2, as expected. Distributing the charge in the final S042- model requires 34 kcal mol-’ of work for
6228 The Journal of Physical Chemistry, Vd. 98, No. 24, 1994
Cannon et al. s-w - 0-w __ w-w_ .
4.0 .
0.0
E
-2.5
-5.0
200.0
0.0
-
-
400.0 600.0 Dynamic stepsnmbda
800.0
1OOO.O
-
Figure 2. Convergence of thermodynamic integrations for neutral SO4 Xe, neutral SO4 S2-04,and S2-04 SO42-. The abscissameasures the number of data gathering steps used in obtaining averages in eq 2.
r (angstroms)
Figure 3. Pair correlation functions for S-O,, O,-O,,
and O,-O,
pairs.
TABLE 5: Comparison of Structural and D ~ M I U ~ C Properties and Relative Partial Molar Volumes MD Simulation and Experiment MD
diffusion constant
cm2s-*)
exwriment
1.5 k 0.6
Figure 4. Sterwview of SO42- and first solvation shell waters. 13.25 3.15 2.68
6.0-14.33 3.1-3.93 2.79-2 .883a3a,34,xis,si
reorganization of the solvent around S042-. Although the 3.3% difference between the experimental value and the simulation value is reasonable, it is important to bear in mind that approximationsare made in both the experimental determination of singleion free energies and free energies of ions from simulation. As discussed above, the experimental values for single ions are difficult to determine and can vary substantially from method to method; the simulations suffer from inaccuracies due to longrange interaction cutoffs,30the cutoff free energy correction from the Born e q u a t i ~ n ,and ~ , ~parametrization ~ and functional form of the solvent and solute models. The close agreement between simulation and experiment may, therefore, be fortuitous. On the other hand, the model S042- parameters were not fit to the experimental free energy data-the two results are independent of each other. It would be desirable to study the process of solvation with an explicitly polarizable water model. However, it was estimated that a single free energy integration of the uncharged SO4to the charge centered model in polarizable water31 alone would require 85 h on a Cray Y-MPusing a 14-A cutoff, even with highly efficient vectorized code. This was calculated assuming that an integration would be performed in 5 1 steps of X with 1000 equilibration and data gathering steps each. It is possible that a less conservative integration could be used which would decrease costs; however, the addition of polarization might add other structural relaxation mechanisms and time scales. The model Sod*- properties are also compared with other dynamic and thermodynamic properties found from experiment in Table 5 . Estimates of the size of the anion and of the surrounding solvation shell are available from X-ray diffraction data on numerous sulfate salt solutions including (NH4)2S0432-35 and ZnS04,36 The diffraction studies show that the first peak in the correlation function of the sulfate solution near 2.8 A is due to a combination of the anion-oxygen, water-oxygen pair correlation and water-oxygen, water-oxygen correlations and a second peak at approximately 3.7-3.8 A is due to the sulfur, water-oxygen pair. The pair correlation functions from the simulation are shown in Figure 3. The first peak in the S-0,
3.0
2.0
-SOe7HPO
A
_ _ _ _ S04-14H20
1.o
0.0,
3.0
4.0
5.0
8.0
7.0
0
r (angstroms)
Figure 5. Plot of composite correlation functions as described in the text.
correlation function is well resolved at about 3.8 A, and a second peak appears near 5.9 A, indicating some water structure beyond the first solvation shell. Likewise, in Figure 3 the O,-O, correlation function shows an obvious peak due to the first hydration shell near 2.7 A and a weaker peak near 4.2 A due to a second solvation shell. On the tailing shoulder of this peak appears a rather sharp peak at approximately 5.0 A due to correlations of the anion oxygens with tightly held first solvation shell waters located on the oppositeside of thesulfate tetrahedron. A single configuration of the first solvation shell is shown in Figure 4. Roughly three water molecules are seen to be coordinated to each sulfate oxygen. This is more precisely quantified below. The intermediate range correlations tend to cancel when the atom pair contributions are combined in order to compare with a hypothetical experimental correlation function in which counterions are absent, as shown in Figure 5. This combined correlation function can be qualitatively compared to the correlation function shown in Figure 2 of the paper by Caminiti ef al.32 The two correlation functions shown in Figure 5 correspond to the combination of the pair correlation functions
The Journal of Physical Chemistry, VoL 98, No. 24, 1994 6229
Sulfate Anion in Water
,:I
S-H
I
I
0.00.0
o . o ' ~ ' ' ' ~ ' ~ ' " " " ' " ' " ' " ' ' ' " ' ' " " ' " ' " 0.0
1.0
2.0
3.0
4.0
5.0 6.0 r (angstroms)
7.0
8.0
9.0
I
10.0
20.0
40.0
80.0
80.0
1w.o
Tim (F4
Figure 6. Pair correlation functions for S-H and 0,-H pairs.
Figure 7. Decay of initial solvationshell as a function of time, normalized to the number of waters in the first solvation shell.
for S-Ow, four O,-Ow, and a first solvation shell water Ow-Ow pair, the last pair correlation being included with statistical weights of 14 and 7 in order to roughly mimic the concentrated solution of the experiment. Notable here is the reproduction of the shape of experimental correlation function a t distances greater than 4 A. Yamaguchi and Lindqvist33 also noticed structuring effects at 4.8 A and out to 6.0-7.2 A in their diffraction study. The diminished effect of the intermediatecorrelations in the weighted sum may explain why some authors believe that SO4" interacts relatively weakly and nonspecifically with water when compared to halogen ions.3S,37Unlike most anions, S042-is thought to be a structure-making ion with respect to the solvent around it1and is generally assumed to interact strongly with water molecules based on obvious energetic considerations and the observed difficulty in dehydrating sulfate crystals.3 The ability of to create structure in the second hydration shell is most likely due to the large charge separation and its tetrahedral geometry. Correlations for S-H pairs and 0,-H pairs are shown Figure 6. Again, intermediate range correlations appear around 4-5 A from the sulfate oxygen and even as far as 6.5 A from the sulfate sulfur. Here, it is somewhat more difficult to interpret the correlation function as far as assigning which hydrogen configurations result in which peaks. The first peak in the OpH correlation function obviously is due to first solvation shell waters in which the bond dipoles point at nearby sulfate oxygens. By integration, we find there are 3.3 hydrogens within 2.4 A of each sulfate oxygen and 6.7 hydrogens within 3.5 A, indicating that the waters are organized around sulfate in a manner analogous to Figure 1b,c and not in the bidentate binding mode of Figure la. Thisisalsoapparent in theinstantaneousconfigurationshown in Figure 4. The number of waters in the first hydration shell of S042- has been reported to be 6.0-14.3 based on analysis of diffraction data.3 All solutions studied by diffraction are relatively concentrated with the exception of one ZnSO, s t ~ d y , ~in6 which various concentrations of ZnS04 were used. Although the precise determination of solvation model parameters for analysis of the diffraction data is difficult for dilute solutions, the number of waters in the first solvation shell of S042-in 1:lOO ZnS04:H20 solutions was determined to he 12.4.36 Carrying out integration of the S-Ow correlation function to the first inflection point at 4.3 A results in 13.25 waters on average in the first solvation shell, and as indicated above there are 3.3 waters in close proximity to each sulfate oxygen. Detailed analysis of individual configurations obtained each picosecond over a 100-ps simulation showed that anywhere between 12 and 15 waters can be present in the first solvation shell of the anion. The number of waters found in the first solvation shell in concentrated solutions in diffraction studies is also dependent on the refinement model~sed.3~ In the initial model used by Caminiti
et to analyze diffraction data of (NH4)2S04 solutions, water molecules were assumed to interact weakly with the anion. For convenience, it was also assumed that NH4+ and water have identical structural properties and that a solution of (NH4)2S04 can be treated as a solution of water and S042-. With these assumptions, a value of 7.6 was obtained32for the number of waters in the first solvation shell of Sod2-. Ina subsequent report,38 the diffraction data reanalyzed without the assumption of structural indistinguishability between NH4+ and water were presented. The conclusion reached in the latter analysis was that two models were consistent with the diffraction data of (NH4)2SO4 in solution: the initial model in which ammonium ion is assumed to be structurally indistinguishable from water and the sulfate has about eight waters in the first hydration shell; a second model in which ammonium is not assumed to be structurally identical with water and S042- has 11 f 3 waters in the first hydration shell. It was concluded that the diffraction studies alone could not discriminate between the models. Other studies of concentrated solutions report a range of 6.0-14.3 waters in the first solvation shell (for a complete list see ref 3). The wide range of values exemplifies the difficulties involved in building solvation models for polyatomic ions, and it would appear that resolution of intermediate range correlations, especially those evident in the individual pair correlations, is critical for developing a sufficient model. The value of the self-diffusion constant obtained from the simulation of 1.5 X 10-5 cm2 s-1 is in the range of results found for other anions by molecular dynamics.2O The relative partial molar volume, AVm,was obtained from the difference between the partial molar volume of water and the partial molar volume of an infinitely dilute solution of S042-as reported by Marcus.1 The simulation value of 0.1 cm3 mol-' is an approximation obtained from the volume difference between the solution used here and a box of 895 SPC/E waters from a molecular dynamics simulation of 20 ps using the same simulation parameters as in the solution. The value obtained here is smaller than that obtained from experimental estimates (see Table 5). The calculated residence time for first solvation shell waters is 23 ps. Exponential decay of the autocorrelation function can be seen in Figure 7. The residence time for waters in the solvation shell of a bulk water has been calculated to be 4.5 ps39 for the MCY water model. For comparison, residence times for F- and C1- have been reported by Impey et ~ 1to be . 20.3 ~ ~and 4.5 ps, respectively. Generally, as the size of the ion increases, the residence time becomes smaller. However, this trend is not followed by S042-due to the magnitude of the charge and the significant charge distribution. The implications of the 23-ps residence time for the thermodynamic integrations is that the simulation time in going from the uncharged SO, to the fully charged final S042- model is much longer than the relaxation
6230 The Journal of Physical Chemistry, Vol. 98, No. 24, 1994 timeof the water environment. This impliesthat the requirement of reversibility in eq 2 should be satisfactorily met since the simulation time for charging and distributing the charges was 288 ps. Conclusion Development of a molecular mechanics potential of sulfate anion has resulted in a model that accurately reproduces relevant experimental solution data from various sources. The model developed here for sulfate shows a large negative free energy of hydration simply due to the solvation of the -2 charge and a significant positivefree energy for distributing the charge around the tetrahedral anion. The anion binds tightly to approximately 13 watersin thefirstsolvationshellinwhichthewatersareoriented with one hydrogen pointing toward a sulfate oxygen and the other hydrogen pointing toward neighboring water oxygens. The solvationshell is relatively long-lived with a characteristic lifetime for each water of 23 ps compared to the characteristic lifetime of the solvation shell of a bulk water of about 4.5 ps. The model suggested here should be useful for a number of further computational studies including solutions involving more complicated solutes where the OPLS potential parameter set is appropriate. Acknowledgment. W.R.C. thanks T. P. Straatsma and J. M. Briggs for useful discussions. This work was supported in part by grants from NSF, the Robert A. Welch Foundation, the Exxon Education Foundation, and the Metacenter Allocation Program of the NSF Supercomputer Centers. W.R.C. is supported in part by the NIH Houston Area Molecular Biophysics Training Program. References and Notes (1) (2) (3) (4) 2503. (5) (6)
Marcus, Y. Ion Solvation; Wiley: New York, 1985. Marcus, Y. Chem. Rev. 1988,88, 1475. Ohtaki, H.; Radnai, T. Chem. Rev. 1993, 93, 1157. Marlow, G. E.; Perkyns, J. S.;Pettitt, B. M. Chem. Rev. 1993.93, Hofmeister, F. Arch. Exp. Pothol. Pharmakol. 1888, 24, 247. Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553.
Cannon et al. (7) Pauling, L. Nature of the Chemical Bond, Cornell University Press: Ithaca. NY. 1966. (8) Chirlian, L.; Francl, M. J . Comput. Chem. 1987, 8, 894. (9) Breneman, C.; Wiberg, K.J. Comput. Chem. 1990, 11, 361. (10) Frisch, M.; Head-Gordon, M.; Trucks, G.; Foresman, J.; Schlegel, H.; Raghavachari, K.;Robb, M.; Binkley, J.; Gonzalez, C.; Defrees, D.; Fox, D.; Whiteside, R.; Seeger, R.;Melius, C.; Baker, J.; Martin, R.; Kahn, L.; Stewart, J.; Topiol, S.;Pople, J. Gaussian, Inc., Pittsburgh, PA, 1990. (11) Jorgensen, W. L. J. Phys. Chem. 1986,90,6379. (12) Jorgensen, W. L.;Briggs, J. M. J. Am. Chem.Soc. 1989,111,4190. (13) Straatsma, T. P.; McCammon, J. A. Annu. Rev. Phys. Chem. 1992, 43, 407. (14) Straatsma,T. P.; McCammon, J. A. J . Chem. Phys. 1991,95,1175. (15) Tembe, B. L.; McCammon, J. A. Comput. Chem. 1984,8, 281. (16) Straatsma, T. P.; McCammon, J. A. J. Comput. Chem 1990,11,943. (17) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R.J. Chem. Phys. 1984,81, 3684. (18) Straatsma, T. P.; Berendsen, H. J. C.; Postma, J. P. M. J . Chem. Phys. 1986,85,6720. (19) Straatsma, T. P.; Zacharias, M.; McCammon, J. A. Chem. Phys. Lett. 1992, 196, 297. (20) Impey, R. W.; Madden, P. A.; McDonald, I. R. J. Phys. Chem. 1983, 87. 5071. (21) Brunne, R. M.; Liepinsh,E.;Otting,G.; WBthrich,K.;vanGunsteren, W. F. J. Mol. Biol. 1993, 231, 1040. (22) Gianturco, F. A.; Coulson, C. A. Mol. Phys. 1968, 14, 223. (23) Marcus, Y. J . Chem. Soc., Foradoy Trans. 1 1987,83, 339. (24) Jenkins, H. D. B. Mol. Phys. 1975, 30, 1843. (25) Ladd, M. F. C.; Lee, W. H. J. Inorg. Nucl. Chem. 1968,30, 330. (26) Ladd, M. F. C.; Lee, W. H.J. Inorg. Nucl. Chem. 1961, 21,216. (27) Conway, B. E. J. Solution Chem. 1978, 7 , 721. (28) Ben-Naim,A. Solwrion Thermodynamics; Plenum: New York, 1987. (29) Straatsma, T. P.; Berendsen, H. J. C. J . Chem. Phys. 1988,89,5876. (30) Madura, J. D.; Pettitt, B. M. Chem. Phys. Lett. 1988, 150, 105. (31) Straatsma, T. P.; McCammon, J. A. Chem. Phys. Lett. 1991, 177, 433. (32) Caminiti, R.; Paschina, G.; Pinna, G. Chem. Phys. Lett. 1979, 64, 391. (33) Yamaguchi, T.; Lindqvist, 0. Acto Chem. Scand. A 1982,36, 377. (34) Caminiti, R. Chem. Phys. Lett. 1982, 86, 214. (35) Caminiti, R. Chem. Phys. Lett. 1982, 88, 103. (36) Licheri, G.; Paschina, G.; Piccaluga, G.; Pinna, G. 2.Noturforsch. 1982,37A, 1205. (37) Magini, M.; Lichen, G.;Paschina, G.; Piccaluga, G.; Pinna, G. %Roy Dffraction of Ions in Aqueous Medio: Hydrotion and Complex Formation; CRC Press: Boca Raton, FL, 1988. (38) Musinu, A.; Paschina, G.; Piccaluga, G. Chem. Phys. Lett. 1981,80, 163. (39) Impey, R. W.; Madden, P. A.; McDonald, I. R. Mol. Phys. 1982,46, 513.
