{'fraction_non_alphanumeric': 0.034356410035347874, 'fraction_numerical': 0.010992326924734891, 'mean_word_length': 4.576682692307692, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Three-dimensional"Mercedes-Benz"model for water

In this paper we introduce a three-dimensional version of the Mercedes-Benz model to describe water molecules. In this model van der Waals interactions and hydrogen bonds are given explicitly through a Lennard-Jones potential and a Gaussian orientation-dependent terms, respectively. At low temperature the model freezes forming Ice-I and it reproduces the main peaks of the experimental radial distribution function of water. In addition to these structural properties, the model also captures the thermodynamical anomalies of water: the anomalous density profile, the negative thermal expansivity, the large heat capacity and the minimum in the isothermal compressibility.


I. INTRODUCTION
Water is the most important fluid on earth. It covers two thirds of the planet's surface and controls its climate. Most importantly, water is necessary for carbonbased organic life being the solvent in most in vivo chemical reactions. Its unique hydration properties drive biological macromolecules towards their three-dimensional structure, thus accounting for their function in living organisms [1]. Water exhibits anomalous properties that affect life at a larger scale. For example, mammals benefit from the large latent heat of water to cool them down through sweating, while water's large heat capacity prevents local temperature fluctuations, facilitating thermal regulation of organisms.
These anomalous properties result from a competition between isotropic van der Waals interactions and highly directional hydrogen bonding (H-bond). A large number of models of varying complexity have been developed and analyzed to model water's extraordinary properties, for reviews see e.g. Refs. [2,3,4,5], but none of the current models can correctly reproduce all physical properties of water. Those model are typically calibrated against experimental data, for example the radial distribution function (RDF) at ambient conditions [6,7], or the temperature of maximum density [8], i.e. T = 3.98 • C. While there is no guarantee that a model optimized to reproduce a given property is able to account for others, adding details increases its quantitative accuracy. For example, TIP5P, which describes water through 5 interacting sites, is typically more accurate [8] than models with 3 or 4 interacting sites. The addition of each interacting site, however, makes the model considerably much more demanding computationally. This is an undesirable feature since a large number of water molecules is required to hydrate even the smallest peptides resulting in a high computational cost. Thus, simple models, such as SPC [9] and TIP3P [10,11], are the most used ones in computational studies of biologically motivated systems. In addition, new models and improvements appear frequently in literature, see e.g., Refs. [12,13,14,15] and references in them.
Coarse-grained models have also been developed and used to study the emergence of water's anomalous properties from its atomic constituents. Both lattice [16,17,18] and continuous models [19,20,21,22] have been applied. Current coarse-grained models cannot, however, be easily used to study hydration of macromolecules, because they can not reproduce the structure of liquid water which is essential in studies of biological systems and molecules [1]. A proper structural description is required since hydration and, in particular, the hydrophobic effect, which is the main driving force for protein folding [23,24], depend on the amount of structural order close to the hydrated molecule versus the amount of order in bulk water. A simplified model that would account for both thermodynamical and structural properties of water, would therefore be highly beneficial in studies related to the hydrophobic effect, protein folding and macromolecules in general.
The main purpose of this work is to introduce a simple but realistic model that reproduces both the main structural and thermodynamic properties of water. To this end, we extend the two dimensional (2D) Mercedes-Benz (MB) model [25] to 3D. In 2D, the MB model has already provided insights into several properties of water: its anomalous thermodynamical behavior [19], hydration of non-polar solutes [26], ion solvation [27], cold denaturation of proteins [28], and the properties of different amino acids [29]. Despite this success, there are several mechanisms which cannot be studied in two dimensions and an extension to 3D is needed.
We show that a previously proposed framework for the 3D MB model [30] does not reproduce the thermodynamical anomalies of water. Here, we extend the model to overcome that problem by making H-bonding dependent on the local environment of atoms by penalizing compact configurations in favor of open-packed ones. With this implementation, structural and thermodynamical properties of are recovered qualitatively. We would like to emphasize that our goal is not to replace atomistically accurate models, such as TIP5P, but rather to provide an alternative for qualitative studies involving water.
The rest of this paper is organized as follows: next, we introduce and discuss the existing 3D MB model, and propose a correction that makes the model suitable to describe thermodynamical and structural properties of water. In the same section, we present the Monte-Carlo scheme and the cooling protocol used in this work. We present our results and comparison to experimental results in Sec. III. In the section entitled results, experimental data are compared qualitatively to our simulations. Finally, we present our conclusions and a discussion in Sec. IV.

A. Mercedes Benz model
In the 3D MB model, water molecules interact explicitly through two types of empirical potentials: H-bonds and van der Waals. H-bonds are directional and account for the tetrahedral structure of water which is described by four arms separated from each other by angles of 109.47 • , see Fig. 1. The energy of H-bonds is minimized whenever arms of adjacent molecules point towards each other. Mathematically if X i represents the position of the i th particle and its four unitary arms, which are denoted by i k (with k = 1, 2, 3, 4), then H-bond interaction between molecules i and j can be written as where and G(x, σ) is an unnormalized Gaussian function The above mathematical description ensures that the intensity of a H-bond is maximized whenever the arms of neighboring molecules are aligned with the vector r ij joining their centers of mass and whenever their distance is equal to R HB . The spherically symmetric van der Waals interactions are approximated by a Lennard-Jones potential: where LJ describes the strength of the interaction and σ LJ is the particle diameter. Then, the total energy describing two MB particles is given by Bizjak et al. [30] studied this model using the following set of parameters: HB = −1, LJ = 1/35 HB , R HB = 1, σ LJ = 0.7, σ R = σ θ = 0.085. They assumed a diamond structure for the model's ground state which is the configuration taken by oxygen atoms when water forms cubicice, i.e. Ice-I c . When tested against simulation, however, this assumption fails and the model can be shown to minimize its energy in an Ice-VII configuration: two interpenetrating diamond lattices with no H-bonds connecting these lattices.
Ice-VII appears to be the optimized ground state for systems trying to maximize their density within a tetrahedral symmetry. It is therefore natural that the 3D MB model of Bizjak et al. [30], whose H-bond term imposes a tetrahedral configuration and the van der Waals term favors compact conformations, has this structure as its ground state. However Ice-VII is not the desired ground state for models of water at ambient pressure such that 3D lattice models for this material have an explicit energetic term penalizing compact configurations of this type [18,31,32]. Figure 2(a) shows a typical Ice-VII ground state obtained by quenching a system of 256 MB particles interacting through the framework of Bizjak et al.. Details about the simulation method and the cooling procedure will be described later in this section. Figure 2(b) shows the density obtained along quenching. The liquid phase has a higher density than ice -as in real water. However, the model does not reproduce the density anomaly of water: the liquid phase does not show a temperature of maximum density below which the density decreases. As a result, the thermal expansion coefficient is never found to be negative in the liquid phase ( Fig. 2(b)). In addition, the model does not reproduce the structure of liquid water (see Fig. 3 in Ref. [30]): the simulated RDF has a non-realistic peak at a distance corresponding to the van der Waals radius.
Despite the problems cited above, the 3D MB model remains an attractive coarse-grained model for water. It does not require calculation of charges, which enables longer simulation times desperately needed in studies of macromolecules. It also holds the promise of being able to provide a qualitatively accurate description of the structure of water due to its tetrahedral nature [33] and water's thermodynamical properties since the model exhibits both open and close packed structures required to described water's anomalous behavior [20]. Next, we de-scribe how the model of Bizjak et al. [30] can be improved to better describe the structure and thermodynamics of water.

B. Corrections to the Mercedes-Benz model
To resolve the above problems, we introduce a term that depends on the local environment of particles. Our approach is inspired by Tersoff-like potentials for covalent materials [34]. This term penalizes H-bonds which are formed in crowded environments through the factor where z i is the coordination of atom i, computed as with the cut-off function defined by [34]: where R and D are chosen as to include the first-neighbor shell only. Note that f (r) decreases continuously from 1 to 0 in the range R − D < r < R + D.
The energy of H-bonds corrected through Eq. (6) becomes: This equation penalizes H-bonds whenever interacting molecules have more than four neighbors. This inhibits the formation of compact tetrahedral phases, e.g. Ice-VII, and favors open-packed tetrahedral phases such as Ice-I. In order to ensure that H-bonds favor chair-like configurations required for diamond structure, we also add a standard potential with three-fold symmetry for dihedral angles. This potential adds an energetic cost to the H-bond between arms m of molecule i and arm n of molecules j, if the dihedral angle formed by the other arms of these molecules is not 60 • : where φ is the strength of the interaction. The term U mn HB (r ij , i m , j n )b(z i ) ensures that the penalty is proportional to the strength of the H-bond. The dihedral angle φ kl describes how the arm k of molecule i aligns with the arm l of molecule j along the vector joining the center of mass of these two molecules. Thus, the total dihedral energy between molecules i and j is Note that because of the dependence on the lo- . This asymmetry has no physical implications since U c HB and U φ possess all the invariance properties required for a potential [34].
We can now write the total potential energy between two water molecules as This model has 10 parameters which were chosen such as to account for a semi-quantitative agreement of the density profile with experiment. We proceeded in two steps to adjust these parameters. First, we chose the values for these parameters such as to produce a density in g/cm 3 [37] that is comparable to experimental values, i.e. about 1 g/cm 3 for the liquid phase and 0.93 g/cm 3 for the ice phase. Only under this condition can the structure of the model be qualitatively similar to real water. Then, we adjusted the parameters such as to obtain a density that is a concave function of temperature with its maximum close to the freezing point. This second condition is the minimal requirement for describing the anomalous properties of water.
The set of parameters calibrated according to the above procedure is given here in reduced units. We report energies and distances in terms of the binding energy | HB | and equilibrium distance R HB of the H-bond. In these units, the three binding energies describing the system are HB = −1, LJ = 0.05 and φ = 0.01. The two distances are R HB = 1 and σ LJ = 1.04/2 1/6 . The two terms controlling H-bond interaction are σ R = 0.1 and σ θ = 0.08, and the three parameters controlling the penalty of crowded environments are υ = 0.5, R = 1.3 and D = 0.2. In this work, temperature is given in units of | HB |/k B , where Boltzmann's constant k B is set to unity. Pressure is given in units of | HB |/R 3 HB . While adjusting the parameters, we found that the behavior of the system is robust upon changing the variables characterizing crowded environments. It is, however, sensitive to the ratio between the binding energy of the van der Waals interaction and the binding energy of the H-bond. This ratio controls the interplay of forces leading to an environment where MB molecules are radially surrounded by their first-neighbors, and forces favoring a tetrahedral distribution of the first-neighbors. The latter favors a high density configuration while the former accounts for a low density one. As opposed to the 2D MB model, we kept the equilibrium distance of the van der Waals interaction comparable to the equilibrium distance of the H-bond such as to avoid artificial peaks in the RDF [30].

C. Numerical simulation method
For numerical simulations, we use the isothermalisobaric (NPT) ensemble to study the thermodynamical properties of a system made of N = 256 MB particles. A Monte-Carlo scheme is used where, at each step, an attempt is made to displace the center of mass and the orientation of particles randomly by a quantity ∆R max and 0.125 rad, respectively. The maximum translational displacement is chosen such as to give an acceptance ratio of 50%. Periodic boundary conditions are used to mimic an infinite system and at every 5 Monte Carlo sweeps, an attempt to rescale the size of the box is made (1 Monte Carlo sweep is equivalent to N attempted steps).
To obtain thermodynamical data throughout the desired range of temperatures, the initial configuration of the system is chosen randomly and equilibrated at the highest temperature (T = 0.17) for 5 × 10 4 sweeps, after which statistics are gathered for the same amount of time. Then, the system is cooled down by ∆T = 0.002 and a similar cycle of equilibration/data gathering is performed. This cooling procedure is repeated until the lowest temperature, i.e. T = 0.11 is reached. At the transition temperature an additional cycle of equilibration/statistics gathering ensured that the system was equilibrated properly. For all the pressures studied here, this protocol was repeated for 10 samples differing by the initial condition. All the quantities reported are the average over those 10 samples and, whenever relevant, the root-mean-square of this average is also shown as the error-bar.
The quantities computed during the simulations were the average potential energy per particle E, the volume per particle V , the heat capacity C P , the compressibility κ T and the thermal expansion coefficient α P . The last three quantities are computed mathematically from the standard fluctuation relations: where H corresponds to the enthalpy of the system. As for the other quantities computed during the simulation, these response functions will be given in reduced units. Thus, C * P will be reported in dimensionless units, κ * T in terms of R 3 HB / HB and α * P in units of −1 HB .

III. RESULTS
In Fig. 3, we provide a qualitative comparison between the properties of bulk water (left panels) and the MB model at P = 0.2 (right panels). The behavior of the MB model follows the trends of water quite accurately: the anomalous density profile (panels on the first row), the negative thermal expansivity (second row), the minimum in the isothermal compressibility (third row) and the large heat capacity (fourth row).
At ambient pressure, water freezes into an open packed configuration called hexagonal-ice, i.e. Ice-I h . This structure is held together by H-bonds which break when ice melts. At this transition, water molecules fill part of the empty spaces, assuming a higher density. In the liquid phase close to the melting temperature, a few openpacked configurations persist -held together by H-bonds. As the system is heated up, those bonds melt gradually removing empty spaces and increasing the density of the system. This reduction of empty spaces occurs until the temperature of maximum density is reached. At this point thermal fluctuations decrease the density of the system with increasing temperature. This behavior has been measured experimentally (Fig. 3(a)) and is captured by the MB model ( Fig. 3(b)): abrupt increase of the density at the melting transition and concave temperature dependence for the density of water with a maximum close to the melting transition.
The thermal expansion coefficient is proportional to the derivative of the volume with respect to temperature α P = 1/V (∂V /∂T ) P . As for most materials, α P decreases upon cooling (Fig. 3(c)) -indicating that the volume of water decreases with temperature. It becomes zero at the temperature of maximum density and negative close to the freezing point. This unusual negative expansivity is reproduced in the model (Fig. 3(d)) and reflects the unusual behavior of water to expand upon cooling below the temperature of maximum density.
The isothermal compressibility measures the tendency of a system to change its volume when the applied pressure is varied: κ T = −1/V (∂V /∂P ) T . For a typical material, κ T decreases upon cooling since it is related to density fluctuations whose amplitude becomes smaller as temperature decreases. This is in contrast with the behavior of water (Fig. 3(e)). The compressibility of water is a convex function of temperature and has a minimum. This anomalous behavior can be explained by noticing that the compressibility is lower for highly packed system than for loosely packed ones since highly packed systems are less susceptible to rearrange their conformation when subjected to a pressure change. Thus, κ T correlates with the volume of the system [19]. Now, since the volume of water is a convex function of temperature, κ T is also convex with respect to temperature for water -see Fig. 3(e). Fig. 3(f) shows that the simulated compressibility is also a convex function of temperature, although the curvature is not very pronounced and its minimum is not as pronounced as in experiments.
Heat capacity, which measures the capacity of a system to store thermal energy (C P = (dH/dT ) P ), is much higher in water than in ice -see Fig. 3(h). This has been explained by the multiple energy storage mechanisms of water as the breakage of van der Waals interactions and H-bonds. The heat capacity of the model presents a much higher variability than real water: close to the transition, C P is much higher than ice and this quantity decreases fast, reaching the same value as ice at about T = 0.15.
In Fig. 4(a), we illustrate schematically the coexistence lines of the solid, liquid and vapor phases of a simple material. At any point along those lines, the free energies of the adjacent phases are equivalent and the Clausius-Clapeyron equation is obtained by equating them: Since ∆h < 0 and ∆v < 0 for the liquid to solid transition of typical materials, (dP/dT ) coex is positive. As a result of this positive slope, a typical liquid freezes when pressure is applied to it -as illustrated by the arrow on Fig. 4(a). On the other side, since water expands upon freezing, ∆v > 0 while the enthalpy difference remains negative (ice has a lower enthalpy compared to water). Thus the coexistence line of the liquid-solid transition has a negative slope, i.e. (dP/dT ) coex < 0. This is illustrated in Fig. 4(b) and leads to the melting of ice when pressure is applied to it. In Fig. 4(c) we show that the model a ) reproduces this anomalous behavior of water. The simulated dependence of the density on temperature is shown for different values of pressure. The freezing temperature shifts to lower values as pressure increases, implying (dP/dT ) coex < 0. Figure 5(a) shows the structure obtained by freezing the MB model. In this configuration, the center of a MB molecule occupies the sites of a diamondlike structure and its arms point to its four neighbors. Note that without a penalty term (Eq. (6)) the empty spaces found in Ice-I can be the stage for the formation of another tetrahedral-like lattice. Thus, this term efficiently shifts the energy of those compact configurations and, in particular Ice-VII, in favor of Ice-I. In Fig. 5(b), the experimental RDF [7] is compared to the ones of the model at different temperature. The second peak of the RDF, commonly referred to as tetrahedral peak [36], is a fingerprint of the tetrahedral geometry of water. It occurs at a distance given by the cosine rule, d 2 = 2R 2 HB −2R HB cos(109.4 • ) ≈ 1.6, much smaller than for a simple liquid [36]. When compared to experiment, the model's RDF has slightly less structure but it peaks at the same position as the experiment -indicating that the average structure of the model agrees well with the experiment.

