Theoretical Studies of Concentrated Aqueous Urea Solutions Using Computacional Monte Carlo Simulation

Article

Theoretical Studies of Concentrated Aqueous Urea Solutions Using Computacional Monte Carlo Simulation

Celso A. Bertran^*a, José J. V. Cirino^b and Luiz C. G. Freitas^c

^a

Instituto de Química, Universidade Estadual de Campinas, CP 6154, 13083-970, Campinas - SP, Brazil

^b

Departamento de Física ¾ ICE, Universidade Federal de Juiz de Fora, 36036-330, Juiz de Fora - MG, Brazil

^c

Departamento de Química, Universidade Federal de São Carlos, CP 676, 13565-905, São Carlos - SP, Brazil

Simulações de Monte Carlo foram realizadas para soluções aquosas concentradas de uréia nas concentrações de 4, 5, 6, 7 e 8 mol L

^-1

no

ensemble

NpT a 298 K e 1 atm. As superfícies de energia potencial foram representadas pela soma dos potenciais clássicos de Coulomb e de Lennard-Jones, nos quais a aditividade de pares foi considerada. Nas simulações a água foi descrita pelo modelo TIP4P e a uréia por dois modelos: o modelo OPLS planar e um modelo no qual a molécula de uréia não é totalmente planar. Os resultados das simulações para os dois modelos mostraram que a uréia não causa perturbações significativas na estrutura da água, na faixa de concentrações estudadas. No entanto, foi observada diminuição do número médio de ligações de hidrogênio entre as moléculas de água. Para as soluções mais concentradas os resultados obtidos com o modelo OPLS sugerem a formação de dímeros de uréia. Observou-se similaridade entre alguns dos resultados das simulações obtidos com os modelos: OPLS planar e não-planar da uréia , nas mesmas concentrações e condições.

Monte Carlo simulations were carried out for concentrated aqueous urea solutions at 4, 5, 6, 7 and 8 mol L

^-1

in the NpT ensemble at 298 K and 1 atm. The potential energy surfaces were represented by the sum of Coulomb and Lennard-Jones potentials in which pair-wise additivity was considered. Water was described by the TIP4P model, while urea was described by the OPLS planar model and also by a proposed non-planar model. The simulation results for both urea models showed that urea does not induce any significant changes in the structure of water over all the concentrations studied. However, it was noted that the number of hydrogen bonds between the water molecules decrease as the concentration increases. The simulation results of the OPLS planar model suggest the formation of urea dimers in the more concentrated solutions. The simulations with the non-planar model showed several features similar to those of the OPLS planar model.

Keywords:

urea, Monte Carlo, liquid simulation

Introduction

Urea is a substance that has been studied for many years due to its important physico-chemical properties, especially in concentrated aqueous solutions. Properties such as protein denaturation,¹ solubility enhancement of hydrocarbons in water,² inhibition of micellar aggregation³ and formation of quasi-ideal solutions with water⁴ have been associated with urea in aqueous solution. Recently, urea was reported to play an important role in the homogeneity of some ceramic materials prepared by sol-gel routes.^5-7

Various experimental studies of aqueous urea solutions were published in the last 30 years, trying to establish a model for urea's action in terms of urea-water interactions, but until now, the published results are not in agreement. Frank and Franks,⁸ for instance, pointed out that the hydrogen bonds between water molecules are destroyed. Infrared studies⁹ showed an increase in the intensity of the O-H stretching band, which also suggests a weakening of the hydrogen bonds between water molecules as a result of urea addition. On the other hand, Stokes⁴ affirmed that the water structure is essentially unchanged and that urea molecules can form dimers in solution. Raman studies¹⁰ did not show any evidence of urea dimers in dilute solutions. More recently, Hoccart and Turrell¹¹ published an infrared study of aqueous urea solutions. The authors concluded that urea does not form dimers in solution, even in more concentrated solutions. They also suggested that the water structure is a result of dipole-dipole interactions instead of hydrogen bonds.

In the field of computational chemistry the theoretical studies of Kuharski and Rossky,¹² Nakanishi,¹³ Tanaka et al.,¹⁴ Cobos et al.¹⁵ and Åstrand et al.¹⁶ showed no influence of urea on the structure of water. These authors suggested that urea was able to enter the structure of water causing little perturbation.

