The density of states and band gap of liquid water by sequential Monte Carlo/Quantum mechanics calculations

The density of states and band gap of liquid water by sequential Monte Carlo/Quantum mechanics calculations

P. Cabral do Couto^I,II; R. C. Guedes^I,II; B. J. Costa Cabral^I,II

^IDepartamento de Química e Bioquímica, Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal

^IIGrupo de Física Matemática da Universidade de Lisboa, Av. Professor Gama Pinto 2, 1649-003 Lisboa, Portugal

ABSTRACT

Electronic properties of liquid water were investigated by sequential Monte Carlo/Quantum mechanics calculations. The density of states (DOS) and HOMO-LUMO gap (E_G) of liquid water have been determined by Hartree-Fock and Density Functional Theory (DFT) calculations. The quantum mechanical calculations were carried out over uncorrelated supermolecular structures generated by the Monte Carlo simulations. The DFT calculations were performed with a modified B3LYP exchange-correlation functional proposed by Abu-Awwad and Politzer which was parametrized to reproduce valence orbital energies in agreement with experimental ionization potentials of the water molecule. We have analyzed the dependence of the DOS and HOMO-LUMO gap on the number of water molecules and on surface effects. Our prediction for E_G is 6:5 ± 0:5 eV in good agreement with a recent experimental prediction of 6.9 eV.

1 Introduction

Liquid water is the most important liquid for life. Reactivity in liquid water is a fundamental process characterizing the behaviour of many living organisms. Therefore, the structure, energetics and electronic properties of liquid water have merited the attention of a large number of experimental and theoretical investigations [1-4]. Some properties of liquid water are fascinating and worth referring: the structure of water is characterized by a complex hydrogen-bond network which leads to a very specific dependence of the density on the thermodynamic state and to the water density anomaly at T=4ºC [5]; the strong polarizability of liquid water, which is related to cooperative effects induced by hydrogen bonding determines its dielectric properties as well as the significant increase of the water molecule dipole moment from 1.85 D in the gas to ~ 2.6 D in the liquid phase [6,7]. The complexity of water makes difficult a first-principle approach to its properties, although some significant contributions have been reported [8-13]. The electronic properties of water, in particular, the density of states (DOS) and band-gap (E_G), have merited the attention of several experimental [14-21] and theoretical [8,21,22] works. In comparison with thermodynamics and structure, electronic properties of liquid water are not very well understood [14,18]. The dominant view is that water can be described as a very large band gap amorphous semiconductor [14].

In the present study we will focus on the electronic properties of liquid water, in particular, on its density of states (DOS) and HOMO-LUMO gap (E_G). To investigate the electronic structure of liquid water we have adopted a sequential Monte Carlo/Quantum mechanical(QM) approach [25,26]. From uncorrelated supermolecular structures generated by the Monte Carlo simulations, QM calculations (HF and DFT) have been carried out to study the electronic properties of liquid water, in particular the density of states and the water band gap. Previous studies of the water electronic properties were based on the Hartree-Fock method [22] and Density Functional Theory [8]. It is known that the HF method is not adequate to describe virtual orbitals. Moreover, Kohn-Sham orbital energies differ from ionization potentials by a factor of about 2 [23]. Recently, Politzer and Abu-Awwad [27,28] proposed a reparametrization of the Becke-3 hybrid exchange-correlation functional that predicts orbital energies in good agreement with experimental ionization potentials. In addition, the modified exchange-correlation functional is able to predict other molecular properties of the water molecule (geometry, dipole moments) in good agreement with experiment [28]. In the present work, this functional will be used to study the electronic properties of liquid water. We have investigated the dependence of the results on the basis set, size of the system, and surface effects. All electrons were included in the calculations. Comparison between our results and data from experiment and other theoretical works indicate that the present approach can provide valuable information on the electronic properties of liquid water.

2 Computational Details

