532
Acta Chim. Slov. 2007, 54, 532–537
Scientific paper
The Three-Dimensional “Mercedes Benz” Model of Water†
Alan Bizjak1, Toma` Urbi~1, Vojko Vlachy1* and Ken A. Dill2
1Faculty of Chemistry and Chemical Technology, University of Ljubljana, 1000 Ljubljana, Slovenia,
2 Department of Pharmaceutical Chemistry, University of California, San Francisco, California, 94143-1204 USA.
* Corresponding author: E-mail: vojko.vlachy@fkkt.uni-lj.si
Received: 12-06-2007
†Dedicated to Prof. Dr. Jo`e [kerjanc on the occasion of his 70th birthday
Abstract
Previously, we and others have explored the “Mercedes Benz” model of water, a simplified 2-dimensional model that captures many of the qualitative features of water. Here, we introduce the 3-dimensional analogue of this model. Water molecules interact through Lennard-Jones attractions and repulsions, and through Gaussian orientation-dependent hydrogen bonding terms. In 3D MB model of water, hydrogen bonding structure is tetrahedral, hence more realistic than in the 2D version. We explore the model using NPT Monte Carlo simulations. Thermodynamic properties such as molar volume, heat capacity, isothermal compressibility, and the thermal expansion coefficient were studied as a function of temperature at fixed dimensionless pressure, P* = 0.12. These results indicate that the model satisfactorily captures the main experimental properties of real water, including the volumetric and thermal anomalies.
Keywords: Simple model of water, Monte Carlo simulation, thermodynamic properties
1. Introduction
Understanding the physical and chemical properties of water and its interaction with ionic and non-polar solutes is a prerequisite for understanding its role in biological processes1–10. Water’s interesting properties are thought to arise from the ability of water to form tetrahe-drally coordinated hydrogen bonds. Among water’s most interesting anomalies are: i) an increase in the volume of water upon freezing, ii) an increased density of liquid water in the range from 0° up to approximately 4° (depending on the pressure), iii) its high heat capacity and surface tension, and iv) its minimum in the compressibility vs. temperature. A good model of water should explain such anomalies.
There have been different approaches to modeling water. One main approach uses atomically detailed simulation models, which aim for realism and therefore include variables describing van der Waals and Coulomb interactions, hydrogen bonding, etc. (reviewed in8). Even though these models are detailed, they have limitations. Some properties differ with different force fields. One
such example has been provided by Soetens and cowork-ers.11 These authors showed that using a polarizable model of water dramatically changes the resulting potential of mean force between two cations in the model water. Also, Hess and Vegt12 presented an extensive study on hydration thermodynamic properties of analogues of 13 amino acid side chains using different force fields. They found that the choice of water model strongly influences the accuracy of the calculated hydration entropy, enthalpy, and heat capacity. Therefore, the inclusion of detail in water models does not necessarily imply accurate physical predictions.
Alternatively, there are simpler models, having less structural detail, often fewer parameters, and typically requiring much less computation. One class of such models has been developed by Nezbeda and coworkers,13–15 who performed computer simulations and applied Wertheim’s associative theories16,17 to calculate thermodynamic properties of the model liquid and its mixture with solutes. Another example is provided by the so–called MercedesBenz (MB) model,18–28 a 2-dimensional version of representation of water. MB model offers the advantages that it
Bizjak et al.:  The Three–Dimensional “Mercedes Benz” Model of Water
Acta Chim. Slov. 2007, 54, 532–537
533
has few parameters, is computationally very efficient, and gives structure-property relationships that are easy to visualize in two dimensions (reviewed in29). Because of the computational efficiency, MB model can explore fluc-tuational properties, such as the heat capacity, and freezing and melting. Such properties tend to be out of reach for all-atom simulations. The MB model has been applied to studies of solvation of nonpolar20 and polar solutes.28 More recently it has been shown by Becker and Collet30 that the MB model captures the main features of hydration of non-charged amino acids.
Nevertheless, even while MB model captures many of the features of water,18–28 it has limitations. For example, it does not capture the properties of alcohol-water mixtures very well,31 in part because such properties are more sensitive to accurate geometries than some others. In the present work, we develop the 3-dimensional analogue of the MB model.
2. The Model
Here, each water molecule is a Lennard-Jones sphere having four arms oriented tetrahedrally. A hydrogen bond between two water molecules is formed when the arm of one molecule is aligned with the arm of an adjacent neighboring water molecule. The angle between the two arms is 109.4° (see Figure 1). The picture shows the molecules with the hydrogen bonding arms oriented in space. Bold arms are before the plane, dashed arms are behind the plane, and the others are in the plane. The in-termolecular axis is denoted by u_i j.
Figure 1: Two molecules of 3D MB water (i and j), separated by the distance rij. Each molecule has four hydrogen bonding arm vectors; _ik and _jl, respectively (k, l = 1, 2, 3, 4).
where rij is the distance between centers of particles i and j, and X_i denotes the vector representing the coordinates and the orientation of the ith particle. The Lennard-Jones part of the potential is given by
(2)
where e LJ and oLJ are the well depth-depth and the contact parameter, respectively. The hydrogen bonding part of the interaction potential is
(3)
k.l    I
where ?_ i is orientational vector of ith molecule and UHk lB describes the interaction between two arms of different molecules
E^fl(rtf,ßi,ßi)-
(4)
u_i j is the unit vector along r_i j and _ik is the unit vector representing the kth arm of the ith particle. In Eq. (4) G(x) is an unnormalized Gaussian function, defined by
(5)
The strongest hydrogen bond occurs when an arm of one particle is co-linear with the arm of another particle and the two arms point in opposing directions. We make no distinction between electron donors and acceptors. Apart from the dimensionality, we otherwise want to keep the 3D MB model as similar as possible to the original 2D model. Hence, the reduced parameters of our 3D model are the same as were used in the original MB model calculations, with the exception of the depth of the Lennard-Jones potential well, which is ?LJ = 1/35?HB, in contrast to the 2D MB model where it was ?LJ = 1/10?HB. This change was needed to maintain the same ratio between strength of LJ and HB potentials in the two different (2D vs 3D) geometries. In total there are five parameters describing the model. Here, we use the values: ?LJ = –1, and ?LJ = 1/35?HB; the width of the Gaussian function ? = 0.085 was chosen to be small enough that a direct H bond is more favorable than a bifurcated H bond. The contact Lennard-Jones distance is 0.7 of the rHB, where we use rHB = 1.
The interaction potential between two water molecules is a sum of a Lennard-Jones potential and a hydrogen bond term,
(1)
3. Monte Carlo Method
To obtain thermodynamic and structural properties of the model water we performed Monte Carlo simulation method in the isothermal-isobaric (NPT) ensemble. At each step, the displacements in the x, y, z coordinates and
Bizjak et al.:  The Three–Dimensional “Mercedes Benz” Model of Water
534
Acta Chim. Slov. 2007, 54, 532–537
in the orientations of the molecules (three Euler angles) were chosen randomly. We used periodic boundary conditions and the minimum image convention to mimic an infinite system of particles. Except for very low temperatures, where we started from a diamond crystal lattice, the starting configurations were selected at random. The diamond-like crystal lattice was formed by the water molecules consisting of two interpenetrating face-centered cubic Bravais lattices, displaced along the body diagonal of the cubic cell by one quarter the length of the diagonal. It can be regarded as a face-centered cubic (FCC) lattice with the two-point basis 0 and (a/4)(x + y + z ). Here a is the length of FCC lattice and x, y, z are the standard unit base vectors.33 Euler angles for first FCC lattice were set to (71/2,71/2,k/2) and for the second FCC lattice (k/2, k/2, 0). Here are some of the simulation details: 1 × 105 moves per particle were needed to equilibrate the system. The statistics were gathered over the next 1 × 106 moves to obtain well converged results. All simulations were performed with N = 500 water molecules, with the exception that when we modeled a crystal lattice, we used 512 particles. To hold the pressure constant, an attempt was made to scale the dimensions of the box after every pass (one pass was equal to N attempted moves). The maximum volume change was calibrated during equilibration simulations. The mechanical properties such as enthalpy and volume were calculated as the statistical averages of these quantities over the course of the simulations34. The heat capacity, Cp*,the isothermal compressibility, k*, and the thermal expansion coefficient, a* are computed from the fluctuations20.
(6)
All our results below are reported in reduced units, normalized to the strength of the optimal hydrogen bond (e.g. T*=kBT/\eHB\, H*=H/\zHB\, V*=V/r3HB, P*=Pr3HB /\zHB\). Similarly, all distances are scaled by the length of an idealized hydrogen bond separation.
4. Results: Properties of 3D MB Model of Water
Figure 2 shows the simulated energy per particle as a function of the temperature. At lower temperatures, the energy is low because of extensive hydrogen bonding of
the waters. Increasing the temperature first leads to melting, which is indicated as a sharp increase of energy. Further increasing of temperature leads to distorted and broken hydrogen bonds, leading to increased energy. This is illustrated by the water-water pair distribution functions calculated in Figure 3. There are two main peaks in this plot. The first peak at 0.7 is a consequence of the Lennard-Jones interaction, and the second one at 1.0 is a result of the hydrogen-bond interaction. At low temperatures, the hydrogen bonding peak is large relative to the van der Waals
Figure 2: The excess energy as a function of temperature at P* = 0.12.
peak. The peak at distance r ~ 1.6, can be understand by applying the cosine rule. Consider a “cluster” of three water molecules connected by hydrogen bonds. The angle between the bonds is 109.4°. The cosine rule then gives: d2 = r2 HB + r%B - 2rHBrHB cos (109.4°), where d is the distance between the centers of the first and the third mole-
Figure 3: Radial distribution functions for two different temperatures at P* = 0.12. Solid line represents radial distribution function at higher temperature (T* = 0.30) and dashed line at lower temperature (T* = 0.20).
Bizjak et al.:  The Three–Dimensional “Mercedes Benz” Model of Water
Acta Chim. Slov. 2007, 54, 532–537
535
cules. Since rHB is set to unity, the cosine rules gives for d the value of 1.63, which is approximately the position of the relevant peak in Figure 3.
An interesting result is shown in Figure 4, which gives the molar volume as a function of the temperature T*. At low temperatures, water is crystalline with a higher molar volume than liquid water, as it should be. If liquid water is heated above its melting temperature, the molar volume first decreases, then ultimately increases as in normal liquids. The temperature at which the molar volume trend changes is at T* ? 0.14. Hence, this result correctly captures the anomaly of the temperature of maximum density in liquid water. The structure and thermodynamics of the ice phase is dominated by hydrogen bonding. The relatively low density of ice is due to the fact that hydrogen bonding is stronger than the van der Waals interactions. Optimal hydrogen bonding is incommensurate with the tighter packing that would be favored by the van der Waals interactions. Ice melts when the thermal energy is sufficient to disrupt and disorder the hydrogen bonds, broadening the distri-
bution of H-bond angles and lengths. Now among this broadened H-bond distribution, the van der Waals interactions favor those conformations of the system that have higher density. Hence liquid water is denser than ice. Heating liquid water continues to further deform hydrogen bonds and increase the density up to the density anomaly temperature. Further increasing the temperature beyond the density anomaly weakens both H bonds and van der Waals interaction, thus reducing the density. As in real water, this high temperature behavior mimics simpler liquids.
Figure 5 shows the predicted dependence of the heat capacity C*p on temperature. Since the heat capacity is defined as Cp = (?H/?T)P, the heat capacity describes the extent to which bonds are broken (increasing H) with increasing temperature. Breaking bonds is an energy storage mechanism. The heat capacity is low in the ice phase because thermal energy at those temperatures is too small to disrupt the H bonds. The reason liquid water has a higher heat capacity than van der Waals liquids have is because water has an additional energy storage mechanism, namely the H bonds, that can also be disrupted by thermal energies.
Figure 4: The molar volume as a function of temperature at P* = 0.12.
Figure 5: The heat capacity as a function of temperature, P* = 0.12.
Figure 6: The compressibility of water as a function of temperature, P* = 0.12.
The compressibility is correlated with the fluid density. As the density increases, the molecules are closely packed and the compressibility decreases. As bonds break with increasing temperature, the fluid density decreases and the compressibility increases (see Figure 6).
Figure 7 shows the thermal expansion coefficient as a function of the temperature. Note that ? = 1/V(?V/?T)P is merely a derivative of the function shown in Figure 4. The model correctly predicts that the thermal expansion coefficient becomes negative at low temperature, and equals zero for the temperature where the model water has a minimum in molar volume.
Bizjak et al.:  The Three–Dimensional “Mercedes Benz” Model of Water
536
Acta Chim. Slov. 2007, 54, 532–537
Figure 7: The thermal expansion of water as a function of temperature, P* = 0.12.
5. Conclusions
We present here the 3-dimensional version of the Mercedes Benz model of water, a model that has previously been studied only in two dimensions. 3D MB water is more realistic than 2D MB water, both in the obvious factor of dimensionality, and also in that hydrogen bonding is tetrahedral in the former. 3D MB model of water proposed here is simple, has few parameters, and is computationally efficient enough that we can study the full temperature dependence of various thermal properties. We find that the model shows many of the key features of real water, including the volumetric and thermal anomalies: its solid is less dense than its liquid, it has a temperature of maximum density in the liquid range, it has a large heat capacity, it has a minimum in the isothermal compressibility, and it has a thermal expansion coefficient that goes negative in the low-temperature liquid range. One advantage of the 3-dimensional model is that it can be quantitatively compared with experimental data. It can be used to explore solvation of nonpolar and ionic solutes in the same way the 2D model has,28 and it may also be applicable for mixtures of alcohols and water, where the 2-dimensional model is problematic.
6. Acknowledgements
This work was supported by the Slovenian Research Agency (Research Program 0103-0201, J1-6653) and by NIH USA Grant GM063592. A.B. was supported through the Young Researcher Program.
7. References
1. D. Eisenberg, W. Kauzmann, in: The structure and properties of water, Oxford University Press, Oxford, 1969.
2. F. Franks, in: Water, a Comprehensive Treatise, Plenum Press, New York, 1972–1980,Vol. 1–7.
3. F. H. Stillinger, Science 1980, 209, 451–457.
4. C. Tanford, in: The hydrophobic effect: formation of micelles and biological membranes, 2nd ed. Wiley, New York, 1980.
5. W. Blokzijl, J. B. F. N. Engberts, Angew. Chem. Int. Ed. Engl. 1993, 32, 1545–1579.
6. G. Robinson, S. B. Zhu, S. Singh, M. Evans, in: Water in Biology, Chemistry and Physics: Experimental Overviews and Computational Methodologies, World Scientific, Singapore, 1996.
7. R. Schmidt, Monatshefte f¨ur Chemie 1993, 132, 1295–1326.
8. B. Guillot, J. Mol. Liquids 2002, 101, 219–260.
9. A. Ben-Naim, Biophys. Chem. 2003, 105, 183–193.
10. L. R. Pratt, Annu. Rev. Phys. Chem. 2002, 53, 409–436.
11. J-C. Soetens, C. Millot, C. Chipot, G. Jansen, J. G. Angyán, B. Maigret, J. Phys. Chem. B 1997, 101, 10910–10917.
12. B. Hess, Nico F. A. van der Vegt, J. Phys. Chem. B 2006, 110, 17616–17626.
13. I. Nezbeda, J. Kolafa, Yu. V. Kalyuzhnyi, Mol. Phys. 1989 68, 143–160.
14. I. Nezbeda, G. A. Iglezias–Silva, Mol. Phys. 1990, 69, 767–774.
15. I. Nezbeda, J. Mol. Liquids 1997, 73/74, 317–336.
16. M. S. Wertheim, J. Stat. Phys. 1986, 42, 459–477.
17. M. S. Wertheim, J. Chem. Phys. 1987, 87, 7323–7331.
18. A. Ben–Naim, J. Chem. Phys. 1971, 54, 3682–3695.
19. G. Andaloro, R. M. Sperandeo-Mineo, Eur. J. Phys. 1990, 11, 275–282.
20. K. A. T. Silverstein, A. D. J. Haymet, K. A. Dill, J. Am. Chem. Soc. 1998, 83, 150–151.
21. K. A. T. Silverstein, K. A. Dill, A. D. J. Haymet, Fluid Phase Equilibria 1998, 120, 3166–3175.
22. N. T. Southall, K. A. Dill, J. Phys. Chem. B 2000, 104, 1326–1331.
23. T. Urbi~, V. Vlachy, Yu. V. Kalyuzhnyi, N. T. Southall, K. A. Dill, J. Chem. Phys. 2000, 112, 2843-2848.
24. T. Urbi~, V. Vlachy, Yu. V. Kalyuzhnyi, N. T. Southall, K. A. Dill, J. Chem. Phys. 2002, 116, 723–729.
25. T. Urbi~, V. Vlachy, Yu. V. Kalyuzhnyi, N. T. Southall, K. A. Dill, J. Chem. Phys. 2003, 118, 5516-5525.
26. T. Urbi~, V. Vlachy, O. Pizio, K. A. Dill, J. Mol. Liquids 2004, 112, 71–80.
27. T. Urbi~, V. Vlachy, K. A. Dill, J. Phys. Chem. B 2006, 110, 4963–4970.
28. B. Hribar, N. T. Southall, V. Vlachy, K. A. Dill, J. Am. Chem. Soc. 2002, 124, 12302–12311.
29. K. A. Dill, T. M. Truskett, V. Vlachy, B. Hribar-Lee, Annu. Rev. Biophys. Biomol. Struct. 2005, 34, 173–199.
30. J. P. Becker, O. Collet, J. Mol. Struc.- Theochem. 2006, 774, 23–28.
Bizjak et al.:  The Three–Dimensional “Mercedes Benz” Model of Water
Acta Chim. Slov. 2007, 54, 532–537
537
31. B. Hribar-Lee, K. A. Dill, Acta Chim. Slov. 2006, 53, 257– 263.
32. A. Ben–Naim, in: Water and Aqueous Solutions, Plenum Press, New York and London, 1974.
33. N. W. Ashcroft, N. D. Mermin, in: Solid state physics, 1976, Thomson Learning, 64–83.
34. D. Frenkel, B. Smit, in: Molecular simulation: From Algorithms to Applications, Academic Press, 2000.
Povzetek
V ~lanku je predstavljena raz{iritev »Mercedes-Benz« (2D) modela vode na tri dimenzionalni model. Molekule vode interagirajo preko Lennard-Jonesovega potenciala in Gaussovih funkcij, ki v modelu predstavljajo vodikove vezi. Slednje so v molekuli razporejene tetraedri~no. S pomo~jo Monte Carlo simulacij pri stalnem tlaku in temperaturi smo izra~unali termodinami~ne lastnosti, kot so prese`na energija, molski volumen, toplotna kapaciteta, izotermna stistljivost in ekspanzijski koeficient. Rezultati ka`ejo, da model dobro opi{e zna~ilne lastnosti vode.
Bizjak et al.:  The Three–Dimensional “Mercedes Benz” Model of Water