In the present work, Monte Carlo statistical simulations were used to investigate the thermodynamics and structure of aqueous urea solutions over a wide range of concentrations. Two urea molecular geometries were used: a planar one derived from crystallographic determinations on solid urea, and the another one, slightly non planar, as predicted by ab initio calculations.

Computational Details

Intermolecular potential function

The potential energy between two molecules a and b is described by the sum of the Coulomb and Lennard-Jones potentials,

where r_ij is the distance between site i on molecule a and site j on molecule b, and q_k is the partial point charges of the sites. The parameters s_ij and e_ij for non-diagonal interactions were obtained by the same combining rules used by Duffy et al.: ¹⁷

Models for urea

Two different models for urea molecules were used in this work. The geometry of the first model was obtained from published crystallographic data¹⁸ and for the other model, the ab initio geometry was calculated at the MP2/MC-311G(d,p) level.¹⁹ While the crystallographic data show an essentially planar structure for urea, the ab initio calculation predicts a non-planar structure.¹⁹

The Lennard-Jones parameters and partial point charges for the planar urea OPLS model were obtained from the work of Duffy et al.¹⁷ For the non-planar model, the same Lennard-Jones parameters were used and the point charges were obtained from ab initio calculation with the CHELPG²⁰ procedure. Both models are summarized in Table 1. For water molecules the TIP4P model²¹ was used.

Monte Carlo Simulations

Thermodynamics and structural properties of concentrated aqueous urea solutions were determined with the use of the Diadorim²² program, which employs a Metropolis Monte Carlo method. Periodic boundary conditions in the NpT ensemble were adopted at 298 K and 1 atm and the minimal image convention was used. Cubic boxes with 360 molecules were used and the ratio between urea and water molecules was adjusted to reproduce the range of concentrations of interest. According to the Metropolis algorithm, new configurations generated in the simulation procedure should be obtained through translations and rotations of a randomly chosen molecule in the simulation box. The ranges permitted for the translations and rotations for the chosen molecule were set to provide an acceptance of about 40-60% of all the motions attempted. For urea molecules, the ranges for translations and rotations were set to 0.011 nm and 11 degrees, respectively. For water molecules the corresponding values were set to 0.015 nm and 15 degrees.

After every 1000 configurations, one volume change of the box was attempted and the range of volume change was set to 0,180 nm³ to provide the same ratio of acceptance for the new generated configurations. The interactions between any pair of sites were calculated within a cut-off radius of 1.1 nm. No long-range corrections to the Lennard-Jones potential were made beyond the cut-off radius, provided that its contribution to the total intermolecular energy was less than 2%.²³ The simulation boxes were equilibrated by the generation of 1´10⁶ configurations and for the averaging procedure, 3´10⁶ additional configurations were generated.

The error bars of the equilibrium averages are calculated according to the expressions:²⁴

This approach assumes that A(t) is a Gaussian function and that <A>_run is the average of property A evaluated in simulations of t_run configurations in the Monte Carlo method.

Results and Discussion

Simulation with OPLS urea model

The Monte Carlo simulations were carried out in the following concentrations: 4, 5, 6, 7 and 8 mol L^-1. This range was chosen because an interesting behavior of the experimental density of the urea solutions was noted as a function of the concentration and is illustrated in Figure 1, where the first derivative is also plotted. The inflexion point near 6 mol L^-1 may be interpreted as a slight expansion of the whole solution as the concentration is raised from 0 to 6 mol L^-1. After this point the solutions appear to contract, which suggests that the solution is less densely packed at the concentration of 6 mol L^-1. This phenomenon may be related to other important properties of aqueous urea solutions.

The Lennard-Jones intermolecular potential parameters of the OPLS urea model, described in Table 1, were obtained from the works of Duffy et al.¹⁷ and Jorgensen and Swenson,²⁵ and the geometric data came from a crystal electron diffraction study.¹⁸

The results for the calculated densities and heats of vaporization (D_vapH) of the studied solutions are summarized in Table 2 and are in good agreement with the experimental values. As can be noted, errors are less than 0.4% for the densities and around 4% for D_vapH. The experimental vaporization enthalpies have been obtained from vapor pressure data²⁵ and the densities of urea solutions at 298 K were taken from reference 26.

