## Química Nova

##
*Print version* ISSN 0100-4042*On-line version* ISSN 1678-7064

### Quím. Nova vol.21 no.6 São Paulo Nov./Dec. 1998

#### http://dx.doi.org/10.1590/S0100-40421998000600004

ARTIGO

**Monte Carlo thermodynamic and structural properties of the TIP4P water model: dependence on the computational conditions **

**João Manuel Marques Cordeiro**

Depto. de Física e Química - Faculdade de Engenharia de Ilha Solteira - UNESP - CP 31 - 15385-000 - Ilha Solteira - SP

Recebido em 7/7/97; aceito em 24/4/98

e-mail: cordeiro@feis.unesp.br

Classical Monte Carlo simulations were carried out on the NPT ensemble at 25^{o}C and 1 atm, aiming to investigate the ability of the TIP4P water model [Jorgensen, Chandrasekhar, Madura, Impey and Klein; J. Chem. Phys., 79 (1983) 926] to reproduce the newest structural picture of liquid water**. **The results were compared with recent neutron diffraction data [Soper; Bruni and Ricci; J. Chem. Phys., 106 (1997) 247]. The influence of the computational conditions on the thermodynamic and structural results obtained with this model was also analyzed. The findings were compared with the original ones from Jorgensen *et al* [above-cited reference plus Mol. Phys., 56 (1985) 1381]. It is notice that the thermodynamic results are dependent on the boundary conditions used, whereas the usual radial distribution functions g(O/O(r)) and g(O/H(r)) do not depend on them.

**Keywords: **Monte Carlo simulation; TIP4P water model; radial distribution function.

INTRODUCTION

It is difficult to over-emphasize the importance of water in chemical and biological processes. The study of its properties has attracted a good deal of attention in spite of the difficulties involved. Experimental results have provided a comprehensive picture of the way water molecules organize themselves in pure water and around particles dissolved in it and supplied a sensitive evaluation of existing theories of aqueous systems^{1}. However there are a lot of obstacles in the theoretical study of the liquid state due to the complexity of the intermolecular forces. It is not feasible in general to calculate statistical mechanics averages by using standard mathematical approaches. With the dissemination of computers, methods such as Monte Carlo (MC) and molecular dynamics have become popular tools in this task^{2-5}. Nevertheless, the success of the investigation is very much dependent on the availability of intermolecular potential functions that yield a reasonable model for the liquids. Consequently, there has been a great deal of interest to develop models that reproduce the water behavior. Among them one can find rigid models^{6-8} and others including polarization effects^{9-11}. In a recent work I have investigated the behavior of some models of water under various thermodynamic conditions and done some considerations about rigid and polarizable models^{12}. As I have noticed, while polarizable potentials surely improve water models, in many ways they are less successful than non-polarization potentials which will continue to be used in water simulations for some time.

In recent years I have undertaken a theoretical investigation of the liquid properties emphasizing the liquid structure and the solute/solvent interaction in aqueous solutions^{13,32}. For water I used the TIP4P model developed by Jorgensen *et al*^{7,14}. This model is, together with the SPC and SPC/E models^{8}, perhaps the most popular model in water simulations. In the quoted papers of Jorgensen *et al*, the authors have shown that the TIP4P model, besides reproducing the thermodynamic properties of water, represents very well the experimental g(O/O(r)) radial distribution function (rdf) obtained from the X-ray diffraction data of Narten and Levy^{15}. The first peak in the g(O/O(r)) function of the TIP4P model is slightly higher than the experimental one, which is a characteristic of several models reported by Jorgensen *et al*^{7}. The g(O/H(r)) and g(H/H(r)) rdfs obtained with the TIP4P model have not been compared with experimental ones. Jorgensen *et al* pointed out that there would be many uncertainties in these radial distribution functions obtained from the X-ray diffraction experiments^{7}. Experimental rdfs for water at 25^{o}C and 1 atm were more recently obtained via neutron diffraction experiments by Soper and Philips^{16}, in which the g(O/O(r)) rdf was estimated by representing the total pair correlation as a Gaussian density distribution for the near-neighbor peak, plus a broader distribution for the rest of the pair correlation function. The peak obtained using this methodology is much higher than that reported by Jorgensen *et al* obtained from the Narten and Levy data^{7}. Since then, some analytical models for water have been proposed by Blum *et al*^{17,18} that represent very well the rdf published by Soper and Philips^{16}. However, in 1994 Soper published a paper in which he used a new formalism to analyze the neutron diffraction data, that does not involve fitting Gaussians, but uses the minimum noise principle to estimate the entire g(r)^{19}. The height of the first peak in the g(O/O(r)) rdf, previously reported as 3.1, decreases in the new data to roughly 2.5. The author declares that in the new work the peak widths and heights have been determined more reliably^{19}. Blum and Degreve do not refer to this work in a recently published paper^{20}. In the present year, Soper and coworkers have made new improvements in the formalism and published site-site pair correlation functions for water from 25 to 400^{o}C^{21}. While in the work of 1994 the three g(r)s (OO, OH and HH) were calculated directly from the measured differential cross sections, each one independently of the others, in the later one, the three rdfs were estimated simultaneously to give the best fit to all three differential cross sections. In this way, systematic errors, which may otherwise cause a positive deviation in one rdf and a negative deviation in another, can be corrected for more accurately^{22}. The g(O/O(r)) plot reported agrees very well with the earlier X-ray results reported by Jorgensen and coworkers^{7}. Thus, the TIP4P model has gained new strenght. By increasing, we have now an accurate experimental g(O/H(r)) to compare with the theoretical one. On the other hand, in many situations it is necessary to investigate systems under computational conditions different from those reported in the original papers of Jorgensen *et al*^{7,14}. Therefore, in the present work we investigated the behavior of the TIP4P model under new circumstances. The structural and thermodynamic results were compared with those reported by Jorgensen *et al*^{7,14} and with experimental ones described in the literature.