Monte Carlo (MC) simulations of water have been carried out in the isobaric-isothermal (npT) ensemble [29] at T=298 K and p = 1 atm. The SPC potential proposed by Berendsen et al. [30] has been adopted to represent the interaction between the water molecules. The number of water molecules was n = 250 and the interactions have been truncated at a cuttoff radius of 9.0 Å. We have carried out 10⁸ MC steps for equilibration. This run has been followed by 1.25 x 10⁹ additional MC steps. Each step involves the attempt to move one molecule of the system. From the configurations generated by the MC procedure we have selected N=40 and N=50 uncorrelated configurations. These configurations can be defined by determining a correlation time over the Markov chain [25,26].

The density of states of liquid water has been determined by performing a quantum mechanical calculation over the uncorrelated MC configurations [25,26]. By adopting such a procedure, it is possible to get reliable results for many average properties using a relatively small number of configurations. Each configuration corresponds to a cluster including explicitely n_W water molecules on a sphere of radius R. Size dependence has been investigated by carrying out calculations with n_W=35 and n_W=50. To minimize surface effects we have assumed that the water molecules surrounding these clusters can be represented by their SPC charge distribution. Thus, no periodic boundary conditions are being applied. The supermolecular structures are nearly spherical clusters that can be embedded in a rigid (nonpolarizable) charge distribution representing the surrounding water molecules. The number of surrounding molecules is represented by n_c. We have taken n_c = 100.

Quantum mechanical calculations were based on Hartree-Fock (HF) and Density Functional Theory (DFT). To analyze the dependence of the results on the basis set, HF calculations have been carried out with the STO-3G, 3-21G(d), 6-31G(d), 6-31G(d,p), and 6-31+G(d) basis sets [31]. Correlation effects have been investigated with the 6-31+G(d,p) basis set. DFT calculations have been carried out with a modified functional exchange-correlation functional (E_xc) proposed by Abu-Awwad and Politzer [28]. This functional is represented by a Becke's three parameter hybrid [32] given by:

Thus, the exchange energy is written as a combination of HF and local Slater terms plus a Becke's gradient corrected contribution. The correlation energy is represented by a local Vosko-Wilk-Nusair [33] expression combined with a non-local gradient correction. By adopting a Lee-Yang-Parr [34] functional for this contribution, Abu-Awwad and Politzer [28] determined the optimal combination of a, b and g for predicting orbital energies. For the water molecule, using the 6-31+G(d,p) basis set, this combination is a = 0.30, b = 0.72, g = 1.0 [28]. We are using this reparametrization of the B3LYP functional to investigate the electronic properties of liquid water. The QM calculations were carried out with the Gaussian-98 program [35]. The Monte Carlo simulations were carried out with the DICE program [36].

3 The electronic density of states in liquid water

3.1 Basis set dependence

The DOS dependence on the basis set for water clusters with n_w = 35 is illustrated in Fig. 1, where HF calculations with different basis sets are compared. The results are quite sensitive to the quality of the basis set. As expected, in comparison with larger basis sets (6-31G(d), 6-31G(d,p) and 6-31+G(d)), calculations based on the STO-3G and 3-21G basis sets leads to a DOS shifted to higher orbital energies, indicating that these energies are not converged at this level. However, predictions based on the larger basis sets (6-31G(d,p) and 6-31+G(d)) lead to similar distribution of orbital energies. We interpret this behaviour as an indication that these basis sets are adequate to study the DOS in liquid water.

Figure 2 illustrates the dependence of the HOMO and LUMO energies on the basis set. With the exception of the STO-3G calculations, they are very similar. The 6-31G(d) and 6-31G(d,p) results are almost identical and cannot be distinguished in Fig. 2.

Table I reports the averaged HOMO-LUMO gap E_G. Results for E_G with the larger basis sets (6-31G(d), 6-31G(d,p), and 6-31+G(d)) are quite similar (within the error bars), illustrating that although the absolute HOMO and LUMO orbital energies are not truly converged, their differences seem to be less dependent on the basis set.

3.2 Surface effects on the DOS

