The 3-Attractor Water Model: Monte-Carlo Simulations with a New, Effective 2-Body Potential (BMW)

According to the precepts of the 3-attractor (3-A) water model, effective 2-body water potentials should feature as local minima the bifurcated and inverted water dimers in addition to the well-known linear water dimer global minimum. In order to test the 3-A model, a new pair wise effective intermolecular rigid water potential has been designed. The new potential is part of new class of potentials called BMW (Bushuev-Muguet-Water) which is built by modifying existing empirical potentials. This version (BMW v. 0.1) has been designed by modifying the SPC/E empirical water potential. It is a preliminary version well suited for exploratory Monte-Carlo simulations. The shape of the potential energy surface (PES) around each local minima has been approximated with the help of Gaussian functions. Classical Monte Carlo simulations have been carried out for liquid water in the NPT ensemble for a very wide range of state parameters up to the supercritical water regime. Thermodynamic properties are reported. The radial distributions functions (RDFs) have been computed and are compared with the RDFs obtained from Neutron Scattering experimental data. Our preliminary Monte-Carlo simulations show that the seemingly unconventional hypotheses of the 3-A model are most plausible. The simulation has also uncovered a totally new role for 2-fold H-bonds.


Introduction
Liquid water is not only the most abundant liquid of our lithosphere, but it plays a significant role in biological processes. Liquid water unusual properties are still unexplained [Robinson 96]. Recent explorations of the supercooled [Belissent-Funel 98], [Star 99] and supercritical regions [Tassaing 98], , [Tassaing 00] have only increased the awareness of the current lack of theoretical understanding. Strangely enough, the fact that water still constitutes a fascinating and unsolved problem is all but ignored by the general public.
Within the 3-attractor water model ( [Muguet 92], [Muguet 95], [Muguet 97], [Muguet 98], [Muguet 99], [Muguet 00]) we postulate that, in the condensed phase, the water-water dimer interaction features two other local minima corresponding to other types of physical interactions, in addition to well known linear H -bond water interaction. The two extra local minima are the inverted water dimer and the bifurcated water dimer.
Since this postulate relates to the dynamics of the water molecules, we prefer to rely on the concepts of attractors. This allows also to distinguish our model from 2-state state models [Robinson 99], which bring with them a notion of thermodynamical phases, and mixtures, which could be misleading. Stricto sensu, this postulate is the only one required to build our model for the liquid dynamics. A stronger postulate, a non-equivalent hypothesis, is that the isolated water dimer does also feature two extra local minima.