DETAILS OF THE SIMULATION

**The Intermolecular Potential Function.**** **As the TIP4P is a rigid model, molecular relaxation effects are not considered. Following usual procedures in force field calculations, the energy E_{ab }between molecules **a** and **b** is represented by a sum of Coulomb and Lennard-Jones potentials centered on the sites, that is

(1) |

where r* _{ij}* is the distance between site

*i*in

**a**and site

*j*in

**b**and q

_{i }and q

_{j}are point charges located respectively on the

*i*and

*j*molecular sites. For each site

*k*, the parameters

*A*and

_{kk}*B*are given by

_{kk}*A*e

_{kk}=4*s*

_{k}*and*

_{k}^{12}*B*e

_{kk = }4*s*

_{k}*, where e*

_{k}^{6}_{ }*and s*

_{k }*are the Lennard-Jones parameters for the*

_{k}*k*-th site.

**Monte Carlo Simulations**** **with Metropolis' importance sampling^{24} and cubic box boundary conditions were carried out on the NPT *ensemble* at 298 K and 1.0 atm, using the DIADORIM code developed by L. C. G. Freitas^{25}. To investigate size effects on the results two samples were simulated, one with 250 and another with 500 monomers in the reference box. Jorgensen and Madura have done the same for other sizes of samples and used a cutoff procedure to calculate the contribution of the Coulombic term to the configurational energy^{14}. In the present work this procedure was not used, but instead the configurational energy was calculated considering all molecules in the box. Also the average was calculated over a significantly larger number of configurations than by those authors. Details of the simulations of the present work and those ones with 216 and 512 molecules in the box, simulated by Jorgensen and Madura^{14}, are recorded in table 1. It is interesting to notice that the cutoff distance used in both cases by Jorgensen and Madura is the same. Starting from an initial distribution of molecules in the central box, a new configuration was generated by randomly translating and rotating a randomly chosen molecule along Cartesian coordinates. In the NPT ensemble, new configurations are also generated by probing the density of the liquid with volume changes. After a volume change, the center of mass coordinates of all molecules in the reference box were scaled in a usual way^{23}. Ranges for translations and rotations of monomers and also volume changes were adjusted to yield an acceptance/trial ratio between 0.4 and 0.45 for new configurations. Statistical uncertainties were calculated from separate averages over blocks of 2x10^{5} configurations. The calculations were carried out on IBM RISC 6000 workstations.

RESULTS AND DISCUSSION

**Thermodynamic Results.**** **The thermodynamic results obtained for the simulations are presented in table 2 compared with those of Jorgensen and Madura^{14} and with experimental data. The origin of the experimental intermolecular energy E_{i} will be discussed below. As one can observe the values obtained in the present work for density, heat of vaporization and intermolecular energy are about 5% greater than those reported by Jorgensen and Madura^{14}. The values reported by these authors are in better agreement with the experimental ones than the results obtained in the present work. Simulations made using the DIADORIM program with cutoff procedure^{26} have shown values in accordance with the results of Jorgensen and Madura^{14}. Then the difference between the present results and those reported by Jorgensen and Madura must be due to the difference in the number of monomers that are interacting and, if this is true, the parameterization of the TIP4P model would be dependent on the cutoff ratio. Previous results obtained by the author indicate that if a liquid is parameterized using a cutoff procedure to calculate the contribution of the Coulombic term to the interaction energy, these parameters do not reproduce the same results if one uses Ewald's sum. This is particularly expected for the TIP4P model for which parameterization was optimized to fit experimental properties using a cutoff procedure^{28}. The experimental intermolecular energy was computed with the equation *DH _{vap }= -E_{i} + P[V(g) - V(l)] = -E_{i }+ RT*, where one assumes that the gas is ideal and that the sum of kinetic and vibrational energies is the same for the liquid and the gas. As discussed by Jorgensen and Madura