Although the OLPS urea model was optimized to reproduce the experimental densities and heats of vaporization of dilute urea solutions, these properties were well reproduced for urea solutions at high concentrations as well as was the behavior of the derivative of the density plot (not shown).

The average interaction energies of urea-urea (E_uu), urea-water (E_uw), water-water (E_ww) and the total energy are plotted in Figure 2. This energy partition is a consequence of the adopted pair additivity approximation, which is represented in equation 6:

The energies E_uu, E_uw, and E_ww show a similar behavior to the interaction energies in a quasi-ideal mixture, that change almost linearly with the molar fraction.

The contribution of E_uu to the total energy (E_T) is very small in the range of concentrations studied, when compared to E_uw and E_ww. E_uu represents only about 3% of the total energy in the least concentrated solution and its contribution to E_T increases almost linearly until the most concentrated one. The interaction energy between water molecules is lowered as the solution becomes more concentrated. In the most concentrated solution we observe that E_uw is almost equal to E_ww suggesting that urea molecules interact with water in almost the same manner as water molecules with other water molecules.

The possible influence of urea on the local water structure can be verified by comparing the radial distribution function (rdf) of pure water to that of the water in urea solutions. The comparison between g_OwOw (r) and g_OwHw (r) is shown in Figure 3. The first peaks in the plots of g_OwOw (r) in Figure 3a show a strong correlation between oxygen atoms near 0.28 nm. Its shape and position are practically unaffected as the urea concentration is changed. The variation of the concentration of urea only causes a small change in the intensity of the peaks of g(r).

The radial distribution functions are defined by equation 7 that shows an inverse dependence of g(r) with r_j, the numerical density of sites j.²⁴

When urea is added to obtain more concentrated solutions, r_j is slightly lowered resulting in an increase of the intensities of peaks in the rdf of water.

The integration of g_OwOw(r) from 0 to the first minimum of the non-normalized g(r) functions, shown in Figure 3a, results in the average number of water molecules in its first hydration shell. For pure water it is almost equal to 5 and for the urea + water mixtures, this number is 4, as is shown in Table 3.

The second peak of the rdf plots, shown in Figure 3a, is correlated with the second hydration shell of water. Only a limited change in the position of this peak was observed, which suggests that urea addition in this concentration range does not disturb the tetrahedral structure of water.

Figure 3b shows the g_OwHw(r) rdf, which is correlated with the hydrogen bonds between water molecules. As observed in the previous case, the peaks maintain their shape and positions and are almost invariant from pure water to the most concentrated solution. However, the integration of g_OwHw(r) shown in Table 3 indicates significant reduction of the number of hydrogen bonds between water molecules, in agreement with experimental data.⁸

The results of g_OwHw(r) rdf and the integration of the first peak indicate that the addition of urea does not cause a significant perturbation in the water structure, although 50% of the water hydrogen bonds were broken. This conclusion, that the water structure in not disturbed when half of the hydrogen bonds are broken, is apparently contradictory. However, Hoccart and Turrell¹¹ suggested that the water structure could be maintained by dipole-dipole interactions.

Figure 4 shows the radial distribution functions g_CuCu(r) and g_OuHu(r) for urea molecules. The first peak of the g_OuHu(r) plots, shown in Figure 4b, indicates that hydrogen bonds between urea molecules are centered at about 0.20 nm and the first peak of the g_CuCu(r), in Figure 4a, shows that the average distance between these sites (located near the center of mass of the molecule and correlated with the distribution of urea molecules around a reference urea molecule) is about 0.45 nm and does not change with concentration.

The integration of the first peaks in Figures 4b and 4a gives the average number of hydrogen bonds between urea molecules and the average number of urea molecules around a reference urea molecule, as given in Table 3. The increase of the average number of hydrogen bonds and the average number of neighbors can be explained by the tendency of urea to form dimers as the concentration is raised. The energetics of urea-urea interactions, shown in Figure 2, support the conclusions based on the radial distribution functions. This result is in agreement with the experimental observations of Stokes⁴ and other theoretical investigations,^15,16 but it is contrary to other experimental studies.^10,11

