J. Phys. Chem. B 2004, 108, 17583-17590
17583
A Theoretical Investigation of the Shape and Hydration Properties of Aqueous Urea: Evidence for Nonplanar Urea Geometry Tateki Ishida† and Peter J. Rossky* Institute for Theoretical Chemistry, Department of Chemistry and Biochemistry, UniVersity of Texas at Austin, Austin, Texas 78712
Edward W. Castner, Jr.* Department of Chemistry and Chemical Biology, Rutgers, The State UniVersity of New Jersey, 610 Taylor Road, Piscataway, New Jersey 08854 ReceiVed: June 19, 2004; In Final Form: August 17, 2004
To model the structure of a urea molecule in aqueous solution including adaption of the solute electronic structure, we have used the reference interaction site model-self-consistent-field (RISM-SCF) method, describing the solute electronic structure at the ab initio level and hydration via an integral equation for the solvent distribution. The RISM-SCF model for aqueous urea gives a clearly nonplanar urea structure, with more than seven waters located within a contact distance defined for hydrogen bonding. The pyramidalization at the urea nitrogen sites is reduced in water relative to the gas phase, but the structure is closer to the gas phase structure than to the planar crystal structure. It is further shown that the solvent produces substantial polarization of the urea solute, with the dipole moment increasing from 4 to 7 D as the geometry and electronic structure are optimized to the aqueous environment. The overall result is to favor a high density of hydrogen bonds between urea and the surrounding solvent. The geometry obtained is intermediate between the previously known nonplanar gas phase structure and well-characterized planar crystal structure. An uncommon hemispherical shape for the electrostatic potential of around the carbonyl of urea is found, which affects the solvation structure.
I. Introduction Despite continuing research interest, a substantial controversy continues about the structure, hydration, and hydrogen-bonding properties of aqueous urea solutions.1,2 Because of the great interest in accurate simulations of protein folding and unfolding in aqueous urea solutions, an understanding of the electronic structure and hydrogen-bonding properties of dilute aqueous urea may be of use in developing more accurate potentials for molecular dynamics simulations. In this paper, we present theoretical work using a combined quantum/classical method to describe the structure and hydrogen bonding of dilute aqueous urea solutions. Using the reference interaction site model-selfconsistent-field (RISM-SCF) method,3-5 we obtain a nonplanar urea structure. Solvent radial distribution functions from these calculations then allow us to evaluate propensities for hydrogen bonding and to predict a greater degree of hydrogen bonding than is found in most other simulations that have used planar urea. Research on the properties of urea spans two centuries: Wo¨hler synthesized urea from the inorganic compounds silver cyanate and ammonium chloride in 18286 and in doing so led to the demise of the vitalism theory, the doctrine that organic compounds could only be produced from living organisms. Urea has great industrial importance, especially in the synthesis of fertilizers and plastics. Concentrated aqueous solutions display a remarkable tendency to render hydrocarbons soluble.7-9 † Present address: Department of Chemistry and Chemical Biology, Rutgers, The State University of New Jersey. * Corresponding authors.
Perhaps the most intriguing property of concentrated aqueous urea solutions is that these solutions can reversibly denature most proteins and some nucleic acids.1 The challenges in understanding urea in aqueous solution are many. A primary difficulty is that several theoretical descriptions of the liquid state of urea have been based on the planar crystal structure, which we show is not a valid starting point for describing urea solutions at infinite dilution. We remove this restriction and use a mixed quantum-classical description of the urea molecule to test its solution structure and hydrogenbonding interactions with the aqueous solvent. In aqueous solutions, urea has strong dipolar and hydrogenbond interactions with both water and itself. Aqueous urea is not thermally stable and reacts via several decomposition mechanisms, adding substantial complexity to experiments on solutions of aqueous urea. Because measurements of the solution properties of urea across a wide variety of physical and chemical environments are challenging, a number of controversies have arisen in the description of the properties and behaviors of urea in water. The molecular structure of urea depends on its environment; the molecular geometries for gas phase and crystalline urea are different. Three relevant conformers are shown in Scheme 1. The lowest energy structure for an isolated gas phase urea molecule is the nonplanar anti conformation, with a C2 point group symmetry.10 Crystal packing forces and geometries cause urea P4h21m to adopt the space group (#113) with local C2V molecular symmetry.11 It is important to note that crystalline urea has eight hydrogen bonds distributed between six neighbor-
10.1021/jp0473218 CCC: $27.50 © 2004 American Chemical Society Published on Web 10/14/2004
17584 J. Phys. Chem. B, Vol. 108, No. 45, 2004 SCHEME 1
ing molecules, with the carbonyl oxygen accepting four hydrogen bonds instead of the two that are usually expected for carbonyls. The fundamental difference between the gas phase and crystalline urea structures is that the crystalline urea has planar nitrogens, i.e., sp2 hybridization for the nitrogen. The gas-phase structure has nitrogens with bond angles closer to tetrahedral, or sp3 hybridization. It is important to note that the N-(Cd O)-N heavy atom framework is planar for both crystal and gas urea structuressonly the urea hydrogens can be nonplanar. Some of the first evidence that urea was nonplanar in the gas phase came from the microwave rotational spectrum, measured by Brown, Godfrey, and Storey in 1975.12 Their analysis of the urea rotational spectrum clearly showed two identical nitrogens and a small but nonzero inertial defect of -0.43 amu Å2, indicating that a small fraction of the urea mass (corresponding to the protons) is out of plane. The gas phase electronic structure calculations of Godfrey, Brown, and Hunter demonstrated that the nonplanar anti (C2) and syn (Cs) urea geometries are lower in energy, with the anti being the most stable.10 Because the X-ray and neutron diffraction data for crystalline urea are of very high quality, and all structures consistently obtain the same planar urea geometry, and because the crystalline environment permits strong intermolecular hydrogen bonding, theoretical chemists have sensibly begun their studies of aqueous urea solutions assuming a planar geometry.13-19 The purpose of this article is to consider the evidence for a planar or nonplanar structure of urea in water; a nonplanar urea geometry could be closer to the gas phase structure than that of the solid. Although only the protons in urea have clearly different geometries when comparing gas vs crystal structures, the presence or absence of planarity in solution can make a substantial difference in the hydration number for urea. The nonplanar structure that we predict for aqueous urea (vide infra) could form up to eight strong hydrogen bonds per molecule. This substantial propensity to form hydrogen bonds is potentially significant for a full understanding of the mechanism of action for aqueous urea in protein denaturation. Nevertheless, we note that the work presented here relates to the structure of urea at infinite dilution in aqueous solution, while the concentrations of aqueous urea used for protein (un)folding studies are very high, often above 8.0 M. For such nearly saturated aqueous urea solutions at room temperature, the mole ratio is less than 4 waters per urea molecule, very far from the dilute limit. In this paper, we discuss a series of electronic structure calculations on gas phase urea and on an aqueous solvent model. We have gained insight into the properties of both gas phase and aqueous urea that have not been previously described. We provide evidence in this article to support the following hypotheses: (1) Urea in aqueous solution is not planar, as in the crystal phase, but rather has NH2 groups with pyramidal geometries, as in the gas phase. (2) The electronic structure of urea in aqueous solution permits the urea carbonyl oxygen to accept not just two, but up to four hydrogen bonds, in analogy with the crystal structure. A primary reason for this unusual behavior is that the electrostatic potential (ESP) for urea is
Ishida et al. hemispherically symmetric about the oxygen atom. (3) In dilute aqueous solution, the average number of hydrogen bonds between water and urea molecules is found to be more than 7. This value is general accord with prior simulation results. Our discussion of the nonplanarity of aqueous urea will begin with several different electronic structure calculations. Our vibrational calculations show that the anti and syn geometries are connected by the planar C2V symmetry saddle surface; the two imaginary frequencies obtained for this transition structure prove that it cannot be the stable ground state. Qualitatively, no changes in geometry or ordering of ground-state energies were seen for calculations at the semiempirical, Hartree-Fock, or MP2 levels of theory. With evidence to suggest that hydrated urea may be nonplanar, we then discuss investigations of the solution structure, using a more complete electronic structure method, the reference interaction site model-self-consistent-field (RISM-SCF) model. Our RISM-SCF calculations include a high-quality electronic structure geometry optimization for urea (DZP basis set) with an applied field implied by a solvent model. We use the classical water model SPC potential. In section II. we outline the methods and theory used in both the gas and solution phase electronic structure calculations. Results are presented in section III. A discussion of the results, and comparison with previous results follows in section IV, with conclusions summarized in section V. II. Methods and Theory Gas Phase Urea at the HF and MP2/6-311++G(d,p) Level. For the study of gas phase urea, the Gaussian 98 program suite (Rev. A.7)20 was used with several levels of theory, up to and including MP2, and with basis sets up to the split triple-ζ with polarization and diffuse functions, 6-311++G(d,p). On increasing the level of theory and size of basis set, the observed geometries, dipole moments, and charges obtained from fits to the electrostatic potential were qualitatively unchanged. The urea geometry was optimized during the electronic structure calculations by imposing molecular symmetry to the conformers shown in Scheme 1. The lowest energy achieved during unconstrained geometry optimization was identical with the constrained C2 symmetry conformer. The RISM-SCF Method Applied to Aqueous Urea. Spatially varying electric fields arising from the water solvent distribution affect the solute electronic structure. Thus, we expect the urea geometry in aqueous solution to differ from either the gas or crystalline structures. Therefore, the description of urea in aqueous solution requires a method incorporating solvent effects and the solute electronic structure in a coupled manner. To model aqueous urea, we have applied the RISM-SCF approach.3-5 In prior work, RISM-SCF has succeeded in describing solvent effects microscopically. In the RISM-SCF description of the aqueous urea solution, interactions between the classical solvent and the quantum solute are included selfconsistently.3-5,21-23 In the present study, we apply the RISMSCF method to investigate the solute urea geometry solvated by the SPC water model, where the optimized structure is obtained from a variational calculation that minimizes the free energy. Details of the RISM-SCF method have been previously published.3-5 Therefore, we only briefly summarize the theoretical methods applied here. The solvation free energy of the system is defined in the RISM-SCF theory as
G ) Esolute + ∆µsol
(1)
Shape and Hydration Properties of Aqueous Urea
J. Phys. Chem. B, Vol. 108, No. 45, 2004 17585
where Esolute is the solute electronic energy in the solvated state and ∆µsol represents the excess chemical potential in solvation process. The solute wave function and solute-solvent correlation functions are calculated under variational conditions. In particular, the correlation functions are estimated by the extended RISM (XRISM) equations24-26 (with the hypernettedchain (HNC) closure relation):
huv ) ωu/cuv/ωv + Fωu/cuv/hvv
(2)
where the asterisks denote the spatial convoluting integral and h and c are matrices for the total and direct correlation functions, respectively. ω represents the intramolecular correlation functions, F is the solvent density, and the subscripts u and v denote solute and solvent sites, respectively. The RISM-SCF method provides a means of creating a solvent optimized molecular orbital set that satisfies the orthonormalization constraint for the optimal molecular orbital set {φi}
∑i 〈δφi|Fi - fi∑λ Vλbλ|φi〉 ) 0
(3)
where Fi is the Fock operator for the unperturbed solute, fi is the occupation number for orbital φi, Vλ is the solvent-induced electrostatic potential, and bλ is the population operator for the site λ. The solvent-induced electrostatic potential Vλ is given by
Vλ ) F
∞ ∑R ∫0 4πr2
qR
hλR(r) r
dr
(4)
where F is the solvent density, the qR are the partial charges are obtained by summing over the contributions from each nucleus and all electrons, and hλR(r) is the site-site correlation function. A useful diagram explaining the connections between the various steps used in the application of the RISM-SCF method is given by Ten-no, Hirata, and Kato in Figure 1 in their 1994 article.4 Several features of the RISM-SCF method clearly distinguish it from other quantum/classical (QM/MM) methods that incorporate a limited number of explicit solvent molecules via combined quantum/classical methods. The RISM-SCF method does not provide any specific snapshots of instantaneous solvent configurations. Rather, all of the structural information about the solute-solvent and solvent-solvent interactions is contained in the pair correlation functions, which are selfconsistently averaged via the integral equation treatment. The accuracy of the RISM-SCF solvation method depends on the quality of the classical force field for the solvent model. Thus, the method provides self-consistent solvent radial distribution functions about the solute and does not require convergence of the model over increasing numbers of solvent molecules as would be the case for a QM/MM method. However, this strength of RISM-SCF is also a limitation, as treatment of multiple solutes with specifically interacting solvents becomes a very difficult problem; i.e., the method applies only for a solute at infinite dilution. After convergence of the RISM-SCF calculation, the excess chemical potential is then given by the Singer-Chandler equation,27 modified with the hypernetted closure relation by Zichi and Rossky4,28
∆µu ) 2πFkT
∑u ∑V ∫0 (huV2 - 2cuV - huVcuV)r2 dr ∞
(5)
Figure 1. Estimated electrostatic potentials for the (gas phase) urea molecule, using MP2/6-311++G(d,p) optimizations. The envelope shows the contour for a density of 0.04 electron/Å3. From top to bottom: (a) C2, (b) C2V, and (c) Cs urea conformers are shown. Note that the sp2 hybridization for the planar C2V geometry eliminates the lone pair density on the nitrogen atoms found for the other two structures.
RISM-SCF Computational Details. The SCF ab initio calculations were performed at the restricted Hartree-Hock (RHF) level in both the gas and solution phase. The Dunning double-ζ plus polarization (DZP) basis set29 was employed. All calculations were carried out under C1 symmetry condition (no symmetry). For the RISM-SCF model for aqueous urea, the analytical energy gradient technique was applied to optimize the solute geometry.5 The partial charges on the solute atomic sites were obtained by an optimization procedure that reproduces the electrostatic potential. The interaction between solvent molecules or between solute and solvent is represented by the pairwise Lennard-Jones plus Coulomb potential, which also accounts for hydrogen-bonding effects. All remaining parameters for both the solute and solvent molecule used for the RISM calculation are adopted from the literature on computer simulation14,15 and are given in Table 1. The simple point charge (SPC)30-like model was used for water solvent; Lennard-Jones parameters for the hydrogen atom (σ ) 1.0 Å and ) 0.056 kcal-mol-1) are included.24 The solvent density of F ) 0.033 34 Å-3 was used. All calculations were performed with T ) 298.15 K. III. Results Gas-Phase Urea by MP2/6-311++G(d,p) and SCF/DZP. We have obtained quantitative agreement with the previous results of Godfrey, Brown, and Hunter10 for the three structures
17586 J. Phys. Chem. B, Vol. 108, No. 45, 2004
Ishida et al.
TABLE 1: Parameters for the Solute and Solvent Molecule atom site
σ (Å)
-1
(kcal mol )
-
q (e )
TABLE 2: Gas Phase Dipole Moments for the Three Urea Conformers
Solute C O N H
3.750 2.960 3.240 1.780
O H
3.166 1.000
µ (D)
0.105 0.210 0.170 0.020 Solvent 0.155 0.056
C2
C2V
Cs
3.46
4.33
4.24
TABLE 3: Summary of Optimized Geometric Parameters for Urea from MP2/6-311++G(d,p) and RISM-SCF/DZP Models for Urea -0.82 0.41
of ground-state urea at the MP2/6-311++G(d,p) level of theory. The C2V conformer found in the crystalline phase is not a stable geometry in the gas phase; normal-mode analysis provides two imaginary frequencies. A saddle surface connects the C2 and Cs symmetries by inversion of the amino protons. In agreement with the Godfrey, Brown, and Hunter analysis, we find that the C2 structure is about 1.2 kcal/mol lower than the Cs structure. It is likely that the C2 and Cs structures do not interconvert by the local maximum energy planar C2V structure but are joined via a nonplanar large-amplitude vibrational motion with a barrier of 0.38 kcal/mol relative to the Cs structure.10 They also showed that the zero-point vibrational energy of the C2 conformation is above the barrier connecting the C2 and Cs geometries. Hence, the appropriate nuclear vibrational state may be delocalized over multiple minima. Masunov and Dannenberg have systematically investigated further details of urea using high-level ab initio calculations for a variety of methods and basis sets.31 They also observed that the zero-point vibrational energy might be above the barrier to planarity. They reported that the path for conversion between the C2 and Cs conversion could be attributed to the pyramidalizing motion of the NH2 groups. Thus, there are multiple previous reports that urea is easily distorted in the gas phase. Here, we should note that the recent simulation result by Kallies32 indicated urea molecule in water solution possesses a nonplanar, i.e., unsymmetrical, structure. It is important to emphasize that the nature of unsymmetrical structures in solution may be rather different from unsymmetrical structures in the gas phase. In the solution phase, the solute is strongly interacting with the solvent and is, additionally, polarized by the solvent. If the solvent does not follow the motion, the large-amplitude vibration implicated in the gas phase for the interconversion of C2 and Cs conformers will entail a large energetic barrier associated with hydrogen bond loss. Alternatively, if the solvent does follow the solute rearrangement, the effective mass associated with the (cooperative) coordinate will be large and entail a correspondingly small zeropoint energy. In either case, in contrast to the possibility in the gas phase, the result in solution will not be well described by a state in which the amino protons are delocalized among multiple minima. The solvent-solute interaction localizes the solute into one minimum at a time, with conventional classical barrier crossing dynamics expected between them. The urea electrostatic potentials for the MP2/6-311++G(d,p) calculations are shown in Figure 1 for the C2, C2V, and Cs conformers, from top to bottom, respectively. Note that the shape of the field surrounding the urea carbonyl oxygen atom is distinctly hemispherical in shape for all of the three conformers. Gas phase dipole moments for each of the three conformers are listed in Table 2. The influence of the electron density from the nitrogen lone pair electrons is clearly evident for the nonplanar C2 and Cs conformers but is absent for the planar C2V geometry because of the conjugation to the carbonyl group. We note that the charges obtained from the fit to the electrostatic potential are robust, in the sense that the 3-dimensional
model R(C1-O2), Å R(C1-N3), Å R(C1-N4), Å A(N3-C1-N4), deg D(O2-C1-N4-H6), deg
RISM-SCF/DZP (no symmetry)
MP2/ 6-311++G(d,p) anti C2 syn Cs (gas) (gas)
anti C2 (gas) (state I)
anti C2 (water) (state III)
1.217 1.390 1.390 113.06 147.40
1.200 1.375 1.375 114.15 151.66
1.235 1.345 1.345 117.55 166.13
1.218 1.386 1.386 114.15 157.67
hemispherical shape of the urea ESP is recovered when the fitted charges based on the 2-dimensional C2V conformer are used to calculate the ESP in the reverse direction. Figure 2a,b provides the atomic numbering scheme used and the optimized ground-state geometries for urea in the gas phase and in water, while Table 3 presents a comparison of several important optimized geometric parameters for both the gas phase and hydrated structures. In the gas phase, the bond length and bending angle calculated are in good accord with experimental values by Godfrey et al.10 As seen in Figure 2a, the SCF/DZP level calculation shows that the urea molecule adopts the anti nonplanar geometry in agreement with the MP2 results. As with experimental X-ray and neutron structures, the heavy atoms (C1, O2, N3, and N4) remain planar with a null dihedral angle though no constraint was used during geometry optimization. The values for the MP2/6-311++G(d,p) bond lengths are about 1% longer than the SCF/DZP values, and the N3-C1-N4 bond angles agree to within 1°. It should be noted here that Godfrey et al. assumed C2 symmetry for the molecular shape in analyzing experimental data. Geometry from RISM-SCF/DZP: Aqueous Model. For urea in water modeled using the RISM-SCF/DZP method, the solute geometry is altered from that in the gas phase. Figure 2b shows that the minimum-energy urea structure in water preserves the gas phase anti conformation. The heavy atom plane defined by C1-O2-N3-N4 is preserved in water. The C1-O2 bond length increases compared to that in the gas phase, arising from electron density polarization due to hydrogen bonds at the carbonyl oxygen. This enhanced C1-O2 bond length enlarges the urea dipole moment, further increasing the magnitude of the solvation energy. The carbonyl bond has less double-bond character, and the electron density spreads to the C-N bonds. As a result, the C-N bond is strengthened and the bond length decreases in aqueous solution. Correspondingly, the N3-C1N4 angle increases in the aqueous phase. The data in Table 3 show that the urea dihedral angle O2C1-N4-H6 increases significantly on going from the gas to aqueous phase. The deviation from planarity for the amino groups is reduced by about 5° in aqueous solution, with the amino groups less pyramidal in water than in a vacuum. From the present results, it is clear however that the urea solute is still nonplanar, but less so in the aqueous than in the gas phase. Masunov et al.31 reported that for geometries optimized at the HF level the urea molecule favored a planar form in increasing uniform electric fields but that optimization at the MollerPlesset second-order (MP2) perturbation level predicted a nonplanar form with increasing electric fields. In the present
Shape and Hydration Properties of Aqueous Urea
J. Phys. Chem. B, Vol. 108, No. 45, 2004 17587
SCHEME 2
study, the nonplanar form of urea persists in water at the HF level, even in the molecular solvent field. Solvation Process and Solute Properties. In Scheme 2, the energy correlation diagram for the RISM-SCF treatment of urea
Figure 2. RISM-SCF/DZP molecular geometry and site assignment for (a) gas phase urea and (b) urea in dilute aqueous solution.
is shown schematically along with our calculated results. State I represents the gas-phase urea anti conformer geometry optimized at the SCF/DZP level, having a dipole moment of
TABLE 4: Solute Site Charges from RISM-SCF Calculations (Units Are Electronic Charges) gas
water
atom site
I
IV
III
II
C1 O2 N3 N4 H5 H6 H7 H8
1.07 -0.68 -1.02 -1.01 0.40 0.40 0.40 0.44
1.08 -0.71 -1.01 -1.00 0.41 0.40 0.42 0.41
1.31 -0.96 -1.14 -1.13 0.47 0.49 0.47 0.49
1.27 -0.87 -1.15 -1.13 0.46 0.48 0.46 0.48
4.03 D. State II represents the urea molecule transferred to the water at the gas phase geometry of state I. The water interactions polarize the urea, increasing the dipole moment to 5.47 D. State III represents the urea molecule in water after reequilibration of the structure to accommodate the urea/water interactions. The relaxed structure is substantially more polarized, having a dipole moment of 6.91 D. At the equilibrium geometry in water, state III, the urea amino groups are still clearly nonplanar, with the value of the dihedral angle being halfway between the gas and crystalline phase conformers. To complete the thermodynamic cycle, we consider state IV, in which the urea molecule with the geometry optimized in water is placed into the gas phase. Removal of the water decreases the dipole moment to 4.88 D. The accuracy of the RISM-SCF/DZP computed dipole moment is worth considering. At state I of Scheme 2, the urea gas phase dipole moment is 4.03 D, which is in reasonable agreement with the experimental value of 3.846 D reported by Brown et al.12 The dipole moment value of 3.46 D from our all-electron MP2/6-311++G(d,p) calculation for the C2/anti conformer is 0.47 D lower than value of 3.93 D from the frozencore MP2 calculation. Previous theoretical work by A° strand et al. using a continuum solvation model reported the aqueous urea dipole moment to be 7.26 D,16 in good qualitative agreement with the RISM-SCF value. Experimental estimates from dielectric relaxation data range from 6.1 to 7.9 D.33 Thus, the RISM-SCF/DZP result agrees well with prior results, reinforcing the validity of other quantities. The charges on the solute sites obtained by fitting atomcentered charges to the ESP are summarized in Table 4. In transfer of urea from gas phase to water (state I to II), the solute charge on the O2 site becomes more negative, and the C1 site charge becomes more positive. The change of the charges on the amino groups found for transferring urea from gas to water is somewhat smaller than for the carbonyl group atoms, but these are also polarized because of interactions with the water.
17588 J. Phys. Chem. B, Vol. 108, No. 45, 2004
Figure 3. Radial distribution functions g(r) between the urea and water calculated from RISM-SCF results: (a) red: urea O2 and water proton, Hw sites; (black) urea O2 and water oxygen, Ow sites; green: urea N3 site and Ow sites; orange: urea C1 and Ow sites; (b) green: urea N3 and water Hw sites; blue: urea C1 and water Hw sites; (c) cyan: urea H5 and water Ow sites; magenta: urea H6 and water Ow sites.
In water, the O2 site charge becomes even more negative (-0.87 f -0.96) as the geometry is reoptimized in response to the electric field of the solvent. At the same time, the charge on the C1 site increases from 1.27 to 1.31, consistent with an increasing polarization of the carbonyl bond. The increase in the C1-O2 bond length in water also contributes to the enhancement of the urea dipole. Because of the increase of the urea dipole moment on transfer from gas phase to water, the difference in the excess chemical potential between state III and II becomes negative. At the same time, the polarized character of the solute electronic structure is enhanced successively in the solvation process toward state III, and thus the solute electronic energy increases as illustrated in Scheme 2. However, this increase in solute electronic energy is necessarily more than compensated by the excess chemical potential, and the total free energy at the equilibrium state III becomes lower than at state II by about 2 kcal/mol. The dipole moment value at state IV, with the equilibrium solvated structure, is larger than that at state I; in the absence of the stabilizing effect of the water solvent, the electronic structure at the state IV has a higher energy than at state I. Solvation Structure. In Figure 3a-c, computed results for solvent radial distribution functions (rdf’s) around urea atoms are shown. The g(O2-Hw) function shown in Figure 3a shows a sharp peak at about 1.70 Å that we assign to strong hydrogen bonding between the urea O2 site and the water hydrogen for the first aqueous solvation shell. The broader peak in the g(O2Ow) rdf peaked at 2.97 Å includes contributions both from
Ishida et al. waters directly H-bonded to the urea O2 site and possibly also from nearby waters that are H-bonded to those in the first coordination layer. In contrast, the urea crystal H-bonding distances from Swaminathan et al. are about ∼2.0 Å for amino proton to carbonyl oxygen and 2.99 Å for N-O distances between urea molecules.11 In Figure 3b, the rdf’s between water proton and the urea C1 and N3 atoms are shown. Two distinct peaks arise for each rdf. For the C1-Hw rdf, the height of the first peak is lower than that of the second peak because the first peak is assigned to van der Waals contact between the C1 and Hw atoms. The second peak can be attributed to the interaction of the C1 site with those water molecules that are also directly adjacent to the O2 site. This interpretation is reasonable considering that the urea C1-O2 bond length in water is 1.235 Å and that the first peak appears at around r ) 2 Å in the rdf for the urea O2-Hw. Both the C1 and Hw sites have positive charges, however, and therefore the position of the second peak is shifted to a longer distance than r ) 4 Å. We can also apply similar arguments to explain the behavior of the N3-Hw rdf. On the basis of the computed result of about 2.25 Å for the distance between the N3 and O2 atoms, the second peak is assigned to arise from the H atom distribution around the solute O2 site. However, since the N3 and the solvent H atoms have charges with opposite sign, the second peak position is found at r < 4 Å. Figure 3c displays the rdf’s for the urea H5 and H6 sites and the water Ow site. These nonequivalent rdf’s show similar features. In agreement with the neutron diffraction assignment from Kameda et al.,34 this is consistent with each proton donating a hydrogen bond to a water molecule at low urea concentration. On comparing the peaks of the H5-Ow and H6Ow rdf’s with crystal H-bonding distances, we find the RISMSCF results from Figure 3c to have maxima at 2.27 and 2.22 Å, whereas the crystal data of Swaminathan et al. show equivalent distances for H5 and a neighboring urea O2 to be at ∼2.0 Å.11 Using the radial distribution functions from the RISM-SCF calculations, we can investigate the configuration of solvent molecules around urea that contribute to hydrogen bonding. We estimate the number of water molecules that are within the cutoff distance for hydrogen bonding using the following formula with an appropriate upper limit r:
N ) 4πFγ
∫0rgRγ(r)r2 dr
(6)
where R and γ represent the solute and solvent sites, respectively, and F is the solvent density. For the rdf’s in Figure 3ac, the hydrogen-bond coordination numbers were evaluated for each case with the integration limit r set to the first minimum of the radial distribution function. We obtain N ) 3.16 from the g(O2-Hw), N ) 1.87 from the g(H5-Ow) and g(H7-Ow) pair, and N ) 2.30 from the g(H6-Ow) and g(H8-Ow) pair, for a sum of 7.33 waters coordinated to urea in dilute aqueous solution. The RISM-SCF rdf’s arise from specific atoms on water molecules, not entire water molecules. Thus, we must exercise caution and differentiate between the number of hydrogen bonds, which are cleanly defined by integration of the rdf’s, and estimates of the hydration number, which is complicated by the fact that individual water molecules may make two H-bonds to the urea. As discussed by several authors,16,18,32 one can envision a single water making a double hydrogen bond to urea in two places. First, consider each of the urea protons H6 and H8 donating H bonds to a single water oxygen Ow, and second,
Shape and Hydration Properties of Aqueous Urea consider a water oxygen accepting a H bond from N3H3- -Ow (or N4-H5- -Ow) with the same water molecule donating an H bond from Ow-Hw- -O2. Because of this possibility for counting doubly-H-bonded water molecules twice in the rdf’s, we must point out that our RISM-SCF calculations are consistent with a hydration number for urea at infinite dilution having a value ranging from 4.3 to 7.3. Considering entropic effects, the most probable value for the hydration number will likely be about six waters per urea. This estimated value for the urea hydrogen-bond number in excess of 7 is similar to the hydrogen-bond numbers for aqueous urea estimated from earlier simulation results. Our prediction of 3.16 waters coordinated to the urea oxygen atom O2 is twice what would be expected for other carbonyl, ether, alcohol, or water oxygens and is larger by about one water molecule than the values of 2.3 estimated by Duffy et al.17 and 2.5 by A° strand et al. in prior simulations.16 We attribute the larger than expected hydration number at the O2 urea oxygen to the nearly hemispherical shape of the electrostatic potential and large partial charge. Note that although this prediction of a hydration number >3 is substantially larger than the value of 1.7 obtained from the recent analysis of neutron diffraction experiments,35 these experiments at high concentration (1 urea:4 water ratio) cannot be compared directly to our calculations that correspond to dilute solutions. The only other solvent molecule that we are aware of that possesses such a hemispherical electrostatic potential and charge density on the oxygen is dimethyl sulfoxide (DMSO).36 The relatively high estimated number of hydrogen bonds of urea in dilute aqueous solution is one of the principal results of this paper. Previous literature references for the urea hydration numbers vary widely. Monte Carlo simulations gave an average of five water molecules coordinated to urea.17 Using their NEMO potential developed for polarizable potential MD simulations, A° strand et al. found 5.7 water molecules surrounding urea by integrating the first peaks of the rdf’s, in analogy with eq 4.16 In their work, we should note that the electronic structure calculations used to obtain the effective charges were constrained to the planar geometry, and the solvent distribution and the polarization effect of the solute electronic structure were not treated in a self-consistent manner.16 In the present study, the RISM-SCF method is fully self-consistent, and polarization effects enhance the magnitude of the solute charges. Therefore, the origin of the difference between our estimate for the number of hydrogen bonds and the literature values is likely the result of these differences in theoretical methods, in addition to differences in solute geometry. IV. Discussion Many theoretical approaches have been applied to the problem of the urea geometry in water. Molecular dynamics simulations with the solute-solvent interaction potential determined by ab initio calculations have been reported by A° strand et al.16,18 It was assumed that planar urea was the most probable conformation in aqueous solution, but it was noted that the amino groups may be somewhat flexible. Mennucci et al. have considered solvent effects during geometry optimizations by employing the polarized continuous model (PCM).37 They showed that geometry changes occur on going from vacuum to water solution for a planar urea structure. A substantial amount of work has been done using neutron diffraction to study the structure of aqueous urea solutions.34,35,38,39 By studying a range of relative isotopic dilutions of H/D and 14N/15N, the intermolecular partial radial distributions
J. Phys. Chem. B, Vol. 108, No. 45, 2004 17589 are obtained. The work of Kameda et al. concludes that for 15 mol % aqueous urea the hydration number for the NH2 groups is 4.1 or such that each amino proton donates a hydrogen bond to neighboring water molecule.34 At this high urea concentration, there are statistically 5.67 waters per urea in solution, so it is not surprising that the hydration number for the urea O atom is 2.0 (rather than the maximum value of 4) because there are not enough waters available for H bonding.34 Other recent work by Soper et al. shows that for 20 mol % aqueous urea 5.7 H bonds to neighbors are found in the analysis35 of the combined neutron diffraction/empirical potential structure refinement method.40,41 The most recent neutron diffraction results35 are again consistent with the N-(CdO)-N heavy atom framework being planar, with the hydrogen atoms slightly out of plane. Though we must use caution in predicting the nature of concentrated aqueous urea solutions based on models for dilute solutions, it seems that the nonplanarity (and hence amino group flexibility) of urea will certainly persist at higher concentrations. The enthalpies for reorientation of urea obtained from temperature-dependent Kerr dynamics experiments differ between 1.3 and 8.3 M aqueous solutions, indicating that different local structures must exist in the solutions.42 The hemispherical shape of the urea electrostatic potential that permits four H bonds in the crystal will still obtain for the case of concentrated aqueous urea, implying that one should expect aqueous solutions to have a urea carbonyl oxygen that accepts more than the traditionally accepted value of two H bonds. Taken together, this collection of results leads us to speculate that the following two principles will describe behaviors of concentrated aqueous urea solutions. (1) Urea molecules in concentrated aqueous solutions tend to maximize the density of hydrogen bonds per urea molecule, making up to about seven H bonds with neighboring water and urea molecules. (2) The unique properties of the urea electronic structure permit it to form H bonds with water over a greater range of angular flexibility relative to bulk tetrahedral water networks because of both the hemispherical shape of the ESP around the oxygen acceptor atom and the flexibility of the amino donor groups. A substantial body of work on simulations of aqueous urea solutions has been done since the early work by Kuharski and Rossky.8,13 Numerous models have been considered,9,16,18,19,32,43-51 with a number of innovative parametrizations and improvements of the potentials.16,18,32,44,46,47,49,50 Though several have relied on ever improving electronic structure calculations, even these have made parametrizations based on urea-water dimers, which will overemphasize multiple H bonds between urea and water. With the exception of the work by Tsai, Gerstein, and Levitt,44 the other simulations all considered only planar urea molecules. It will be interesting to learn whether further improvements in the parametrization of potentials for molecular simulations can be made with incorporation of a nonplanar urea model, and atom-centered charges that reproduce the unique hemispherical shape of the electrostatic potential around the urea carbonyl oxygen. V. Conclusions In the present study, we applied the RISM-SCF method to the case of urea aqueous solutions at infinite dilution. The solute optimized geometry, dipole moment, solute site charge, and the energy correlation diagram in the gas phase and solution were calculated. We summarize our results as follows: 1. The solute geometry of urea has a nonplanar form in the gas phase, associated with pyramidalization at the N sites.
17590 J. Phys. Chem. B, Vol. 108, No. 45, 2004 2. In water, urea adopts a nonplanar structure with the amino protons out of the central heavy atom plane formed by the C1, O2, and N3/N4 atoms. The shape of urea in dilute aqueous solution has a nitrogen pyramidalization between that of the gas phase structure and the planar crystal structure. 3. There is significant solvent-induced polarization of the urea molecule in water. The estimated value of the urea dipole moment in aqueous solution is 6.91 D and agrees reasonably with other theoretical estimates. The H-bonding of water to the urea O2 site draws electron density from the C1dO2 bond, leading to an increased bond length, increased charge on the urea oxygen atom O2, and an enhanced dipole moment relative to the gas phase model. 4. The excess chemical potential of solvation of -7.45 kcal/ mol for urea in water exceeds the unfavorable increase in electronic structure energy of 5.30 kcal/mol; the net solvation free energy is thus -2.15 kcal/mol in the RISM-SCF model. 5. The solvent distributions about urea are illustrated by the radial distribution functions that show a strong agreement between the estimated density of hydrogen bonds around the O2 urea site and the experimental results. Our estimation is that, on average, more than seven water molecules form H bonds with urea in dilute solution. On the basis of the results presented here, it seems reasonable to suggest that classical force fields for molecular dynamics simulations of aqueous urea should be reconsidered. In particular, the charges used in a classical potential need to reflect the substantial solvent-induced polarization of the urea dipole and should reproduce the unusual hemispherical symmetry of the ESP surrounding the urea oxygen O2 site. Also, incorporation of both flexible, pyramidalized nitrogens with nonplanar protons will increase the hydrogen-bond donating capacity of urea relative to a rigid planar model. New potentials incorporating these design criteria may provide substantial improvements in modeling accuracy in studies of both protein folding and solvation of hydrocarbons in concentrated aqueous urea solutions. Acknowledgments. P.J.R. acknowledges support from the R.A. Welch Foundation and Texas Advanced Research Program. E.W.C. acknowledges partial support from the donors of the Petroleum Research Fund and the NSF. References and Notes (1) Creighton, T. E. ProteinssStructures and Molecular Properties, 2nd ed.; W.H. Freeman and Co.: New York, 1993. (2) Tirado-Rives, J.; Orozco, M.; Jorgensen, W. L. Biochemistry 1997, 36, 7313-7329. (3) Ten-no, S.; Hirata, F.; Kato, S. Chem. Phys. Lett. 1993, 214, 391396. (4) Ten-no, S.; Hirata, F.; Kato, S. J. Chem. Phys. 1994, 100, 74437453. (5) Sato, H.; Hirata, F.; Kato, S. J. Chem. Phys. 1996, 105, 15461551. (6) Wo¨hler, F. Ann. Phys. Chem. (Leipzig) 1828, 88. (7) Frank, H. S.; Franks, F. J. Chem. Phys. 1968, 48, 4746-4757. (8) Kuharski, R. A.; Rossky, P. J. J. Am. Chem. Soc. 1984, 106, 57945800. (9) Wallqvist, A.; Covell, D. G.; Thirumalai, D. J. Am. Chem. Soc. 1998, 120, 427-428. (10) Godfrey, P. D.; Brown, R. D.; Hunter, A. N. J. Mol. Struct. 1997, 413, 405-414. (11) Swaminathan, S.; Craven, B. M.; McMullan, R. K. Acta Crystallogr. B 1984, 40, 300-308. (12) Brown, R. D.; Godfrey, P. D.; Storey, J. J. Mol. Spectrosc. 1975, 58, 445-450. (13) Kuharski, R. A.; Rossky, P. J. J. Am. Chem. Soc. 1984, 106, 57865793.
Ishida et al. (14) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G. S.; Profeta, J.; Weiner, P. J. Am. Chem. Soc. 1984, 106, 765. (15) Jorgensen, W. L.; Swenson, C. J. J. Am. Chem. Soc. 1985, 107, 569-78. (16) A° strand, P.-O.; Wallqvist, A.; Karlstro¨m, G.; Linse, P. J. Chem. Phys. 1991, 95, 8419-8429. (17) Duffy, E. M.; Severance, D. L.; Jorgensen, W. L. Isr. J. Chem. 1993, 33, 323. (18) A° strand, P. O.; Wallqvist, A.; Karlstro¨m, G. J. Chem. Phys. 1994, 100, 1262-1273. (19) Idrissi, A.; Sokolic, F.; Perera, A. J. Chem. Phys. 2000, 112, 94799488. (20) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, Rev. A.7. Revision A.7 ed.; Gaussian, Inc.: Pittsburgh, PA, 1998. (21) Ishida, T.; Hirata, F.; Sato, H.; Kato, S. J. Phys. Chem. B 1998, 102, 2045-2050. (22) Ishida, T.; Hirata, F.; Kato, S. J. Chem. Phys. 1999, 110, 39383945. (23) Ishida, T.; Rossky, P. J. J. Phys. Chem. A 2001, 105, 558-565. (24) Hirata, F.; Rossky, P. J. Chem. Phys. Lett. 1981, 83, 329-334. (25) Hirata, F.; Pettitt, B. M.; Rossky, P. J. J. Chem. Phys. 1982, 77, 509-520. (26) Hirata, F.; Rossky, P. J.; Pettitt, B. M. J. Chem. Phys. 1983, 78, 4133-44. (27) Singer, S.; Chandler, D. Mol. Phys. 1985, 55, 621-625. (28) Zichi, D. A.; Rossky, P. J. J. Chem. Phys. 1986, 84, 1712. (29) Dunning, T. H.; Hey, P. J. dzp-Modern Electronic Structure Theory; Plenum: New York, 1977. (30) Berendsen, H. J. C.; Postma, J. P. M.; van Gusteren, W. F.; Hermas, J. Intermolecular Forces. In Intermolecular Forces; Reidel: Dordrecht, 1981. (31) Masunov, A.; Dannenberg, J. J. J. Phys. Chem. A 1999, 103, 178184. (32) Kallies, B. Phys. Chem. Chem. Phys. 2002, 4, 86-95. (33) Pottel, R.; Adolph, D.; Kaatze, U. Ber. Bunsen-Ges. Phys. Chem. 1985, 79, 278. (34) Kameda, Y.; Naganuma, H.; Mochiduki, K.; Imano, M.; Usuki, T.; Uemura, O. Bull. Chem. Soc. Jpn. 2002, 75, 2579-2585. (35) Soper, A. K.; Castner, J.; E. W.; Luzar, A. Biophys. Chem. 2003, 105, 649-666. (36) Wiewio´r, P. P.; Shirota, H.; Castner, J. E. W. J. Chem. Phys. 2002, 116, 4643-4654. (37) Mennucci, B.; Cammi, R.; Cossi, M.; Tomasi, J. J. Mol. Struct.: THEOCHEM 1998, 426, 191. (38) Finney, J. L.; Soper, A. K.; Turner, J. J. Phys.: Condens. Matter 1989, 156-157, 151-3. (39) Turner, J.; Finney, J. L.; Soper, A. K. Z. Naturforsch., A 1991, 46, 73-83. (40) Soper, A. K. J. Mol. Liq. 1998, 78, 179-200. (41) Soper, A. K. Mol. Phys. 2001, 99, 1503-1516. (42) Shirota, H.; Wiewior, P. P.; Castner, E. W., Jr. Manuscript in revision. (43) A° strand, P.-O.; Wallqvist, A.; Karlstro¨m, G. J. Phys. Chem. 1994, 98, 8224-8233. (44) Tsai, J.; Gerstein, M.; Levitt, M. J. Chem. Phys. 1996, 104, 94179430. (45) Mountain, R. D.; Thirumalai, D. J. Am. Chem. Soc. 2003, 125, 1950-1957. (46) Mountain, R. D.; Thirumalai, D. J. Phys. Chem. B 2004, 108, 68266831. (47) Weerasinghe, S.; Smith, P. E. J. Phys. Chem. B 2003, 107, 38913898. (48) Weerasinghe, S.; Smith, P. E. J. Chem. Phys. 2003, 118, 59015910. (49) Sokolic, F.; Idrissi, A.; Perera, A. J. Chem. Phys. 2002, 116, 16361646. (50) Sokolic, F.; Idrissi, A.; Perera, A. J. Mol. Liq. 2002, 101, 81-87. (51) Idrissi, A.; Cinar, E.; Longelin, S.; Damay, P. J. Mol. Liq. 2004, 110, 201-208.
