SciELO - Scientific Electronic Library Online

 
vol.21 issue6Reactions of Cp'(CO)2MnPPh2H with CH3COCl and CH3S(O)2Cl in THF/triethylamine: evidences of the first complex stabilisation of a phosphorus-homologous of the sulfonamidesTriterpenic esters from himatanthus sucuuba (Spruce) Woodson author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

Share


Química Nova

Print version ISSN 0100-4042On-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 25oC 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 systems1. 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 task2-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 models6-8 and others including polarization effects9-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 models12. 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 solutions13,32. For water I used the TIP4P model developed by Jorgensen et al7,14. This model is, together with the SPC and SPC/E models8, 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 Levy15. 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 al7. 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 experiments7. Experimental rdfs for water at 25oC and 1 atm were more recently obtained via neutron diffraction experiments by Soper and Philips16, 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 data7. Since then, some analytical models for water have been proposed by Blum et al17,18 that represent very well the rdf published by Soper and Philips16. 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 reliably19. Blum and Degreve do not refer to this work in a recently published paper20. 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 400oC21. 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 accurately22. The g(O/O(r)) plot reported agrees very well with the earlier X-ray results reported by Jorgensen and coworkers7. 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 al7,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 al7,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 Eab between molecules a and b is represented by a sum of Coulomb and Lennard-Jones potentials centered on the sites, that is

(1)

where rij is the distance between site i in a and site j in b and qi and qj are point charges located respectively on the i and j molecular sites. For each site k, the parameters Akk and Bkk are given by Akk=4eksk12 and Bkk = 4eksk6 , where ek and sk are the Lennard-Jones parameters for the k-th site.

Monte Carlo Simulations with Metropolis' importance sampling24 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. Freitas25. 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 energy14. 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 Madura14, 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 way23. 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 2x105 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 Madura14 and with experimental data. The origin of the experimental intermolecular energy Ei 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 Madura14. 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 procedure26 have shown values in accordance with the results of Jorgensen and Madura14. 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 procedure28. The experimental intermolecular energy was computed with the equation DHvap = -Ei + P[V(g) - V(l)] = -Ei + 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 Madura14, there is little error in making this approximation. The heat capacity Cp, the isothermal compressibility k and the thermal expansion a have been computed from the standard fluctuation formulas23. The rates of convergence for these properties are slower than for the energy or density, especially k and a29. The computed configurational contribution to the Cp 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 et al7. 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 Madura14 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 Philips16 is presented for comparison. As one can observe, there is a good agreement between the g(O/O(r)) obtained by Soper et al21 and that obtained by Jorgensen et al7 using the Narten and Levy experimental data14. The behavior verified for the TIP4P model in the present work is the same as that verified earlier by Jorgensen et al7,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 al7.

 

 

 

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 al7 and is a little greater than the experimental one21. 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 account30-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 Blum18 is in line with the fact that their water model fits very well the old data for the rdf of Soper and Phylips16, which according to more recent data must be overestimated21.

 

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.; Mol. 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 ]

Creative Commons License All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License