Figures 5a and 5b show the g_OuHw(r) and g_HuOw(r) radial distribution functions between urea and water sites. The results of the integration (n_OuHw and n_HuOw) of these rdf (Table 3) show that urea forms hydrogen bonds with water and the average number of these bonds is about 5 in the most dilute solutions, in good agreement with the work of Duffy et al.¹⁷ Although the first peak of the g_HuOw(r) rdf is less than unity, the integration shows that there are hydrogen bonds in this site of urea.

The average distance between the centers of mass of urea and water, indicated by the g_CuOw(r) rdf in Figure 6, is about 0.4 nm . The corresponding peak of the 6 mol L^-1 solution is broader than the others, indicating that, at this concentration, the water molecules are, on the average, more distant from urea molecules, characterizing a less densely packed solution, in very good agreement with the behavior of the experimental and calculated densities shown in Figure 1.

Simulation with non-planar urea model

In this section, we studied a proposed model for urea in which the hydrogen atoms are out of the plane of the heavy atoms, with dihedral angles of cis and trans hydrogen atoms of about 13 and 150 degrees, respectively.¹⁹ The same concentrations of the previous urea model were used, as well as the Monte Carlo simulation conditions. The intermolecular Lennard-Jones potential parameters were transferred from the planar urea model and the charges were derived from electrostatic potential grid²⁰ (CHELPG) calculations using the software Gaussian 94.²⁹ The non-planar urea model, described in Table 1, has a lower dipole moment (3.2 D) than either the OPLS model (4.9 D) or the experimental value (5.7 D).³⁰

The calculated densities and vaporization enthalpies for non-planar urea solutions are shown in Table 2. The results for densities are very good, but we obtained some large deviations in the calculations of the enthalpy of vaporization. In this non-planar model, the less densely packed effect of the 6 mol L^-1 solution was not observed.

The partition of the total configurational energy was also analyzed in this model of urea by the pair-wise additivity approximation. The total average configurational energy (E_T) and the contribution of the individual pair interaction energies of urea-urea (E_uu), urea-water (E_uw) and water-water (E_ww) to E_T show a similar behavior of these partial interaction energies as does the OPLS urea model. Thus, non-planar urea aqueous solutions are also quasi-ideal solutions.

The similarity observed for the interaction energies between the OPLS and the non-planar models was also observed for the g_OwOw(r) and g_OwHw(r) rdf, for the average number of water molecules in its first hydration shell and for the same significant reduction of the number of hydrogen bonds between water molecules.

Figures 7a and 7b present the g_CuCu(r) and g_OuHu(r) rdf for non-planar urea sites, where the most significant differences between the non-planar and OPLS urea models are noted . In the non-planar urea model, there is essentially no interaction between urea molecules. The integration of the first peak of g_OuHu(r) rdf shows that the average number of hydrogen bonds (n_OuHu) is 0.1 in the concentration range studied.

The first peak of g_CuCu(r) shows some splitting that is not to be expected and may be due to the small basis set used to optimize the intermolecular potential.

The integration of the first peak of the radial distribution functions g_OuHw(r) and g_HuOw(r) results in 3 hydrogen bonds between urea and water while there are 5 hydrogen bonds in the OPLS urea model. The difference of the dipole moment between the two urea models may be the cause for the reduction in the number of hydrogen bonds.

Conclusions

Thermodynamic and structural results from liquid state simulations showed that urea has a tendency to form dimers when the concentration of the solutions is above 6 mol L^-1.Moreover, its aqueous solutions exhibit a quasi-ideal solution behavior in the range of concentrations studied. No significant perturbations of the water structure were evident from the data of the radial distribution functions, even in a mixture with a concentration of 8 mol L^-1. The number of hydrogen bonds between water molecules was observed to decrease with the addition of urea; however, the whole structure of water was kept unchanged. Most of the results of the simulations for the non-planar urea solutions were similar to those of the OPLS urea model simulations. There is a difference in the description of the interaction between two urea molecules that was not well described by the intermolecular potential used, which is presumably due to the basis set chosen to minimize the structure and to calculate the charge distributions over the atoms.

The results of the simulation with the two urea molecule models show that urea does not disturb the structure of water, even in a highly concentrated solution. Urea appears to enter the water structure as if it was a large water molecule.

Acknowledgements

The authors are grateful for financial support from CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico ) and FAPEMIG (Fundação de Amparo à Pesquisa do Estado de Minas Gerais).

Received: January 3, 2001