^{14}, there is little error in making this approximation. The heat capacity

*C*, the isothermal compressibility

_{p}*k*and the thermal expansion a have been computed from the standard fluctuation formulas

^{23}. The rates of convergence for these properties are slower than for the energy or density, especially

*k*and a

^{29}. The computed configurational contribution to the

*C*has been summed to 3R (the classical kinetic energy contribution of the translation and rotation from the monomers) giving the value shown in table 2. This was the methodology also used by Jorgensen

_{p}*et al*

^{7}. Considering the problems resulting from this approximation and from the way of calculating these properties, the obtained values are in good agreement with the experimental ones. The present results are in better agreement with the experiments than those of Jorgensen and Madura

^{14}and this can be attributed to the number of configurations averaged, which is considerably larger in the present work.

**Radial Distribution Functions.** Figs. 1 and 2 show the g(O/O(r)) and g(O/H(r)) rdf obtained in the present simulations compared with experimental results. The g(O/O(r)) peak obtained by Soper and Philips^{16} is presented for comparison. As one can observe, there is a good agreement between the g(O/O(r)) obtained by Soper *et al*^{21} and that obtained by Jorgensen *et al*^{7} using the Narten and Levy experimental data^{14}. The behavior verified for the TIP4P model in the present work is the same as that verified earlier by Jorgensen *et al*^{7,14}. The height of g(O/O(r)) peak is overestimated compared to the experimental one. The integration of the first peak up to the minimum point gives a value of 5.0 for the O-O coordination number, in agreement with the value reported by Jorgensen *et al*^{7}.

As one can see in figure 2, the TIP4P model also produces an overestimation in the amplitude of the first peak of the O/H correlation compared with the experimental one. This seems to be in accord with the results of the simulations that show that the heat of vaporization obtained with the TIP4P model is greater than the experimental value. Therefore, it looks like the Lennard-Jones parameters are overestimated (at least the one of oxygen). For the second peak the agreement with the experiment is very good. Nothing was found in the literature for comparing the g(O/H(r)) obtained with the TIP4P model with experimental results. By integrating the first peak up to the minimum point a value of 4.0 for the O-H coordination number was obtained. This value agrees with that reported by Jorgensen *et al*^{7} and is a little greater than the experimental one^{21}. Some analyses suggest that the height of the OH peak is not a reliable measure of the degree of hydrogen bonding in water, so that geometric or energetic aspects, or both, should be taken into account^{30-31}. Since the average O-O coordination number per monomer is about 1 unit greater than the average number of hydrogen bonds, one can conclude that one not-hydrogen-bonded molecule is trapped in the first hydration shell due to the intermolecular forces. It would be interesting to try to imagine the effect on the properties of water of this extra molecule bonded in the field of the water molecule. The distortion in the tetrahedral field arising from this molecule must result in a decrease in the melting point of water with clear effects on the development of life. In Fig. 3 we show the g(O/O(r)) and g(O/H(r)) rdfs for water obtained with the TIP4P model for boxes with 250 and 500 molecules. As one can observe, the plots are exactly equal meaning that the g(r)s obtained are not dependent on box size.

CONCLUSIONS

The TIP4P parameterization looks like being dependent on the number of molecules that interact in the simulation box. The parameterization of the TIP4P model was optimized to reproduce the experimental properties of water using a cutoff procedure to calculate the contribution of the Coulombic term to the configurational energy. The model super-estimates the experimental amplitude of the rdf for water as compared with the more recent results. The greater amplitude in the O/O rdf and of the first peak of the O/H rdf is in accord with the greater heat of vaporization yielded by the TIP4P model in the present simulations. Perhaps this is a consequence of attributing a too big e Lennard-Jones parameter to the oxygen atom. Considering the thermodynamic results and rdf plots, there is no difference between the simulation with 250 and 500 molecules in the simulation box. So one can do the simulation using the smaller number of molecules — one does not loose information and the simulation is faster. The number of configurations used in the present work yield an acceptable convergence of the isothermal compressibility and thermal expansion coefficients. Finally, the disagreement between the experimental value of the mean molar energy and the theoretical one obtained by Degrève and Blum^{18} is in line with the fact that their water model fits very well the old data for the rdf of Soper and Phylips^{16}, which according to more recent data must be overestimated^{21}.