Our main aim is to investigate the electronic properties of liquid water. However, our uncorrelated supermolecular structures generated by the Monte Carlo simulations are water clusters where surface effects are certainly significant. Therefore, it is necessary to discuss how these effects can be minimized in order to allow a comparison with bulk water. Here, we have adopted the strategy used for predicting the average dipole moment of the water molecule in liquid water [7]. SPC charges, representing n_c surrounding water molecules were added to the quantum system (in the present case, the water cluster with n_w molecules). Surface effects are illustrated in Fig. 3, where the DOS for free clusters with n_w = 35,50 are compared with those for the clusters embedded in the point-charge distribution of n_c = 100 water molecules. The results show that when the charge distribution, representing the surrounding water molecules is included, the distribution of the orbital energies or the DOS is significantly modified. In comparison with free clusters, the DOS of the embedded clusters is characterized by a distribution of the orbital energies in a narrower range. Although we have not included a larger number of polarizing charges, we believe that the main effect is well described. Minor changes should be expected with an increasing number of polarizing charges and it seems reasonable to assume that the DOS of these systems can be used to discuss the electronic properties of the bulk phase. Moreover, as we could expect, the importance of surface effects on the orbital energies is essentially related to the polarization of the water molecules close to the surface of the cluster. To illustrate this effect we have calculated the difference of the electronic density in free and embedded clusters. Isosurfaces representing this quantity are shown in Fig. 4. The main changes of the electronic density concern the water molecules close to the surface. No significant changes are observed in the electronic density of the inner molecules.

Quite interestingly, the DOS for free clusters with n_w=35 and 50 are very similar, indicating that the size of the present systems can adequately describe the orbital energy distribution of water clusters.

Figure 5 illustrates the dependence of the HOMO and LUMO energies on the number of water molecules in the cluster (n_W) and on the inclusion of surrounding point-charges (or surface effects). In keeping with the modifications of the DOS, these energies are not very dependent on the number of water molecules (n_W), but are quite sensitive to the presence of polarizing charges. We note that the presence of embedding charges shifts the LUMO to higher energies and the HOMO to lower energy values, illustrating that E_G is very sensitive to surface effects.

The results for E_G are summarized in Table II. We can observe that at the HF level this quantity changes from 11.7 eV to 14.6 eV when we move from the free cluster including n_w=35 to the cluster embedded in n_c = 100.

3.3 Correlation effects: HF versus DFT calculations

Figure 6 represents the DOS for the water clusters from HF and DFT calculations. In comparison with HF results, DFT calculations shift the occupied orbital energies to higher values, whereas virtual orbital are shifted to lower values. This is very clear for the DOS describing core states where the peak of the orbital energy distribution in the presence of polarizing charges changes from -560 eV to -527 eV. The same effect can be observed in the valence DOS where DFT results are shifted to high energies by ~ 10 eV relative to HF predictions. In addition, in comparison with the HF results, the difference between the top of the valence band and the bottom of the conduction band (that also includes negative LUMO eigenvalues not present in the HF calculations) is significantly reduced. Thus, as it should be expected, introduction of correlation effects, contribute to reduce the band gap. This effect (reduction of E_G by introduction of correlation effects at the DFT level, is illustrated in Fig. 6) is opposite to the effect induced by the addition of embedding charges, which contributes to increase the gap. By adding point charges (or minimizing surface effects), the top of the valence band is shifted to lower energies and the bottom of the conduction band is shifted to higher energies.

Figure 7 shows the behaviour of the HOMO and LUMO energies from DFT/6-31+G(d,p) calculations with and without embedding charges. In keeping with the behaviour of the DOS, introduction of embedding charges increase the HOMO-LUMO gap. In the free clusters the average HOMO energy is -6.1 ±0.5 eV. By adding charges it is shifted to -7.5 ± 0.4 eV. At the same time, the average LUMO energy is -2.4±0.4 eV in the free water clusters and changes to -0.96 ±0.3 eV when charges are added.