IV. CONCLUSION
In this work, we have constructed a simple but realistic model for water based on the Mercedez-Benz approach [19]. At low temperature the model freezes forming Ice-I and it reproduces the main peaks of the experimental RDF of bulk water. In addition to these structural properties, the model reproduces the density anomaly of water: ice has a lower density than water and the density of water is a concave function of temperature, with a maximum close to the freezing point. Also, the slope of the solid-liquid coexistence curve is also found to be negative, in agreement with experiments.
In the 2D MB model, the H-bond interaction favors environments having 3 first-neighbors at a distance R HB which competes with the van der Waals interaction that favors 6 neighbors at a distance 0.7R HB . This competition is the underlying physics of the model that accounts for the density anomaly of bulk water. In the 3D MB model, H-bond and van der Waals interactions have the same equilibrium distance and the density anomaly re-sults from an energy penalty for crowded environments. Without this penalty, the system would solidify into a compact Ice-VII configuration. With the penalty term, Ice-VII conformations compete with an open packed diamond-like structure. The competition between these interactions is the underlying mechanism that leads to the density anomaly of the system.
The MB model for water is based on local interactions which are much faster to compute than usual models that uses long-range Coulomb forces. We believe that this model will provide new insights into water mechanisms related to molecular hydration. In particular, investigations of the hydrophobic effect are being undertaken with this model.