Design of the Potential
One ideal theoretical approach would be to be able perform truly ab initio molecular dynamics or even Monte-Carlo simulations, while solving the electronic Schrödinger equation for the whole simulation box. This would take care of all n-body effects. Then, assuming the Born-Oppenheinmer approximation is valid, a quantum dynamics motion of the nuclei would be performed. This approach would take into account the ZPE, that we have s hown [Muguet 00] cannot be neglected, and all quantum vibrational and even rotational levels. Of course, the method chosen in order to solve the electronic Schrödinger equation would have be accurate enough to reproduce the potential energy surface of a single dimer. Considering the state of the art (that we are aware of), we set forward the challenging view that there is no current ab initio method, which is good enough to reproduce the water dimer PES. One proof i s that current ab initio methods have consistently failed to reproduce the ammonia dimer PES. We have shown that one reason for this failure is the basis set superposition error (BSSE) [Muguet 95]. We proposed a new method to correct the BSSE and it appears that with our new correction method the ammonia dimer PES is improved and in agreement with experimental observations. Yet we believe a convenient and correct method to correct the BSSE has yet to be found. For the ammonia dimer, it is a question of finding the correct global minimum, for the water dimer, it is a question of finding extra local minima, which is much harder and challenging. It must be stressed that performing nuclear classical dynamics, even on a correct n-body ab initio PES would not be sufficient as we have shown the crucial importance of ZPE effects. Some remedies have been proposed to be introduced in a classical MD scheme [Alimi 92]. Another possible approach is to implement path integral representations of the hydrogen atoms, which amounts, in practice, to representing the proton by a necklace. One problem that proves difficult to handle is the quantum correlation between protons and related decorrelation effects [Chatzidimitriou-Dreismann 97]. It is also interesting to notice that the methodical and well-known effort by Clementi and coworkers [Clementi 83] to develop effective water potentials (MCY) with the help of ab initio calculations have not been successful. While waiting to find the correct electronic ab initio PES and then to develop sophisticated MC or MD methods with nuclear quantum dynamics that do not yet exist, a completely different approach is to develop semiempirical effective 2-body potentials. We are proposing a new class of effective potential called BMW potentials (an acronym for Bushuev-Muguet-Water potentials) that is built by modifying existing empirical potentials. This version 0.1 ( BMW v. 0.1) has been designed by modifying the SPC/E empirical water potential. It is a preliminary version suited only for Monte-Carlo simulations.
This alternative empirical approach is less rigorous but it is computationally feasible. We may fit selected parameters of the potential functional while trying to minimize differences between calculated and experimental data. Non-additive interactions may be accounted "effectively", which means in the sense of a statistical average. This "effective" approach is not necessarily very sensitive to structural details since it implies angular averaging. According to statistical physics principles; macroscopic properties are integral quantities of the Hamiltonian, and a rather weak dependence of liquid properties on the details of the effective dimer PES (Potential Energy Surface) is observed. There are many available potentials for liquid water [Robinson 96] and most of them predict a selected set of experimental properties correctly, but none of them appears to be able to predict correct properties for a wide range of state parameters. Some calculated properties are significantly different from experimental data. It must be stressed also that even the functional form of the potential is unknown..
By construction the BMW potentials includes at least two extra local minima in addition to linear dimer (LD) (cf. Fig 1). A first local minima, favored at short intermolecular distances [Muguet 00] is the Inverted Dimer (ID) whose symmetry group is C i with an inversion center at the mid-point between oxygen atoms. At intermediate intermolecular distances, we find the well-known global minimum linear H -bond dimer (LD) whose symmetry group is C s with a mirror plane. At larger intermolecular distances become favored the bifurcated dimer (BD). Sometime the term "bifurcated bond" is used for describing a proton engaged in two hydrogen bonds. We prefer to call this situation a twofold bond, in order to avoid any confusion. The ID and BD generates non-tetrahedral symmetry clusters. Experimental results are pointing out that there is a non-perfect random tetrahedral H -bonded network in liquid water. At a first glance, it could appear that the existence of a large number of ID and BD structures might be contradicted by Xray and neutron scattering data. At near ambient condition radial distribution functions (RDFs) of oxygen atoms has second pick near 4.5 Å that corresponds to the second neighbors from given atom in tetrahedral network. We must, however, remember that by definition where d loc (r) is local number density of atoms in thin spherical layer on the distance r from fixed atom, and d is mean number density of atoms averaged on whole molecular ensemble. When d loc is computed, we average the density both for angles and distances. Different atomic spatial distributions may give the same radial local density profile. The second peak of RDF g oo has much less height than the first one, and therefore we should expect many defects within the tetrahedral network. The nontetrahedral bonds (ID, BD, twofold bonds) might play a significant role and influence all water properties on micro-and macroscopic levels. One of the goals of this work is to check the consequences of the hypotheses of the 3-A model both in terms of thermodynamical and structural quantities.