The results for the water HOMO-LUMO gap (E_G) are summarized in Table III where we also report data from several experiments and theoretical investigations. Experimental values range from 6.9 eV [21] to 8.9 eV [19]. Our HF/6-31+G(d,p) prediction (13.8 ± 0.4 eV) clearly overestimates the band gap and illustrates the need to introduce correlation effects. It is significantly different from the HF prediction for cubic ice (7.8 eV), which was based on pseudo-potential theory [22]. The low E_G value predicted by Laasonen et al. (4.6 eV) [8] reflects the inadequacy of their exchange correlation functional to predict absolute values of the orbital energies. The present DFT/6-31+G(d,p) result for E_G (6.5 ± 0.4 eV) is close to a recent prediction by Coe et al. (6.9 eV) [21], which was based on the value of the OH^• radical hydration energy. Very recently [24], we have verified that the Coe et al. semi-empirical prediction for this quantity (0.369 eV) was in excellent agreement with the Monte Carlo simulation result for the OH^• radical hydration energy (0.379 eV) [24]. This agreement and the present findings support the Coe et al. prediction and suggest that the energy gap of liquid water can be smaller than it has been indicated by several experimental works.

4 Conclusions

This work reports theoretical results for the density of states (DOS) and the HOMO-LUMO gap of liquid water. We have adopted a sequential Monte Carlo/Quantum Mechanics approach. Quantum mechanics calculations were based on HF and DFT functional theory. DFT calculations were carried out with a modified exchange-correlation functional proposed by Abu-Awwad and Politzer [28]. By using uncorrelated supermolecular strucures, the present approach was able to provide information on the DOS and HOMO-LUMO gap in water clusters. Surface effects were minimized by the introduction of a charge distribution representing the surrounding water molecules and the results indicate that the electronic properties of the embedded clusters can be compared with the electronic properties of liquid water. Our estimation for the HOMO-LUMO gap (E_G) is 6.5 ± 0.4 eV. Although this value is smaller than most of the experimental data, which are in the 7-9 eV range it is close to a recent prediction by Coe et al. (6.9 eV) [21].

Our main conclusion is that the present approach can be of great interest to study the electronic properties of liquids. Concerning its present application to liquid water, we believe that further studies with other functionals able to predict orbital energies in agreement with experiment, and with larger samples could be important for a better understanding of the electronic properties of this fascinating liquid.

Acknowledgments

R. C. Guedes and P. C. do Couto gratefully acknowledge the support of Fundação para a Ciência e a Tecnologia. (PhD grants PRAXIS/XXI/BD/15920/98 and PRAXIS/XXI/BD/6503/2001). This work was partially supported by the Sapiens program of the FCT, Portugal (Grant No. POCTI/43315/QUI/2001). The authors would like to thank Sylvio Canuto and Kaline Coutinho for many interesting discussions and advice.

Received on 29 September, 2003