Published on the web: March 18, 2002

FAPESP helped in meeting the publication costs of this article.

• 1. Brandts, J. F.; Hunts, L. J.; J. Am. Chem. Soc 1967, 89, 4826.
• 2. Wetlaufer, D. B.; Malik, S. K.; Stoller, L.; Coffin, R. I.; J. Am. Chem. Soc 1964, 86, 508.
• 3. Shick, M. J.; J. Phys. Chem 1964, 68, 3585.
• 4. Stokes, R. H.; Aust. J. Chem 1967, 20, 2087.
• 5. Bertran, C. A.; Macedo, M. I. F.; Oliveira, M. A. S.; Thim, G. P.; J. Eur. Ceramic Soc. 2001, 21, 759.
• 6. Bertran, C. A.; Thim, G. P.; Lima, P. T.; An. Ass. Bras. Quim 1996, 45, 142.
• 7. Lima, P. T.;Bertran, C. A.; Quim. Nova 1998, 21, 608.
• 8. Frank, H. S.; Franks, F. J.; J. Chem. Phys 1968, 48, 4768.
• 9. Holtzen, A.; Emerson, M. F.; J. Phys. Chem 1969, 73, 26.
• 10. Hoccart, X.; Turrell, G.; J. Chem. Phys 1993, 99, 8498.
• 11. Hoccart, X.; Turrell, G.; J. Mol. Struct 1995, 349, 41.
• 12. Kuharski, R. A.; Rossky, P. J.; J. Am. Chem. Soc 1984, 106, 5786.
• 13. Nakanishi, K.; Chem. Soc. Rev 1993, 177.
• 14. Tanaka, H.; Touhara, H.; Nakanishi, K.; Watanabe, N.; J. Chem. Phys 1984, 80, 5170.
• 15. Tanaka, H.; Nakanishi, K.; Touhara, H.; J. Chem. Phys 1984, 82, 5184;
• 16. Hernández-Cobos, J.; Ortega-Blake, I.; Bonilla-Marín, M.; Moreno-Bello, M.; J. Chem. Phys. 1993, 99, 9122.
• 17. Åstrand, P. O.; Wallqvist, A.; Karlstrom, G.; Linse, P.; J. Chem. Phys 1991, 95, 8419.
• 18. Duffy, E. M.; Severance, D. L.; Jorgensen, W. L.; Isr. J. Chem 1993, 33, 323.
• 19. Swaminathan, S.; Craven, B. M.; McMullan, R. K.; Acta Cryst. B 1984, B40, 300.
• 20. Cirino, J. J. V.: PhD Thesis, Universidade Estadual de Campinas, Brazil, 1999.
• 21. Brenneman, C. M.; Wiberg, K. B.; J. Comput. Chem 1990, 11, 361.
• 22. Jorgensen, W. L.; J. Chem. Phys 1979, 71, 5034.
• 23. Freitas, L. C. G.; Programa Diadorim, Universidade Federal de São Carlos, Brazil, 1992.
• 24. Freitas, L. C. G.; Cordeiro, J. M. M.; Garbujo, F. L. L.; J. Mol. Liq. 1999, 79,1.
• 25. Allen, M. P.; Tildesley, D. J.; Computer Simulation of Liquids, Oxford: London, 1987, p. 191.
• 26. Jorgensen, W. L.; Swenson, C. J.; J. Am. Chem. Soc 1985, 107, 569.
• 27. Perman, E. P.; Lovett, T.; Trans. Faraday Soc 1926, 22, 1.
• 28. Egan Jr., E. P.; Luff, B. B.; J. Chem. Eng Data 1966, 11, 192.
• 29. West, R. C., ed.; Handbook of Chemistry and Physics, 58^th ed., CRC Press: New York, 1978.
• 30. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, T. K.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Cioslowski, J.; Stefanov, B. B.; Nanayakkara, A.; Challacombe, M.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Head-Gordon, M.; Gonzalez, C.; Pople, J. A.; Gaussian 94, Revision D.2, Gaussian, Inc, 1995.
• 31. Gaumann, T.; Helv. Chim. Acta 1958, 41B, 1965.

Publication Dates

• Publication in this collection
  25 June 2002
• Date of issue
  2002

History

• Received
  03 Jan 2001
• Accepted
  18 Mar 2002