Potential and simulation details
So far, to our knowledge, current existing empirical water potentials treat BD and ID as stationary points which are transition states but not local m inima. However, we expect [Muguet 00] that MC or MD simulations, even with those potentials, should predict populations of BD and ID simply because the intermolecular Oxygen-Oxygen distance is a temperature and pressure dependent distribution. They should be present, but to our knowledge, nobody has carried out a systematic statistical analysis concerning the presence of ID and BD over a wide range of state parameters. We suspect that the best semi empirical potentials are in fact those who predict, albeit still inaccurately, relatively large populations of ID and BD unknowingly. Consequently we took one of the most popular pair additive potential: the SPC/E potential [Berendsen 87] as the reference potential for further qualitative improvements. The SPC/E potential is able to correctly reproduce vapor pressure, but underestimates the critical temperature of water. On the other hand, the SPC/E model accurately predicts the critical temperature, but underestimate the vapor pressure by more than a factor of two. Our first idea is to approximate the PES near stationary points with the help of additional Gaussian functions. We suppose that the shape of PES in these regions should like a Gaussian-shaped well in a six dimensional space and can be approximated by smooth functions. At this exploratory stage of our investigation, we are not making any further explicit hypotheses concerning the heights of the energy barriers between LD, BD, and ID. A barrier simply results from the intersection of the Gaussian well with the unmodified SPC/E region. We do not use any special functions for barrier approximation.
We underline that this simple procedure might create some discontinuities of the first derivative of the PES. Therefore this preliminary version ( v. 0.1) of the BMW potential is suited for MC simulations but not for MD simulations.
We used rigid molecules with geometrical parameters taken from SPC/E model. The lengths of the O-H bonds were 1 Å and the HOH angle was 109.47°. Universal interactions were described by the Lennard-Jones (12-6) potential where ε= 0.6277 kJ/mol, σ = 3.18 Å. These values are only slightly different from original SPC/E ones (ε= 0.649 kJ/mol, σ = 3.166 Å). Point charges are located on H atoms with a -0.405 e charge, and O atoms with a +0.81 e charge (to compare with -0.4238 e and 0.8476 e respectively for SPC/E).
We used the following Gaussian function for the interaction between two molecules m and n where the adjustable parameters c 1 =0.29092 kJ/mol, c 2 =4.4 Å , and c 3 =0.15 Å 2 , r ij .-are the distances between atoms of two molecules (i=1,2 for atoms H, i=3 for atoms O).
where z i , z j, are unit vectors directed along dipole moments, x i , x j are unit vectors that perpendicular to dipole moments (see Fig. 2 where the other empirical parameter is c 4 =7.1128 kJ/mol.
For ID, optimal geometrical parameters are r 33 = 2.8 Å, r 13 and r 31 = 2.28 Å. Respectively. Then we compute functions s2 and s3 If the product of this functions was greater than 0.01 we regard that ID exist with the energy where c5=8.368 kJ/mol.
The energy dependence on O-O distance for optimal geometry for linear, bifurcated and inverted dimers for BMW and SPC/E potentials are shown in Figure 3.
Monte Carlo simulations for bulk liquid water include 512 molecules in a cubic elementary cell. Periodic boundary conditions are applied in all three directions. The simulations are carried out in an isobaric-isothermal (N,p,T) ensemble in a wide range parameters of state. Standard Metropolis sampling was performed. A potential cutoff of 10 Å is chosen to reduce computational effort. Longrange correction term for the Lennard -Jones interaction was less than 0.3 kJ/mol and was not included in the total energy. We suppose that second order corrections should give a systematic shift on thermodynamic functions and are weakly dependent on state points. New configurations were generated by randomly translating selected monomers in all three Cartesian directions and by random rotations. At 298 K the ranges were 0.15 Å for the translation and 0.15 radians for the rotation. These choices yielded acceptance rates of 45%. Volume moves were made every 1024 accepted configurations and involved scaling all the intermolecular distances. The maximum length of cell change was 0.12 Å that gave acceptance rates of 20%. The lengths of Markov's chain were exceeded 200 millions of configurations.

Results and Discussion
We have chosen the same series of thermodynamical points as those studied in the Neutron Scattering experiments by Soper and coworkers [Soper 97], [Botti 98], [Soper 00]. These papers present different sets of radial distribution functions (RDF). For the most part, however, RDFs are obtained from the same experimental data but analyzed with two different methods. Pair correlation curves from the 1997-1998 and 2000 papers have been established from almost the same set of S(q) experimental data which correspond to the spectra in the momentum space. There is a little more data in the 2000 paper. The transformation from the momentum space to the coordinate space is not as simple, as one would imagine at a first glance, because of truncation. A simple Fourier Transformation does not suffice. We are not going to delve here with this highly technical topic. What matters is that the pair correlation functions were not reconstructed in the same way. For the 1997 set of data, the pair function curves were constructed by a "maximum entropy" fitting method. While for the 2000 set of data, a new EPSR (Empirical Potential Structure Refinement) method was used. In this new method, an "empirical potential" is added to the "reference potential" in order to reproduce the s(q) momentum spectra. The SPC/E empirical potential was chosen as a reference potential. The respective Fourier transform RDFs reconstructions are differing in many respect.
The 1997 curves are more structured than the 2000 ones. One must stress that the SPC/E potential features only the linear water dimer as a mininum and that a significant empirical potential had to be added to the reference SPC/E potential in the EPSR procedure which shows that the reference potential is far from being adequate. We hypothesize that the inadequacy and the bias brought by the reference potential is going to artificially smear out interesting curves features. Extreme caution must be exercised when comparing "experimental" radial distribution functions. For one supercritical condition (Run 673b), it is also possible to compare with the results from the experiments of Bellissent-Funel and coworkers [Tassaing 98], which were obtained in very similar experimental conditions (380°C instead of 400°C). The temperature dependence of the molecular number density and potential energy is compared with the corresponding experimental curves i n Figures 4 and 5 at P=0.1 MPa. The standard deviations are 3·10 -4 Å -3 for number densities and 0.3 kJ/mol for potential energies. Computed water properties agree very well with experimental data. Only small little systematic shifts are observed. The relative error does not exceed 3% for the density and 4.5% for the energy. The ability of the BMW potential in reproducing correctly experimental densities for some state points can be checked in Table 2. It means that calculated thermal expansion coefficient, heat capacity, isothermal compressibility coefficient are quite close to experimental values in a wide temperature range.   We have compared BMW thermodynamic predictions with data obtained from a comparative study [Jedlovsky 99] that used a set of polarizable and non-polarizable potentials. The data are summarized in Table 3. Polarizable potentials are not predicting accurate densities (d) for all state points, especially for 573b and 673b. Flexible potentials based on rigid SPC and TIPS models give errors ~10 -30 % for state points close to liquid-vapor coexisting curve. At ambient conditions, configurational potential energies are close to e xperimental data for all models. Polarizable potentials predict energies that are too high at other state points. Among the non-polarizable potentials, BMW is the most effective in reproducing densities. The BMW potential predicts the best thermodynamical quantities for almost all points.
We hypothesize that polarizable models are exaggerating the strength of the linear H -bond through cooperativity and therefore they severely quenches the populations of ID and BD in an even more radical fashion than conventional non-polarizable models . In the supercritical regime and its vicinity, computed and experimental densities of water with corresponding run labels are reported in Table 4. Significant differences with experimental data are only noticed for the first three points. The maximum deviation is observed for state point 573b, which is very close to the experimental liquid -vapor coexistence curve. It is only on 11 K lower then boiling temperature at 10 MPa. For other states, deviations are minimal Table 4 Density BMW potential a Standard deviations are given in parentheses.
We have used the RDFs reported in the neutron scattering studies NDIS by Soper and coworkers [Soper 97], [Botti 98], [Soper 00]. NDIS-1997 and NDIS-2000 were used for comparison with computed RDFs. The atom-atom RDFs (g ij (r)) for five points are shown in Figures 6, 7 and 8 concerning the low temperature regime (graphics are shifted along the y-axis for clarity). Our simulated results are in much better agreement with the latest RDFs (NDIS-2000) than with older 1997 RDFs. They seems quite satisfactory, but honestly, one should not give too much credit to this good agreement, since often differences between the two "experimental" curves exceed their respective differences with our BMW computed curves. As stressed before, the "experimental" RDFs are not raw experimental data. Neutron Scattering experiments are unusually complex and costly "big science" experiments. The raw experimental a re obtained in the momentum space. A first series of possible uncertainties affect the collection of the data: the reactor initial beam coherence quality, the temperature distribution of the moderator, subtraction of the scattering due to the container vessel and many other technical issues. A second source of errors lies within the extraction of the partial structure factors. It relies on the implicit assumption that the structure of light and heavy water are absolutely identical. A third source of errors lies within the mathematical procedure itself when transforming data from momentum space to spatial coordinates. Fourier transformations are not accurate enough due to the truncation of the momentum space and they are quite sensitive to the data noise. A more detailed analysis of those delicate problems as well as a more trustful direct comparison between BMW computed partial and total structure factors with momentum experimental data will be presented elsewhere.
Within the framework of the 3A-model, some unconventional metastable clusters or chains of ID and BD are expected. Not surprisingly, our BMW potential MC simulations are producing unusual molecular configurations that may be metastable in water and aqueous solutions. These configurations may play a significant role in triggering chemical and physical processes. Needless to say, those clusters might account for native microstructural defects or micro heterogeneities of liquid w ater [Bushuev 97], [Bushuev 99].
A rather striking example of one of those transient clusters is shown in Figure 12. It consists of six molecules and has the form of chain. The sequence of intermolecular interactions is : ( 1-BD-2-ID-3 ID-4-BD-5-ID-6 ).

Figure 12.
A selected water cluster formed by BD and ID interactions.
It is interesting to notice that H atoms belonging to water molecules 3 and 4 (labeled respectively A and B) may also be viewed as forming twofold H -bonds. Proton A is located between the oxygen atoms of water molecules 2 and 4. Similarly proton B lies between the O atoms of water molecules 3 2 3 1 A B 5 6 4 and 5. It is therefore quite possible that twofold H -bonds are providing yet another stabilization interaction to this cluster. In fact, in constrained crystalline structures, occurrence of the 2-fold H-bond is quite common. In the present case, it is the ID/BD structure that provides the constraint. This cluster is likely going to b e the future topic of our ab initio investigations, but only when we are feeling confident enough that proper theoretical methods are available for correcting the BSSE and taking into account the ZPE. Within the context of the 3-A model, two-fold bonds may play a complementary and totally new role by further stabilizing BD/ID clusters. The BMW MC simulation has brought this hindsight. One should also expect that these types of weakly bonded clusters are going to be entropy favored. The non-tetrahedral BD / ID clusters may lead to the formation of a significant population of two-fold 3-body H-bonds.
In connection with high temperature / high pressure liquid water structure, where a significant number of ID are expected, it has been proposed [Bassez 99] that it is the specific structure of supercritical liquid water or high pressure water that triggered prebiotic synthesis reactions. The 3-A model brings therefore some support to the idea that the chemistry in ocean hydrothermal vent systems might constitute the cradle of life.
A detailed (ID, LD, BD) population analysis of the configurations produced in our BMW MC simulations, in comparison to standard SPC/E MC simulations will be presented elsewhere.
Our future studies should concern two fronts: A first issue is to a get a detailed understanding of the water dimer and n-mer PES with ab initio methods. This line of study might require a significant amount of time since new theoretical developments are needed. Without waiting, yet an another assault will be attempted by designing BMW potentials and performing Monte-Carlo and Molecular Dynamics simulations.

Conclusions
According to the precepts of the 3-attractors water model, a new class (BMW of effective 2-body potentials for classical simulation has been proposed. A first version BMW v. 0.1 has been designed starting from the SPC/E potential to be used in Monte-Carlo simulations. Our (N,p,T) Monte-Carlo simulations show that thermodynamical predictions by the BMW potential seem to be the best in our comparisons. Concerning structural predictions in relationship with Neutron Scattering experiments, the BMW potential RDFs predictions appear very satisfactory, but until a comparison can be made between computed and measured structure factors, extreme caution is warranted.
There is no positive proofs in natural sciences, contrary to mathematics. A model may be considered as valid until falsified. Our Monte-Carlo simulations have not highlighted any element allowing one to falsify the 3 -attractor water model. On the contrary, considering the comparative performance of the BMW potential in regards to the SPC/E potential, the 3 -attractor water model appears as a most plausible model. Furthermore, Monte-Carlo simulations have uncovered a totally new role for 2-fold H-bonds.
Despite the preliminary nature of our simulations, results are obviously very encouraging. At the same time, because of the preliminary nature of our simulations, there is much room for improvement, while carrying out a cautious assessment of both theoretical and experimental data