• [1] Water a comprehensive treatise, Vols 1-4, Ed. by F. Franks, Plenum, New York (1976).
• [2] M. C. R. Symons, Acc. Chem. Res. 14, 179 (1981).
• [3] T. Head-Gordon and G. Hura, Chem. Rev. 102, 2651 (2002).
• [4] J. V. Coe, Int. Rev. Phys. Chem. 20, 33 (2001).
• [5] G. Franseze and H. E. Stanley, Phys. A 314, 508 (2002).
• [6] Y. Tu and A. Laaksonen, Chem. Phys. Lett. 329, 283 (2000).
• [7] K. Coutinho, R. C. Guedes, B. J. Costa Cabral, and S. Canuto, Chem. Phys. Lett. 369, 345 (2003).
• [8] K. Laasonen, M. Sprik, M. Parrinello, and R. Car, J. Chem. Phys. 99, 9080 (1993).
• [9] M. Sprik, J. Hutter, and M. Parrinello, J. Chem. Phys. 105, 1142 (1996).
• [10] P. L. Silvestrelli and M. Parrinello, Phys. Rev. Lett. 82, 5415 (1999).
• [11] P. L. Silvestrelli and M. Parrinello, J. Chem. Phys. 111, 3572 (1999).
• [12] M. Boero, K. Terakura, T. Ikeshoji, C. C. Liew, and M. Parrinello, Phys. Rev. Lett. 85, 3245 (2000).
• [13] M. Boero, K. Terakura, T. Ikeshoji, C. C. Liew, and M. Parrinello, Journal of Chemical Physics, 115, 2219 (2001).
• [14] F. Williams, S. P. Varma, and S. Hillenius, J. Chem. Phys. 64, 1549 (1976).
• [15] B. Baron, D. Hoover, and F. Williams, J. Chem. Phys. 68, 1997 (1978).
• [16] D. Grand, A. Bernas, and E. Amouyal, Chem. Phys. 44, 73 (1979).
• [17] T. Goulet and J.-P. Jay-Gerin, J. Phys. Chem. 93, 7532 (1989).
• [18] T. Goulet, A. Bernas, C. Ferradini, and J.-P. Jay-Gerin, Chem. Phys. Lett. 170, 492 (1990).
• [19] A. Bernas, C. Ferradini, and J.-P. Jay-Gerin, Chem. Phys. 222, 151 (1997).
• [20] A. Bernas, C. Ferradini, and J.-P. Jay-Gerin, J. Photochem. and Photobio. A: Chem. 117, 171 (1998).
• [21] J. V. Coe, A. D. Earhart, M. C. Cohen, G. J. Hoffman, H. W. Sarkas, and K. H. Bowen, J. Chem. Phys. 107, 6023 (1997).
• [22] G. P. Parravicini and L. Resca, Phys. Rev. B 8, 3009 (1973).
• [23] X. Hua, X. Chen, and W. A. Goddard III, Phys. Rev. B 55, 16103 (1997).
• [24] P. Cabral do Couto, R. C. Guedes, B. J. Costa Cabral, and J. A. Martinho Simőes, J. Chem. Phys. 119, 7344 (2003).
• [25] S. Canuto and K. Coutinho, Adv. Chem. Phys. 28, 90 (1997).
• [26] K. Coutinho and S. Canuto, J. Chem. Phys. 113, 9132 (2000).
• [27] P. Politzer and F. Abu-Awwad, Mol. Phys. 95, 681 (1998).
• [28] F. Abu-Awwad and P. Politzer, J. Comput. Chem. 21, 227 (2000).
• [29] D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic, San Diego, 1996).
• [30] H. J. C Berendsen, J. P. M. Postma, W. F. van Gusteren, and J. Hermans in Intermolecular Forces, edited by B. Pullman (Reidel: Dordrecht, The Netherlands, 1981), p 331.
• [31] W. J. Hehre, L. Radom, P. v.R. Schleyer, and J. A. Pople, Ab initio Molecular Orbital Theory (John Wiley&Sons, New York, 1986).
• [32] A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
• [33] S. H. Vosko, L. Wilk, and M. Nusair, Canadian J. Phys. 58, 1200 (1980).
• [34] C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
• [35] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, V. G. Zakrzewski, J. A. Montgomery, R. E. Stratman, J. C. Burant, S. Dapprich, J. M. Millan, A. D. Daniels, K. N. Nudin, M. C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G. A. Peterson, P. Y. Ayala, Q. Cui, K. Morokuma, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. Cioslowski, J. V. Ortiz, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe, P. M. W. Gill, B. G. Johnson, W. Chen, M. W. Wong, J. L. Andres, M. Head-Gordon, E. S. Repogle, and J. A. Pople, GAUSSIAN-98, Gaussian Inc., Pittsburgh, PA, 1998.
• [36] K. Coutinho and S. Canuto DICE: A general Monte Carlo program for liquid simulation (University of Săo Paulo, Brazil, 2000).

Publication Dates

• Publication in this collection
  27 Apr 2004
• Date of issue
  Mar 2004

History

• Received
  29 Sept 2003