ACKNOWLEDGMENTS

The author is grateful to Dr. Alan K. Soper for the experimental g(O/O(r)) and g(O/H(r)) data and the profitable discussions.

REFERENCES

1. Soper, A. K.; *J. Phys.:Condens. Matter.*, **1997**. *50*, 2717. [ Links ]

2. Gibson, K. D.; Sheraga, H. A.; *J. Comp. Chem.* **1990**,** ***11*, 468. [ Links ]

3. Jorgensen, W. L.; *Chemtracts - Organic Chemistry* **1991**, *4*, 91. [ Links ]

4. Strastsma, T. P.; McCammon, J. A.; *Annu. Rev. Phys. Chem.* **1992**, *43*, 407. [ Links ]

5. Ladanyi, B. M.; Skaf, M. S.; *Annu. Rev. Phys. Chem.* **1993**, *44*, 335. [ Links ]

6. Stillinger, F. H.; Rahman, A.; *J. Chem. Phys.* **1974**, *60*, 1545 [ Links ]

7. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L.; *J. Chem. Phys.* **1983**, *79*, 926. [ Links ]

8. Berendsen, H. J. C.; Grigera J. R.; Straatsma, T. P.; *J. Phys. Chem.* **1987**, *91*, 6269. [ Links ]

9. Brodholt, J.; Sampoli, M.;Vallauri, R.; M*ol. Phys.* **1995**, *86*, 149. [ Links ]

10. Ahlstrom, P.; Wallqvist, A.; Engstrom, S.; Jonsson, B.; *Mol. Phys.* **1989**, *68*, 563. [ Links ]

11. Telleman, O.; Jonsson, B.; *Mol. Phys. ***1987**, *60*, 193. [ Links ]

12. Cordeiro, J. M. M.; unpublished results.

13. Freitas, L. C. G.; Cordeiro, J. M. M.; *J. Mol. Struct*.* (Theochem)* **1995**, *335*, 189. [ Links ]

14. Jorgensen, W. L.; Madura, J. D.; *Mol. Phys.* **1985**, *56*, 1381. [ Links ]

15. Narten, A. H.; Levy, H. A.; *J. Chem. Phys.* **1971**, *55*, 2263. [ Links ]

16. Soper, A. K.; Phylips, M. G.; *Chem. Phys.* **1986**, *107*, 47. [ Links ]

17. Blum, L.; Vericat, F.; Bratko, D.; *J. Chem. Phys.* **1995**, *102*, 1461. [ Links ]

18. Degrève, L.; Blum, L.; *Physica A* **1996**, *224*, 550. [ Links ]

19. Soper, A. K.; *J. Chem. Phys.* **1994**, *101*, 6888. [ Links ]

20. Blum, L.; Degrève, L.; *Mol. Phys.* **1996**, *88*, 585. [ Links ]

21. Soper, A. K.; Bruni, F.; Ricci, M. A.; *J. Chem. Phys.* **1997**, *106*, 247. [ Links ]

22. Soper, A. K.; personal communication.

23. Allen, M. P.; Tildesley, D. J.; *Computer Simulations of Liquids* Oxford University Press, Oxford, 1987. [ Links ]

24. Metropolis, N. A.; Rosembluth, W.; Rosembluth, M. N.; Teller, A. H.; Teller, E.; *J. Chem. Phys.* **1953**, *21*, 108. [ Links ]

25. DIADORIM fortran code by L. C. G. Freitas, Theoretical Chemistry Lab., Chemistry Dept. - UFSCar - S. Carlos - SP - Brazil. [ Links ]

26. Freitas, L. C. G.; *J. Mol. Struct. (Theochem)* **1993**, *282*, 151. [ Links ]

27. Lide, D. R., editor-in-chief; CRC *Handbook of Physics and Chemistry*, 73RD Edition, CRC Press Inc, 1992. [ Links ]

28. Jorgensen, W. L.; personal communication.

29. Jorgensen, W. L.; *Chem. Phys. Letters* **1982**, *92*, 405. [ Links ]

30. Chialvo, A. A., Cummings, P. T.; *J. Chem. Phys.* **1994**, *101*, 4466. [ Links ]

31. Gorbaty, Yu. E. , Kalinishev, A. G.; *J. Phys. Chem*. **1995**, *99*, 5336. [ Links ]

32. note added in proof: Besides the work mentioned above, I have a paper accepted by the J. Mol. Liq. in collaboration with L. C. Gomide Freitas and F. Garbujo (Theoretical studies of liquids by computer simulations: the water-acetone mixture) and another one published in *Int. J. Quant. Chem. ***1997**, *65*, 709. [ Links ]