J. Phys. Chem. 1994,98, 8224-8233
8224
Molecular Dynamics Simulations of 2 m Aqueous Urea Solutions P.-0.Astrand,*J A. Wallqvist,$9s and G. Karlstromt Theoretical Chemistry and Physical Chemistry 2, Lund University, P.O. Box 124, 221 00 Lund, Sweden Received: May 1I , 1994”
Results from molecular dynamics simulations of 2 m aqueous urea solutions are presented. A recently derived set of ab initio intermolecular potentials with explicit atomic polaritabilities is employed. In the simulated systems consisting of 10 urea molecules and 277 waters, many different urea-urea complexes are found, in contradiction to other recent studies where empirical force fields have been used. These differences are due to the inability of the empirical force fields to properly assign the global urea dimer energy minimum. Four simulations with different initial starting configurations were carried out in order to ascertain whether the systems are capable of describing the dynamics of urea dimerization processes. Since we found complex life times and orientational decay times that are longer than the simulation time (88 ps), it was not possible to achieve equilibrium with respect to the number and character of urea-urea complexes. However, comparison with neutron scattering experiments gives that the partial radial distribution functions Gp~(r), GHH(r), and GH(c,o,N)are insensitive to the degree of urea complexation and therefore also to the intermolecular potential adopted in a simulation. Future neutron scattering experiments are proposed, where the presence of urea-urea complexes can be determined with certainty.
1. Introduction
Aqueous urea solutions have interesting and exceptional properties;urea is highly soluble in water,’ increases the solubility of hydrocarbons in water,2 inhibits micellar aggregation,3 and denaturates proteins.435 Several, mutually exclusive, thermodynamic modeW9 have been proposed to describe these properties, and a number of molecular dynamics (MD) simulation^^^^^ of aqueous urea solutions have been aimed at investigating these models. In a model proposed by Frank and Franks (the FF model),9 urea is supposed to act as a “water structure breaker”, altering the structure of bulk water. This model is supported in several experiments such as an NMR experiment,20 an X-ray measuremenLZ1 calorimetric experiment~,2~.~~ and a recent Raman spectroscopy studyOz4On the other hand, in two recent MD simulationsl5.16 the bulkwater region has been found to be almost indistinguishable from pure water with regard to structural and dynamical properties. This result is also corroborated by earlier simulations of Kuharski and RosskylO and Tanaka et al.I2 but contradicted by Cristinziano et al.14 However, Boek and Brielsl6 investigated the small truncated octahedron used by Cristinziano et al.14and found artifacts arising from the small system size and periodicboundary conditions. Urea causes no apparent disruption of water structure according to an earlier calorimetricexperiment25 and recent neutron scattering experiment^.^^^^^ In another model, the SKSS model, first proposed by Schellman,6and later refined by Kresheck and Scheraga7 as well as by Stokes,8it is assumed that formation of urea dimers and oligomers plays an important role for explaining the behavior of aqueous urea solutions. That urea forms dimers and higher order aggregates is found in osmotic pressure measurements.28 The MD simulations are, however, not conclusive concerning the stability of urea dimers. In two simulations, 14.16 where the GROMOS p ~ t e n t i a lcombined ~ ~ ~ ~ ~ with the intermolecular interaction parameters by Hagler et al.31 (the HHL potential) has been employed, it has been found that the urea dimer is not stable in aqueous solutions. On the other hand, Tanaka et al.13 + Theoretical Chemistrv
Frederick, MD 0 Abstract published in Advance ACS Abstracts, July 15, 1994.
0022-3654/94/2098-8224$04.50/0
and Herndndez-Cobos et al.I9 found hydrogen bonding between urea molecules. This difference is probably due to the large difference in urea dimer potentials used in the simulations. The energy minimum of the ab initio urea dimer potentialsI3J9are about -20 and -15 kcal/mol, respectively, while the energy minimum of the GROMOS/HHL potential is -10.3 kcal/mol. In a previous work,32where we constructed NEMO potentials” for the urea-water systems, we found that the urea dimer has a minimum potential energy of -21.9 kcal/mol, indicating a strong possibilityfor finding urea complexesin aqueous solutions. Apart from our work,32 the previous urea dimer ab initio potentialsl3J7J9 have been constructed with the supermolecular approach. The disadvantagesof this approach have been discussedelsewhere.34.32 Even in the very recent ab initio potentials,17J9 the basis set superposition errors (BSSE) and dispersion energies are not included. The BSSE could be expected to be large for the rather small basis sets employed in these w o ~ k s , and ~ ~ Jthe ~ dispersion energy in the urea dimer potential minimum has been estimated to -6 kcal/mol.Q The properties of urea as a a-solute to a hydrophobic species in an aqueous solution can be explained by direct interactions between urea and the hydrophobic species.35 This model is supported in Raman difference ~pectroscopy~~ and fluorescence ~pectroscopy~~ experiments. The main argument against this model is that these properties of urea are general and should therefore be independent of the specific urea-solute interaction. This is rationalized by Nilsson et al.,38who considers the balance between all types of molecule-molecule interactions in a lattice model and concludes that this is a natural result if the urea-water interaction is slightly unfavorable compared to urea-urea and water-water interactions. The radial distribution functions obtained from a MD simulation can be directly compared with the data from neutron scattering experiments, which yield detailed information about the solution structure. As far as we know there are two such experiments on aqueous urea solutions. The first experiment by Finney and T ~ r n e r ~ presents ~ . 3 ~ the total nitrogen scattering function, GN(r),for a 2 m urea solution. This measurement has recently been extended to include 7 and 14 m urea solutions.& The second experimentz7determines the GHH(r)and GH(o,c,N)(~) radial distribution functions for a 10 m urea solution. The data analysis of the first experiment is, however, questioned by Boek 0 1994 American Chemical Society
Molecular Dynamics Simulations of Urea Solutions
TABLE 1: Repulsion Parameters for the Urea-Water System (au)' atom pair
f
00
199.4903 7.4199 199.0090 3.1390 201.2366 0.8489 11.4278 0.1949
OH NO NH
co
CH HO HH n = 16.
g 2.2014 2.0566 2.1297 1.6931 2.5467 1.3472 2.2656 1.1612
h 2.4471 2.0290 2.8229 2.0739 2.5805 1.3327 2.1850 0.1616
The same notation is used as in ref 32.
System # 21 __.._. -
-
0.25
#3 #4
and BrieW who attribute oscillations in the experimental G&) to truncation ripples in the Fourier transform. In this work we employ our recently derived potentials for the urea-water systems32in an MD simulation of a 2 m urea solution. We investigate the presence of stable urea dimers and oligomers as well as the influence of urea aggregation on a structural and dynamic properties of the solution. The results are also compared with neutron scattering data.26.27 In a molecular simulation one always has to ensure that the system is studied long enough to achieve proper time averages. Therefore we run several different simulationsoriginating from independent starting configurations. This problem is especially cumbersome for a system like a 2 m urea solutionwhere one can expect strongly bonded urea complexes with long lifetimes and relaxation times. This problem must be even more severe in simulations of larger biomolecules (e.g. proteins), where the interactions are similar to the ones treated here and the interesting part of phase space can be orders of magnitudes larger. The proper treatment of long-timecorrelations is one of the most important problems to solve in the field of molecular simulations of hydrogen bonded liquids.
2. Method A NEMO potential33is derived from intermolecular perturbation theory, where the interaction between two molecules is regarded as the perturbation. Consequently, the interaction energy can be calculated from the properties of the individual molecular wave functions. In our approach we use the HartreeFock approximation to describe each molecule. A London-type41 dispersion term is also added to the interaction energy. The detailed calculationsof the potentials employed here are described el~ewhere.3~ A preliminary set of repulsion parameters for the urea-water system were employed here compared to ref 32, and they are given in Table 1. One of the main features of a NEMO potential is that sphericalatomicdipole polarizabilitiesare adopted explicitly in the simulations,thus naturally including electrostatic many-body effects. The molecular dynamics simulations were carried out with the MOLSIM package.42 The simulated system consisted of 10 urea molecules and 277 water molecules enclosed in a cubic box of length 20.9 A, resulting in a 2 m urea solution at experimental density. Both urea and water were regarded as rigid bodies, and urea was assumed to be planar. A standard geometry was employed for the water molecule,33and an X-ray crystallography geometry was chosen for ureaS43Periodic boundary conditions were used together with a center-of-mass interaction cutoff at 8.5 A. The simulation was performed in the NVT ensemble, where the temperature was kept at 300 K by velocity rescaling. The equations of motion for the rigid molecules were integrated using quaternionsu together with a velocity Verlet algorithm.45 A first-order predictor was used to calculate the induced dipole moments.& A full iterativesolution of the induced dipole moments was generated every fifth time step in order to avoid any drift of the dipole moments. The time step used in the simulation was 0.8 fs. Two different initial systems were generated by randomly selectingmolecular positions and orientations in an extended unit
-
.
1
-'
-1.25 1 . -
20
30
40
50
60
70
80
90
100
KPS)
Figure 1. Instantenous urea-urea pair interaction energy for the four different simulations as a function of simulation time.
cell to avoid overlaps. These systems were then slowly contracted to the desired density. An equilibration period of 16 ps followed the initial configuration assembly, and subsequent production runs were carried out for 88 ps during which configurations were sampled every 0.04 ps for further analysis. These simulations will be referred to as system 1 and system 2. Since systems 1 and 2 yielded very different results in terms of atomic urea-urea radial distribution functions, their final configurations were heated to 900 K during an 8-ps period and again equilibrated for 16 ps at 300 K and another production run was carried out for 48 ps on each system. These new simulations will further on be referred to as system 3 and system 4. The total simulation time was approximately 820 CPU h on an IBM RS/6000 Model 550 computer.
3. Results and Discussion For all systems studied in this work a different average number of urea molecules involved in urea complexes was established. This indicates that slow relaxation processes present in the simulations allow urea dimers to survive for longer times than each simulation itself lasts. In order to sample all possible urea complexes (different dimers, trimers, etc.) and to achieve an ensemble average over these states would require an MD simulation that is at least a few orders of magnitude longer than in this work. This means that we do not have a complete equilibrium description of the entire phase space encompassed by an aqueous urea solution. However, during the simulated time period, there appears to be a pseudoequilibrium of that particular state for each of the systems studied. Thus, we have chosen to compare all four simulations for the properties we are interested in and examine how sensitive these properties are to urea complex formation, though we cannot say anything certain about the probability of complex formation per se. 3.1. Energetics. The potential energies of the different simulations are given in Table 2. The total potential energy of the systems, Um., are similar for the different simulations,whereas the individual components of the potential energy varies more. The largest variation is in the urea-urea pair energies whose instantenous values as a function of simulation time are given in Figure 1. Here it is seen that the urea-urea pair energy of system 1 is always lower than the corresponding energies of the other
8226 The Journal of Physical Chemistry, Vol. 98. No. 33, 1994
Astrand et al.
TABLE 3: Average Number of Urea Molecules Involved in Urea Comdexes system 1 system 2 system 3 system 4 (n) 5.8 0.4 0.2 2.3 (nh I .4 0.4 0.2 2.3 003
(n),
(n)s
3.7 0.7 0.002
0.05
simulations. This indicates that we do not have a conversion between the distributions of urea complexes for the simulated systems. The fluctuation of the urea-urea pair energy reveals time scales that are on the same order as the simulation itself. If we consider the sum of Vu. and Uuw, i.e. the total potential energy of urea in the aqueous environment, there is much less variation. This is indicative of the strong coupling of urea to its polar surroundings,regardlessoftheexact nature (waterorurea) of the local environment. It is alsonoticeable that the polarization energy is virtually constant for all the systems and amounts to almost half of the total potential energy. 3.2. Formation of Urea Complexes. The formation of urea complexes has been analyzed in several different ways. The definition of a urea dimer is arbitrary, and in this analysis we adopt a hydrogen bond distance criterion. Thus, we use one or several oxygen-hydrogen distances, ROH,to classify whether a urea dimer is present or not. In this work we have normally used ROH= 2.3 .&to be the maximum distance allowed for a hydrogen bond between two urea molecules, in order to classify the configuration as belonging to an urea-urea complex. In Table 3 we present the average number of urea molecules that are engaged in a urea complex, (n), regardless of the nature of the urea configuration. One can see that these numbers are very different for the various simulations and range from a low 296toa high60%ofallureamoleculesheingengagedinacomplex. Systems 1 and 4 have substantially more urea aggregation than the other systems. A curious fact is that system 3 is generated from system 1 and system 4 from system 2, respectively, assuring that there is little memory left of the original configuration from the annealing procedure. The average number, ( n ) , can also be divided into contributions arising fromdifferentoligomers, (n),,,, where m indicates an urea dimer, trimer, etc. If there was only one dimer present at all sampled time points, ( n ) zwould equal 2; if one dimer is not always present, ( n ) z is less than 2; and if thereis morethanonedimerpresent, (n)zislargerthan2. System 1 encompasses many different kinds of oligomers, including pentamers, but trimers predominate. In theother systemsdimers are the most prevalent complexes. This is in contradiction with the work by Hernandez-Cobos et aI.,l9 where only dimers were found. A snapshot of a urea pentamer is shown in Figure 2 in order to illustrate the nature of urea complexes. In Figure 2 it is recognized that one urea molecule hinds to three other urea molecules. It binds cyclically to two of them (a planar cyclic complex is shown in Figure 3), and it binds to the third urea molecule with a single bond from the oxygen to the cis-hydrogen of the other molecule (cis- and trans-hydrogen are defined in Figure 3). The fifth urea molecule is also bonded with a cyclic bond to one of the urea molecules. The cyclic bond between the oxygens and cis-hydrogens of two urea molecules is the energetically most stable configuration,”zand it seems that cyclic bonds, as theoneshown in Figure 3, are required to formstableoligomers in solution. One urea molecule is fully hydrated by nine surrounding water molecules.ls Given the number density of urea molecules, resultingin 1 urea molecule/28 watermolecules, it is theoretically possible to construct a system that contains fully hydrated, singly solvated urea molecules. In our case this corresponds to systems 2 and 3 which have very little urea association. Even though the mwt attractive urea-urea interaction is four times as large as that between two water molecules, thereis littledifferencebetween
Flyre 1. Snapshot in stereo of a urea pentamer connected by well developed hydrogen bonds. The shell waters surrounding the urea molecules are excluded for clarity from the figure.
rmns
FlyreB Cyclicallybonded urea dimer. This configuration corresponds to the global energy minimum of the urea dimer.
TABLE 4 RM/* 2.3
Lifetimes of Urea-Urea Hvdroeen Bonds system 1 system 2 system 3 system 4
(fnie)lps
2
0.2 2
0.2 1
4 28
a, .R is explained in the text. The dimer is prcsent during the entire simulation.
the total urea-urea and urea-water interactionsof the four systems studied. What two urea molecules lose in potential energy by dissociating is made up for by the increase in the many more possible urea-water attractive interactions. In Table 4 we present the average life times and the maximum life time of a urea-urea hydrogen bond. A urea-urea hydrogen bond is considered to be created at ROH= 2.3 and to be dissociated at a distance of . , R These data are given in Table 4 for several choices of R,.. . One can note that independently of Rmaithe life times are considerably larger for systems 1 and 4. In thesesimulations wealso find complexesthat persist during theentiresimulation. Thelife timesofthe hydrogen bondsdepend on the size of the aggregate as well as on the type of hydrogen bond developed between the urea molecules. In Table 5 we classify the different types of hydrogen bonds between urea molecules found in the simulations according to their geometry. The cyclically bonded complex (see Figure 3) dominates in systems 1 and 4 but is not present in systems 2 and 3. Other types of bonds that are present are single bonds to cisor trans-hydrogens. Notable is that the linear complex, where a urea oxygen binds to two trans-hydrogens, is not present in our simulationsof urea solutions. This is in contradiction to the work of Hernandez-Cobos et al.,’9 where linear complexes were also found. The linear complex is the only type of complex that is present in urea crystals” and is therefore used to construct the global energy minimum of empirical pair potentials. It is as far as weknowunknownwhethercrystalbondsarepresent insolution or not, but the difference between the two simulations is a very
The Journal of Physical Chemistry, Vol. 98, No. 33, 1994 8227
Molecular Dynamics Simulations of Urea Solutions
TABLE 5
Division of Urea Dimer Hydrogen Bond Types (%) system 2
system 1
ROHIA cvclicallv double bonded dbuble bbnd to cis- and trans-hydrogen double bond to two trans-hydrogens single bond to cis-hydrogen single bond to trans-hydrogen
RoHIA
ROHIA
= 2.3
= 2.8
51.1
54.2 0.2
48.8 0.1
43.8 1.8
system 3
ROHIA
= 2.3
= 2.8
RoHIA = 2.3
0.3 0.3 1.o 49.3 49.1
80.2 19.8
system 4
RoHJA = 2.8
89.1 10.9
71.3 28.7
RoHIA = 2.3
RoHIA = 2.8
78.8
74.7 0.1
21.0 0.2
19.6 5.6
20
System# 1
a)
-
# 2 _._._.
# 3 ....... # 4 ..... 15
Astrand
2t
-
f IO
5
0
0 '
0
1
2
3
4 rlA
5
7
6
81
0
8
#
2
#3 #4
n
2
3
4 dA
6
5
I
b
? .
I
System# 1
1
I1
_.._.
System# 1
b)
-
# 2 -.-_.
#3 #4 Astrand
---
2
v
M
1
01 0
~
1
2
3
4 rIA
5
6
1
0
8
Figure 4. Atomic radial distribution functions for urea-urea: (a) O-Hh;
2
3
4 rlA
~~
5
6
1
3 r
(b) 0-H,,,. The remaining urea-urea atomic radial distribution functions are available from the authors on request.
strong argument for using ab initio derived potentials, where no bias of the solid or liquid phase is present. Using an effective potential derived from crystal data, such as the GROMOS/HHL potential,30.31 will not allow for even a qualitative description of solution properties. 3.3. Atomic Radial Distribution Functions. A selection of intermolecular urea-urea radial distribution functions (RDF) is presented in Figure 4 in order to illustrate the difference between the four different simulations. In Figure 4a we show the O-H,i, RDF. The peak at 1.85 8, stems from the cis-hydrogen involved in the hydrogen bond and the peak at 5.75 8, from the cis-hydrogen on the other side of the same urea molecule. It is also clearly seen that ordered urea complexes are more prevalent in systems 1 and 4. The 0-H,,, RDF is shown in Figure 4b. The peaks at 3.55 and 5.35 A for systems 1 and 4 stem from the trans-hydrogens in a complex similar to that of Figure 3. For systems 2 and 3 the RDFs are broad and featureless with a maximum at 3.95 8,. For these systems the nearest urea molecules are mainly situated outside the first hydration shell. The radial distribution functions for urea-cis-hydrogen, ureatrans-hydrogen, and urea-oxygen relative to the water-oxygen are presented in Figures 5a-c. The H&-O, RDF and the second
1
2
I
System# 1 #2
Astrand Bock
-
1
.-..--- ~
c,
I 0
8
I
2
, 3
4 rlA
5
6
7
8
Figure 5. Atomic radial distribution functions for urea-water: (a) &Ow;(b) H-4,; (c) Ou-Ow. The remaining urea-water atomic radial distribution functions are available from the authors on request. The results from a previous simulation of one urea molecule in water is labeled Astrand.l5 In our previous work's the radial distribution functions for H h and H- in Figure 7 have unfortunately been exchanged with each other. In Figure 5c the results from Bock and BrieW are also given.
and third peaks in the H,,,-O, RDF differ between the four different simulations, but the first peak in the H,,s-O, RDF is
8228
The Journal of Physical Chemistry, Vol. 98, No. 33, 1994
Astrand et al. System# 1 # 2 ._.._. # 3 ....... #4 .... Finney Astrand - - -.
2.5
3
Bo& .......
2 4
2
1
-0
0
1
2
3
4 r/A
5
6
8
1
Figure 6. Atomic radial distribution function for water oxygen-water oxygen. The remaining water-water atomicradial distributionfunctions are availablefrom the authors on request. The results from our previous urea simulation on a dilute solution are labeled AstrandI5and from a simulation of pure water is labeled Wallqvist$* respectively. 4
-
all shell -----. bulk .......
3
G 2
I
0 0
1
2
3
4
5
6
7
8
dA
Figure 7. Atomic radial distribution function for water oxygen-water oxygen. All denotes that all water molecules have been included, bulk denotes that only water molecules belonging to the bulk region are included, and shell denotesthat only water molecules belonging to the shell region are included. The first minimum in water oxygen-urea atom radial distribution functions are used to divide the water molecules into bulk and shell water molecules. The results are shown for system 1; identical results are obtained for systems 2 4
identical in all four simulations. From this we can conclude that urea molecules preferably bind to the &hydrogens of other urea molecules, not to trans-hydrogens which are engaged in hydrogen bonding to the surrounding waters. Comparison with our previous simulation15 gives that the first peak of the 0,-0, RDF in Figure 5c has been shifted from 2.65 A outward to 2.75 A in this work. On the other hand, the structure around the urea hydrogens is more pronounced in this work than p r e v i o ~ s l y . This ~ ~ is a consequence of the difference in the potentials employed in the different simulations. Comparison with the simulation by Boek and Briels16 indicates that the water structure around the urea molecules is much less pronounced in their work, which is due to the different urea-water potentials employed in their simulation. In Figure 6, which shows the O,-O, RDF, one can see that the water structure is almost identical in all four simulations. The first peaksituatedat 2.70Aand thesecondat4.35Aareidentical toour previous simulationsof a dilute urea solution15and a NEMO simulation of pure water.48 In Figure 7 we show the RDFs for bulk water and for water molecules in the first hydration shell around the urea molecules. These graphs are indistinguishable from each other for the different systems studied, indicating the similar nature of urea hydration regardless of the urea complex formed. As in a previous workI5 the shell region is defined as the
0
1
2
3
4 r/A
5
6
1
8
Figure 8. Total nitrogen radial distribution function, h ( r ) The results from a neutron scattering experiment by Finney and Turner26is labeled Finney, and simulation results from our previous simulation on a dilute solution15 and the work by Bcek and BrieW are labeled Astrand and Boek, respectively.
water molecules whoseoxygen atom lies within the first minimum of the urea-atom water-oxygen RDFs. It is clearly noticed that the water molecules in the shell region maintain their liquid water structure, which promotes the local nature of urea solvation. 3.4. Comparison with Neutron Scattering Experiments. In principle it is possible to determine detailed structuralinformation of complex liquids from neutron diffraction studies. With utilization of isotope substitution, partial radial distribution functions may be evaluated from the difference between two such experiment^.^^ Such a partial radial distribution function is written as
where
(3) q is the concentration of atom i, and bi its scattering length. M' denotes a different isotopeof element M . Dividing by the constant term 2 c ~ ( -b b~~ fgives )
which is then normalized to 1 a t large r. Finney and Turner26 have measured the total nitrogen radial distribution function, G ( r ) , in a 2 m urea solution by substituting I4Nto 15Nand using deuterated urea and water samples. Thisquantity can be directly evaluated from atomic RDFs calculated in MD simulations and experimental scattering lengths.50 It should be noted that while we have simulated a 2 m urea/HzO solution, the experimentz6 is carried out on a 2 m urea[D4]/DzO solution, which yields a slightly different concentration. The calculated GN(r), presented in Figure 8, appear identical for systems 1-4 and could therefore not be used to differentiate the extent of urea complexation. This is in agreement with the experiments by Turner et a1.N where they find similar GN(r) over a range of different urea concentrations (2.0-14.0 m). The maximum standard deviation in the calculated G ( r ) is as small as 0.02. In the experiment conducted on the 2.0 m urea solution26 there are several peaks arising from nitrogen water correlations. Intramolecular correlations a t 2.25 and 3.25 A obscure the intermolecular contributions. Truncation ripples from the Fourier transform of the neutron scatteringdata
The Journal of Physical Chemistry, Vol. 98, No. 33, I994 8229
Molecular Dynamics Simulations of Urea Solutions 2
2
-
Total G(r) contributions from urea -----. contributions from water .......
t
15
-
81
01 0
I
1
2
A'
I
3
4 r/A
5
6
I
I
0
8
2
1
3
4
5
6
8
7
r/A
Figure 9. Total nitrogen radial distribution functions, G&), divided into contributions from urea and water molecules. The results are shown for system 1; identical results are obtained for systems 2-4.
TABLE 6 Atom-Pair Contributions to the Total Nitrogen Radial Distribution Function, C d r ) peak positionlbr
atom pair
2.85 3.35-3.65 4.05 4.45 4.85-4.95
N-OW, N-O,, N-Hs N-Hw, N-N, N-C, N-HtraM N-OW N-HtraM N-O,, N-O,, N-N
have been shown by Boek and Briels16 to give spurious results similar to the experimentalcorrelationsat 2.67 and 3.75 A, leaving the correctly interpreted intermolecular GN(r) as a rather broad and featureless peak centered roughly at 3.5 A. In our work the intramolecular contributions are obviously absent since we use rigid molecules. In general our distributionsare rather featureless, which is in agreement with the recent simulation data by Hernlndez-Cobos et al.19 We do find a small shoulder at 2.85 A arising mainly from N-O correlations in the RDF of ureaurea and urea-water. The main peak is located at 3.55 A compared to 3.65 8, in our previous work's and compared to 3.45 A in the work by Boek and Briels.I6 Intramolecular flexibility, allowed for in the simulation by Boek and Briels,16 shifts the broad maximum at 3.5 inward to shorter distances. We also find a peak at 4.95 A, as found elsewhere,lgbut a comparison with experimental data is difficult since the analysis of the experiments is truncated at 5 A. It is not meaningful to make the retransformation of our GN(r) according to Boekand Briels,16 since the truncation ripples will not appear if there are no intramolecular correlations present. The simulated GN(r) of system 1 is divided up into contributions from water and urea in Figure 9. The main contributions stem from urea-water interactions, but an interesting feature is that the peaks of the urea dimer contribution coincide with the peaks from the ureawater contribution (see Table 6 ) . This is an additional proof that urea is accommodatedin liquid water without extensivedisruption of the local or long range ordering of water molecules. We have also compared with another experimentz7where the GHH(r) and GH(c,o,N)partial RDFs were measured in a 10 m urea solution by using different fractions of DzO/HzO and urea[D4]/ urea. These partial RDFs are calculated as
which are normalized to 1 at large r. As for GN(r) we find very similar GHH(r), in Figure 10,and GH(c,o,N),in Figure 12,partial
Figure 10. Partial radial distribution function, &H(r), given for the four simulated systems 1-4, a 10 m urea solution from the experimental data of Finney et aL2' The hydrogen-hydrogen radial distribution functions is also given for two pure water systems, a simulated system32 employing the model potential used here and the experimental data set of Soper and Phillips.51
-
Total G(r) urea-urea contributions --..-. urea-water contributions ....... water-water contnbutions 1.5
-
8 1
0.5
0
2
1
3
4 rIA
5
6
7
8
Figure 11. Partial radial distribution functions, GHH(r), divided into contributions from urea and water molecules. The results are shown for system 1; identical results are obtained for systems 2-4. 2
-
System# 1 # 2 ...... ....... 1.5
81 h
0.5
01 0
1
w
I
I
2
3
4 rIA
5
6
7
8
,
Figure 12. Partial radial distribution function, GH(c,o,N)(~),given for the four simulated systems 1-4; a 10 m urea solution from the experimental data of Finney et al.27 Theoxygen-hydrogen radial distribution functions is also given for two pure water systems, a simulated system32 employing the model potential used here and the experimental data set of Soper and Phillip~.~'
RDFs for the four different simulations. However, here it is mainly water-water correlations that contribute to the RDFs, as is evident from the different contributions to the RDF given in Figures 1 1 and 13, and consequently they contain very little
8230 The Journal of Physical Chemistry, Vol. 98, No. 33, 1994
Astrand et al. 7 .
-
1.5
Total G(r) urea-urea contributions ..---. urea-water contributions ...... water-water contributions
-
0
'
1
2
3
4 r/A
5
6
7
#4
'
8
Figure 13. Partial radial distribution functions, GHH(c,o,N)(~), divided into contributionsfrom urea and water molecules. The results are shown for system 1; identical results are obtained for systems 2-4.
0
1
2
3
4 r/A
5
6
7
8
Figure 14. Total carbon radial distribution function, Gc(r). 7 .
TotalG(r) contributions from urea -----.
information about urea complexation. Even though the expericontributions from water ....... was carried out in a 10 m urea solution, we note in Figure 10 that there is a quantitative agreement for GHH(r),except for the peak at 2.5 A. We find the first peak a t 2.4 A equal to the experimental pure water peak5l (as well as the first peak of the s i m ~ l a t e dpure ~ ~ , ~water ~ system using the model potential employed here), compared to 2.5 A in the experimental G H H ( ~ )The . ~ ~height of the first peak is noticeably larger here than in the diffraction experiments of pure wateP1 and in urea s0lution,2~but similar to the M D simulation of NEMO We find the second peak at 3.75 A in close resemblance to the e ~ p e r i m e n t a one l ~ ~at 3.80 A. The GH(c,o,N)RDF given in Figure 12 is also in good agreement 0 1 2 3 4 5 6 7 8 with the available experimental data.Z7 As previously discussed, r/A the dominant contributions stem from water-water correlations, Figure 15. Total carbon radial distribution functions,G&), divided into as deduced from Figure 13. The main impression is that the contributionsfrom urea and water molecules. The results are shown for calculated GH(c,o,N)is very similar to the experimental go&) for system 1; identical results are obtained for systems 2-4. pure water31 as well as for NEMO water.48 The first peak is situated at 1.80 A compared with 1.75 8, in the e ~ p e r i m e n t . ~ ~ at 1.90 A, and urea-oxygen peaks in the range 1.80-1.85 A. If the first GHH(r) peak is resolved in the same way, it is noticed The experimental peak is however much broader and asymmetric that the main contributions are a H,-H, peak a t 2.4 A and the which can be rationalized by considering the intramolecular H-(C,O,N) correlations in urea which are present in the H,i,-H, and Htrans-Hwpeaks a t 2.5 A compared with the first experiment but not in our simulations. For pure water the first gHH(r) peak in pure water a t 2.4 A. This would mean that both peak in g&r) is experimentally found at 1.85 A.51 However, the 0,-H, hydrogen bond distances are slightly longer than in for the two MD simulations on NEMO water this peak was found 0,-H, and that the H,-H, distances are slightly longer than the at 1.80 32 and 1.85 A,'@respectively. Since it is most consistent 0,-H, distance, which is in contradiction with the suggested to compare with the water simulation where we have adopted interpretati0n.2~9~~ An equally goodinterpretation of the increased repulsion between the water-hydrogen and the -NH2 group in exactly the same water potential as here,32we are unable to support urea, compared with water dimer interactions, would simply be the conclusion by Finney et al.27that the first GH(c,o,N)peak is shifted slightly inward in a urea solution compared with pure that the -NH2 group has a different orientational dependence of the electrostatic properties compared to a water molecule. water. The second peak is relatively similar in all cases, and its position ranges from 3.25 to 3.35 A. It has been arg~ed~~952 However, we would like to stress that the uncertainties in the that calculated RDFs are of the order 0 . 1 A and that the arguments the structure of a hydrogen bonded dimer could be determined both elsewherez7 and here are based on very subtle differences by a balance between the electrostatic attraction between a in the RDFs. One example of this is the difference between the hydrogen bonded proton and a water-oxygen and a repulsion two MD simulations on NEMO ~ a t e r , ~where * . ~ the ~ first go&) between a water-hydrogen and the hydrogen bonded proton, as peak is shifted from 1.85 to 1.80 A. Even though the inclusion well as the atom bonded to the hydrogen bonded proton (a nitrogen atom in the case of urea). If the hydrogen bond distance is reduced, of the damping term on the dispersion energy and the addition of two extra terms to the expansion of the repulsion energy were the increased repulsion could be balanced by a rotation of the necessary in order to obtain a reliable urea-urea p~tential,'~ it proton-accepting water molecule from a tetrahedral to a more only amounts to a refitting of the dispersion and repulsion this planar configuration. In the experiment by Finney et is noticed in urea solutions since the position of the first GH c . 0 , ~ ) parameters to the same potential surface for the water dimer. This emphasizes the importance of developing accurate potentials peak has been shifted from 1.85 A in pure water to 1.75 and for molecular simulations in order to make a quantitative the position of the first GHHis shifted from 2.4 A in pure water interpretation of experimental data. to 2.5 A. The conclusions from our simulation results are however quite different. Even though we find the first GH(c,o,N)peak at We have also calculated the total carbon radial distribution function, G&), by employing eqs 1-3. In Figure 1 4 one can see the same position as the first goH peak for pure water (1.80 A), that it is almost identical in all four simulations, and in Figure the contributions to the first GH(c,o,N)peak stem from an 0,-H, 15 it is noticed that the major contributions to Gc(r) stem from peak a t 1.80 A, an Ow-Htrans peak a t 1.85 A, an Ow-Hcispeak
8,
The Journal of Physical Chemistry, Vol. 98, No. 33, 1994 8231
Molecular Dynamics Simulations of Urea Solutions
TABLE 7: Diffusion Coefficientsa ( lo9 m*/s) urea water urea complexes free urea bulk water shell water
system 1
system 2
system 3
system 4
Astrand et aleb
Tanaka et a1.f
expt
0.2 0.8 0.2 0.3 0.8 0.6
0.4 0.9 0.4 0.3 0.9 0.8
0.4 0.8 0.3 0.4 0.8 0.7
0.2 0.8 0.2 0.3 0.8 0.7
0.5
1.o 2.0
1.4c 2.6d
0.5 1.o 0.7
(1/6)[A($)(t)/At] in the interval 1 S t I7 ps. See ref 15. See ref 12. See ref 53. Center-of-mass distances are used to divide the molecules into different groups: ruu = 5.5 A; rw = 4.5 A.
TABLE 8: Orientational Decay Times (ps)’
System# 1 # 2 _..... # 3 .......
system 1 system 2 system 3 system 4
I
0‘ 0
1
2
3
4 rIA
5
IO
6
7
8
shell water
free urea
urea complexes
in-plane dipole out-of-plane in-plane dipole out-of-plane in-plane dipole out-of-plane in-plane dipole out-of-plane in-plane dipole out-of-plane
13 16 10 13 14 10 14 17 11 14 16 11 14 15 11
15 21 [20] 12 13 15 9 (10) 14 16 11 14(15) 17 11 23 31 22
19 39 35 57 63 41 52 89 49 46 56 36 29 81 24
102 139 75 47 70 43 40 73 45 69 104 64
a Several different criteria for dividing the water molecules into shell and bulk water molecules have been employed: (1) minima in g(r) of urea-atom and water-oxygen (as in previous work’s); (2) minima in all atomic urea-water g(r) (given in parentheses ( ) if they are different from 1); (3) center-of-mass distance between urea and water (given in brackets [ ] if they are different from 1). A center-of-mass distance (ruu = 5.5 A) is used to divide the urea molecules into free urea and urea complexes. See ref 15.
#4
-
bulk water
~
System# 1 # 2 _..... #3
b)
p
previous workb
axis
5 -
0
1
2
3
4 rlA
5
6
7
8
Figure 16. Atomic radial distribution functions: (a) g ” r ) and (b) gccw.
urea-water interactions. The first peak is situated at 2.85 A and the second at 4.05 A. Gc(r) for urea solutions has not been measured by neutron diffraction experiments, but an experiment on methanol has been reported.40 The coordination numbers are, however, too low by a factor of more than 2, but carbon isotope substitution may in the near future provide more experimental information on the structure of aqueous urea solutions. An interesting experiment to determine the degree of urea complexation would be to investigate GN(r) or Gc(r) at different concentrations and in line with the GH(r) experimentz7 separate G N ( or ~ G c W into mN(r) and GN(H,C,o)(r) or gccW and GqH,N,O)(r), respectively. The g” and gcc RDFs are, as seen in Figure 16, highly dependent on the number of urea oligomers in the solution. 3.5. Dynamics. Diffusion coefficients and the decay times extracted from the orientational autocorrelation functions are presented in Tables 7 and 8. The diffusion coefficients are, as in earlier simulations with NEMO potentials,15.48 much lower than the experimental diffusion coefficients.12.53 They are also significantly smaller than in a previous MD simulation of one urea in water,lSdue to the changes in the intermolecular potentials of the urea-water system.32 The diffusion coefficients are also
different for the different simulations, indicating that they are dependent on the degree of complexation in the system. One conclusion that may be drawn is that the diffusion coefficients are somewhat smaller for shell water, than for bulk water, implying that the hydration structure around urea is more rigid than the structure in pure water. However, the difference is only about lo%, which further confirms that urea is well accommodated in water. The orientational decay times given in Table 8 are in general much larger for urea than for water. Some of the relaxation times for urea are even larger than the simulation times, which emphasizes the fact that an adequate sampling of the phase space of an aqueous urea solution would require a simulation time of several nanoseconds. Shell water molecules have different decay times compared with our earlier simulation,ls where the decay time for shell water was about twice as big as for bulk water. This points to the fact that the dynamics of a system is very sensitive to the potential used in the simulation. 3.6. Induced Dipole Moments. The molecular induced dipole moments (IDM) are given in Table 9. The IDM for urea is on the average 1.95 D, which can be compared with the static dipole moment of 5.31 D. The largest contribution, 70%,to the IDM is from the highly polarizable (1.43 A3) oxygen atom, with the nitrogen and hydrogens contributing equally to the remaining 30%. Regardless of the environment of the urea molecule, the individual atom contributions remain the same, and thus the total IDM of urea is the same in all our simulations. Consequently the electrostatic environment developed around a urea molecule in solution is similar for all our studied systems. Comparison of the electrical field magnitude &/a) between urea-oxygen and water-oxygen exhibit a high degree of similarity, confirming that the urea oxygen atom is in an electrostatic environment similar to that of a water oxygen atom. That may mean that it will be
8232
Astrand et al.
The Journal of Physical Chemistry, Vol. 98, No. 33, 1994
TABLE 9: Molecular Induced Dipole Moments (D)* all urea cyclically double bonded single bond to cis-hydrogen single bond to trans-hydrogen free urea all water shell water bulk water
system 1
system 2
system 3
system 4
1.94(0.29) 1.84(0.26) 1.93(0.26)
1.96(0.28)
1.97(0.28)
1.96(0.21) 1.89(0.24) 1.96(0.28) 0.93(0.19) 0.90(0.19) 0.94(0.19)
1.95(0.20) 1.87(0.18) 1.97(0.28) 0.94(0.20) 0.92(0.20) 0.95(0.19)
1.91(0.28) 1.92(0.24) 1.97(0.22)
2.01 (0.28) 0.94(0.19) 0.91(0.19) 0.94(0.19)
1.90(0.27) 0.94(0.19) 0.91(0.20) 0.95(0.19)
pure water*
0.84
,IThe standard deviation is given in parentheses. See ref 48. ROH= 2.3 A is used to divide the urea molecules into different groups. The minimum in water oxygen-urea atom RDFs is used to divide the water molecules into different groups. hard to decide whether urea dimers are present or not in a Raman spectroscopy experiment, where the polarizability derivative with respect to the nuclear coordinates is measured. If the dominant effect on d a / d r originates from the electric field, then the recent experiment by Hoccart and T ~ r r e 1 may 1 ~ ~not be inconsistent with the results presented here. The IDM of water is 0.94 D compared with the static dipole moment of 2.04 D. It is considerably larger than the IDM for pure water in a previous NEMO ~ i m u l a t i o nwhich , ~ ~ is due to a reparametrization of the p~larizabilities.~~ Water molecules belonging to the bulk region have a slightly larger induced dipole moment than water molecules belonging to the shell region. The standard deviations of the molecular IDMs, which describes the fluctuations of the dipole moment, are also given in Table 9. They are on the average 0.2 D for water and 0.3 D for urea. The dipole moment fluctuations for complexed urea is consistently smaller than for free urea.
4. Conclusion Molecular dynamics simulationsare a powerful tool for studying structure, thermodynamics, and dynamics of condensed matter but require that accurate intermolecular force fields are employed and that the phase space is adequately sampled. In this work we use a more accurate set of potentials than in previous simulations of aqueous urea solutions. It is important to note that some of the previous investigationsl4,l619 used a urea-urea potential that is not in qualitative agreement with ab initio calculations of the interaction strength. Extensive calculations show that the urea dimer can be stabilized with up to 22 kcal/mol upon formation of the cyclic dimer relative to the isolated species, compared to 10-13 kcal/mol in the GROMOS/HHLS0J1and OPLS potentials.” Even though the OPLS parameters are generated from quantum chemical calculations, they only give a minimum potential energy of -12.52 kcal/mol.17 The time scales associated with the urea complexation process cannot be adequately sampled in a single molecular dynamics calculation due to a lack of phase-space sampling. Decay times of complexation on the order of 100 ps found here imply sampling times on the order of nanoseconds. The four different systems studied do not represent the complete phase space but only an undetermined part of it. The lackof conversion between systems at 300 K is also symptomatic of long time scales in the problem. Dimer formation in a concentrated solution is very likely d u e to the spatial proximity of molecules and the large interaction energy possible upon dimerization. There are a number of possible urea dimer complexes that have energies considerably more attractive than a single water hydrogen bond, thus compensating the unfavorable entropy term associated with arranging this configuration. The urea-water interactions used here have been shown to result in a negligible perturbed water shell, with very little water structuring associated with urea solvation, regardless of the urea concentration. This is consistent with the SKSS where formation of urea dimers and polymers were used to explain thermodynamic data for aqueous urea solutions. As pointed out by Finer et the SKSS model would require urea-urea bonds with considerably longer life times than those
normally ascribed to hydrogen bonds, which is certainly found here. Furthermore, the conclusion that urea molecules do not disturb the local water ordering discredits the FF model.9 Comparison between results from MD simulations and neutron scattering data could yield a stringent test for the potential employed in the simulations. On the other hand, detailed information from a simulation can also aid in the understanding of the experimental data and its interpretations. Computer simulations can be used as a time-saving tool to probe which neutron scattering experiments can provide the most useful structural information about an aqueous solution. In this work it is clearly noticed that a measurement ofglyN(r)or g&) would give an undisputed experimental answer to which degree urea forms complexes in aqueous solutions.
Acknowledgment. P.-O.A. thanks Dr. P. Lime for assistance with running the MOLSIM package. References and Notes (1) Ellerton, H. D.; Dunlop, P. J . Phys. Chem. 1966, 70, 1831-1837. (2) Wetlaufer, D. B.; Malik, S.K.;Stoller, L.; Coffin, R. I.J.Am. Chem. SOC.1964,86, 508-514. (3) Brunig, W.; Holtzer, A. J . Am. Chem. SOC.1961, 83, 4865-4866. (4) Brandts, J. F.; Hunt, L. J. J . Am. Chem. Soc. 1967,89,4826-4838. (5) Makhatadze, G. I.; Privalov, P. L. J. Mol. Biol. 1992,226,491-505. (6) Schellman, J. A. Comput. Rend. Trav. Lab. Carlsberg Ser. Chim. 1955, 29, 223-229. (7) Kresheck, G. C.; Scheraga, H. A. J. Phys. Chem. 1965,69, 17041706. ( 8 ) Stokes, R. H. Aust. J . Chem. 1967, 20, 2087-2100. (9) Frank, H. S.;Franks, F. J. Chem. Phys. 1968, 48,4746-4757. (10) Kuharski, R. A.; Rossky, P. J. J . Am. Chem. SOC.1984,106,57865793. (1 1) Kuharski, R. A.; Rossky, P. J. J . Am. Chem. SOC.1984,106,5794 5800. (12) Tanaka, H.; Touhara, H.; Nakanishi, K.; Watanabe, N. J . Chem. Phys. 1984, 80, 5170-5186. (13) Tanaka, H.; Nakanishi, K.; Touhara, H. J . Chem. Phys. 1985,82, 5184-5191. (14) Cristinziano, P.; Leu, F.; Amodeo, P.; Barone, V. Chem. Phys. Len. 1987. 240.401405. (IS) Ahrand, P.-0.; Wallqvist, A,; Karlstram, G.; Linse, P. J . Chem. Phys. 1991, 95, 6395-6396. (16) Boek, E. S.;Briels, W. J. J. Chem. Phys. 1993, 98, 1422-1427. (17) Duffy, E. M.; Severance, D. L.; Jorgensen, W. L. Isr. J . Chem. 1993, 33. 323-330. (18) Duffy, E. M.; Kowalczyk, P. J.; Jorgensen, W. L. J . Am. Chem. SOC. 1993, 115, 9271-9275. (19) Hernhdez-Cobos, J.; Ortega-Blake, I.; Bonilla-Marin, M.; MorenoBello, M. J. Chem. Phys. 1993, 99, 9122-9134. (20) Finer, E. G.; Franks, F.; Tait, M. J. J. Am. Chem. SOC.1972. 94, 4424-4429. (21) Adams, R.; Balyuzi, H. H. M.; Burge, R. E. J. Appl. Crystallogr. 1977, 10, 256-261. (22) Bloemendal, M.; Somsen, G. J . Am. Chem. SOC.1985, 107, 34263431. (23) Piekarski, H.; Somsen, G. Can. J. Chem. 1986,64, 1721-1724. (24) Hoccart, X.;Turrel, G. J. Chem. Phys. 1993, 99, 8498-8503. (25) Subramanian, S.; Sarma, T. S.;Balasubramanian, D.; Ahluwalia, J. C. J. Phys. Chem. 1971, 75, 815-820. (26) Finney, J. L.; Turner, J. Electrochim. Acta 1988, 33, 1183-1 190. (27) Finney, J. L.; Soper, A. K.; Turner, J. Physica B 1989, 156-257, 151-153. (28) Jakli, G.; van Hook, W. A. J. Phys. Chem. 1981, 85, 3480-3493. (29) van Gunsteren, W. F.; Berendsen, H. J. C. Groningen molecular simulation library, 1987.
Molecular Dynamics Simulations of Urea Solutions (30) Hermans, J.; Berendsen, H. J. C.; van Gunsteren, W. F.; Postma, J.
P. M. Biopolymers 1984, 23, 1513-15I8.
(31) Hagler, A. T.; Huler, E.; Lifson, S . J . Am. Chem. SOC.1976, 96, 5319-5327. (32) Astrand, P.-0.; Wallqvist, A.; Karlstrbm, G. J. Chem. Phys. 1994, 100, 1262-1273. (33) Wallqvist, A,; Karlstrbm, G. Chem. Scr. 1989, 29A, 131-137. (34) Astrand, P.-0.; Wallqvist, A.; Karlstrbm, G. J. Phys. Chem. 1991, 95, 8419-8429, (35) Roseman, M.; Jencks, W. P. J. Am. Chem. SOC.197597,631440. (36) Mizutani, Y.;Kamogawa, K.; Nakanishi, K. J. Phys. Chem. 1989, 93, 5650-5654. (37) Sarkar, N.; Das, K.; Nath, D.; Bhattacharyya, K. Chem. Phys. Lett. 1992, 196,491-496. (38) Nilsson, S.;Piculell, L.; Malmsten, M. J. Phys. Chem. 1990, 94, 5149-51 54. (39) Finney, J. L.; Turner, J. Ann. N . Y.Acad. Sci. 1986,482, 127-142. _ _ (40) - _ Turner, J.; Finney, J. L.; Soper,A. K. Z . Naturforsch 1991, 46a,
73-83.
(41) Gray, C. G.; Gubbins, K. E. Theory of molecularfluids; Clarendon: Oxford. U.K.. 1984: Vol. 1. (42) Linse, P.; Wallqvist, A. MOLSIM 1.2, Lund University: Sweden, 1991.
The Journal of Physical Chemistry, Vol. 98, No. 33, 1994 8233 (43) Andrew, E. R.; Hyndman, D. Discuss. Faraday Soc. 1955,19,195200. (44) Allen, M. P.; Tildesley, D. S. Computer simulations of liquids; Clarendon: Oxford, U.K., 1987. (45) Swope, W.C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. J . Chem. Phys. 1982, 76, 637-649. (46) Ahlstrbm, P.; Wallqvist, A.; Engstrbm, S.;Jbnsson, B. Mol. Phys. 1989,68, 563-581. (47) Dovesi, R.; Causa, M.; Orlando, R.; Roetti, C.; Saunders, V. R. J. Chem. Phys. 1990, 92, 7402-7411. (48) Wallqvist, A.; Ahlstrem, P.; Karlstrbm, G. J. Phys. Chem. 1990,94, 1649-1656. Erratum in: J . Phys. Chem. 1991, 95, 4922. (49) Soper,A. K.; Neilson, G. W.; Enderby, J. E.; Howe, R. A. J . Phys. 1977, ClO, 1793-1801. (50) Dore, J. C. In NATO Advanced Study Institute on ‘Molecular Liquids”; NATO: Florence, Italy, 1983. (51) Soper, A. K.; Phillips, M. G. Chem. Phys. 1986, 107,4740. (52) Savage, H. F. J.; Finney, J. L. Nature 1986, 322, 717-720. (53) Krynicki, K.; Green, C. D.; Sawyer, D. W. Discuss Faraday SOC. 1978,66, 199-208.
