density functional theory

Chemical Physics Letters 460 (2008) 466–469

Contents lists available at ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

Electronic properties of liquid water by sequential molecular dynamics/density functional theory Claude Millot a, Benedito J. Costa Cabral b,* a

Équipe Chimie & Biochimie Théorique, UMR CNRS-UHP 7565, Nancy University, University Henri Poincaré, F-54506 Vandoeuvre-lès-Nancy Cedex, France Departamento de Química e Bioquímica, Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal and Grupo de Física Matemática da Universidade de Lisboa, Av. Professor Gama Pinto 2, 1649-003 Lisboa, Portugal

b

a r t i c l e

i n f o

Article history: Received 20 May 2008 In final form 17 June 2008 Available online 24 June 2008

a b s t r a c t Electronic properties of liquid water were analysed by a sequential molecular dynamics (MD)/density functional theory approach. MD simulations are based on a polarisable model for water. Emphasis was placed on the prediction of the water dipole moment, liquid state polarisability, ionisation potential (IP), and vertical electron affinity. The dipole moment of the water molecule in liquid water is not dependent on the number of molecules included in the quantum mechanical calculations. The polarisability of the water molecule in liquid water is 4% lower than its gas phase value. The IP of liquid water (9.7 ± 0.06 eV) is in good agreement with recent experimental data. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction Water is characterized by a complex hydrogen bonding (HB) network with unique properties, which are strongly associated with cooperative effects [1]. These effects are responsible by anomalous properties including the dependence of the density on the thermodynamic state and the water density anomaly at T = 4 °C. They also lead to a significant increase of the water molecule dipole moment from 1.85 D in the gas to 2.6 D in the liquid phase [2]. Although electronic properties of water [3–9] are of crucial importance for a better knowledge of chemical reactivity in solution, they are not very well understood [7]. Two fundamental processes strongly associated with the electronic structure of the HB network are ionisation [7] and attachment of an excess electron to water [10,11]. It is known that molecular clusters can bind an excess electron when their total multipole moments are above a critical value [12]. This feature leads, in some cases, to magic numbers that can be related with the total dipole moment of the aggregate [11]. Several works discussed the structure and electronic properties of water anionic clusters [11]. Most of these studies were carried out for low energy aggregates in their most stable configurations, although an investigation on excess electron localization sites in neutral water clusters at 200 and 300 K has been recently reported [13]. For liquid water, in contrast with low energy configurations, it should be expected that thermal effects will induce fluctuations in the structure of the HB network, which can be correlated with the local charge distribution. In this Letter, we report a sequential molecular dynamics (MD)/density func-

* Corresponding author. E-mail address: [email protected] (B.J. Costa Cabral). 0009-2614/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2008.06.059

tional theory (DFT) approach to the electronic properties of liquid water, with emphasis on the evaluation of the dipole moment, electronic polarisability, ionisation potential, and vertical electron affinity. The Letter is organized as follows. Initially, we present the technical details of the theoretical approach. DFT results for the average dipole moment of the water monomer in water and the total average dipole moment of water aggregates with a different number of water molecules (nw) embedded in the charge distribution of polarisable water are then reported. Finally, results for the ionisation potential and vertical electron affinity of liquid water are presented and compared with experimental data.

2. Computational details The MD simulation has been carried out in the (NVE) microcanonical ensemble with the Niesar–Corongiu–Clementi (NCC) model [14] for polarisable water with N = 343 molecules in a cubic box with periodic boundary conditions at a density q = 0.997 g cm3 and average temperature T = 298 K. With the NCC model the interaction energy is calculated as a sum of two-body contributions (short range and electrostatic interactions) and a many-body polarisation term. The electrostatic interactions are modelled through three partial charges per molecules located on the hydrogen atoms and an extra site on the C2v axis and the polarisability is represented through the introduction of two local bond anisotropic dipolar polarizabilities. Long-range electrostatic interactions for both permanent charges and induced dipoles were taken into account by the reaction field method [15] with a molecular-based spherical cutoff of 10.8 Å and the dielectric constant of the surrounding medium equal to infinity. At each time step, the induced moments were obtained by an iterative process initiated with the

467

C. Millot, B.J. Costa Cabral / Chemical Physics Letters 460 (2008) 466–469

induced moments of the previous time step. On average, 5.4 iterations were necessary to get all induced dipole components converged with an accuracy of 0.001 D, and to avoid the convergency test 7 iterations have been done at each time steps. The molecules are considered as rigid bodies with an OH bond length of 0.9572 Å and an HOH angle of 104.52°. The integration of the equations of motion was carried out with a quaternion based algorithm proposed by Fincham [16] with a timestep of 2 fs. To avoid energy drift due to the cutoff, velocities have been rescaled every 100 time steps. After equilibration, 106 steps (2 ns) were generated. One hundred configurations (saved every 20 ps) were then selected for the calculation of the electronic properties. Successive configurations generated by MD are strongly correlated and when the property of interest involves a high computational effort, the use of uncorrelated structures is of crucial importance for evaluating averages over a relatively small number of representative configurations [17]. The statistical error of average properpffiffiffiffiffi ties will be reported as r= nd , where r is the standard deviation and nd is the number of data. For each selected configuration, water aggregates of different sizes nw = 2, 3, 5, 6, 8, 10, 12, 13, 15, 20 were defined. They include a central molecule (the first one) and the closest nw1 water molecules. For the DFT calculations, the aggregates were embedded in the charge distribution of the closest one hundred surrounding water molecules. The charge distribution of the embedding molecules is represented by the permanent charges of the NCC model as well as by the charges representing the induced dipole moments on each water monomer. For the quantum mechanical calculations, the self-consistent induced dipole moments associated with the local O–H bond polarisability on each water molecule were reproduced by placing two point charges of (+) 0.2 a.u., separated by a variable vector centred at a distance of 0.43296 Å from the oxygen atom and along the O–H bond. The introduction of embedding charges mimics the presence of the surrounding water molecules and no periodic boundary conditions were applied for evaluating electronic properties in the liquid phase. The importance of embedding was recently discussed [18] and only results for embedded aggregates will be presently reported. For the different aggregates, dipole moments were evaluated by using charges fitted to the electrostatic potential (ESP charges). ESP charges were calculated by the Breneman method [19]. For a cluster with nw water molecules the first ionisation potential (IP) was calculated as the energy difference between the neutral and ionized aggregates. The vertical electron affinity (VEA) was calculated as the energy difference between the neutral and the anionic aggregates in the geometry of the neutral system. Therefore, no geometry relaxation is involved in the DE calculations of the IPs and VEAs. DFT calculations were carried out with the modified Perdew-Wang functional (MPW1PW91) proposed by Adamo and Barone [20], and the Dunning’s aug-cc-pVDZ and d-aug-cc-pVxZ basis sets (x = D,T) [21]. DFT calculations were performed with the GAUSSIAN-03 suite of programs [22].

3. Results and discussion 3.1. Liquid state dipole moment The MPW1PW91/aug-cc-pVDZ dipole moment of the isolated water molecule in the NCC geometry is 1.87 D, in very good agreement with the experimental value (1.85 D). The average dipole moment of the water monomer (l) and the average total dipole moment for each aggregate with nw water molecules (lT,nw) are reported in Table 1. A comprehensive discussion on the evaluation of the average monomer dipole moment l in water and its relationship with the dielectric constant e0 was reported by Soetens

Table 1 Average dipole moment l of the water monomer in liquid water, total average dipole (lT,nw) for aggregates with nw water molecules from sequential molecular dynamics/ DFT calculations nw

l

lT,nw

1 2 3 5 6 8 10 12 13 15 20

2.49 ± 0.01 2.54 ± 0.02 2.56 ± 0.02 2.56 ± 0.01 2.56 ± 0.01 2.54 ± 0.01 2.56 ± 0.02 2.56 ± 0.02 2.57 ± 0.02 2.56 ± 0.02 2.54 ± 0.01

2.46 ± 0.01 4.34 ± 0.06 5.59 ± 0.13 7.09 ± 0.21 7.57 ± 0.25 8.45 ± 0.31 9.26 ± 0.36 9.96 ± 0.39 10.50 ± 0.42 11.18 ± 0.48 13.00 ± 0.56

Other values

2.75 ± 0.02a; 2.47b; 2.59 ± 0.14c; 2.66–2.69, 2.92d; 2.95–3.00e; 2.64 ± 0.17f; 2.74–2.48g

DFT calculations were carried out at the MPW1PW91/aug-cc-pVDZ level. Values in D. a Average dipole moment from the NCC polarisable model. b From Delle Site et al. [25]. c From Coutinho et al. [26]. d From Chalmet and Ruiz-López [27,28]. e From Silvestrelli and Parrinello [29]. f From Tu and Laaksonen [30]. g From Kongsted et al. [31].

et al. [23]. In that study, examination of different polarisable potentials of the literature indicates that to reproduce correctly e0, a model would lead to l in the range 2.5–2.6 D, supporting an idea proposed by Sprik [24]. Here, two aspects should be stressed. Firstly, the average dipole moment of the water monomer in liquid water is not significantly dependent on the number of water molecules explicitly included in the quantum mechanical calculations. For example, l is 2.54 ± 0.02 D (nw = 2), 2.56 ± 0.02 D (nw = 15), and 2.54 ± 0.01 D (nw = 20). In addition, the average values are quite close to other theoretical estimates of the average water dipole moment in liquid water [25–28,30]. In contrast with previous estimates [26], which were based on effective charge distributions in liquid water, the electrostatic field of the surrounding water molecules now includes explicitly polarization effects. Hybrid Quantum Mechanics-CCSD/Molecular Mechanics studies for one QM water molecule embedded in a cluster of 127 MM water molecules, leaded to dipole moments of 2.74 D and 2.48 D when the polarisability of the surrounding MM molecules was or was not taken into account, respectively [31]. The average dipole moment based on the polarisable NCC model (2.75 D) is 0.2 D above the quantum mechanical prediction. The average total dipole moment of the aggregates lT,nw exhibits significant fluctuations and increases from 2.46 ± 0.01 D (nw = 1) to 13.00 ± 0.48 D (nw = 20). 3.2. Liquid state polarisability The electronic polarisability of the water molecule in liquid water has been the subject of previous studies, which indicated that the average polarisability of the water molecule is 79% lower relative to its gas phase value [32]. This issue is relevant for a correct parametrisation of polarisable force-field models of water. In the NCC model the isotropic polarisability of the water molecule [a ¼ 13 ðaxx þ ayy þ azz Þ] is 9.29 a.u., which is lower than the experimental gas phase value (9.78 a.u.). The polarisability of the isolated water molecule is 9.67 a.u. at the MPW1PW91/d-aug-ccpVDZ level. If we assume that the liquid state polarisability can be estimated by using only one QM water molecule, as in the case of the dipole moment, our prediction for the average polarisability of water is 9.31 ± 0.02 a.u., which is 4% lower than the gas phase value. The above prediction for the liquid phase polarisability

C. Millot, B.J. Costa Cabral / Chemical Physics Letters 460 (2008) 466–469

VEA /eV

0.5 0 -0.5 -1 -1.5

IP /eV

(9.31 a.u.) is very close to the value of the NCC model (9.29 a.u.) [14]. Calculations with the smaller basis set (aug-cc-pVDZ) lead to gas and liquid phase polarisabilities of 9.11 and 8.77 ± 0.01 a.u., respectively. However, the liquid to gas phase decrease of a (4%) is similar to the one observed with the d-aug-cc-pVDZ basis set. The polarisability of water was calculated by carrying out a sum-over-states (SOS) and its gas phase value was also confirmed by a finite field calculation. The calculations with the d-aug-ccpVDZ basis set included 265 states (Ns) and the final reported value (9.67 a.u.) for the gas phase polarisability was converged within 103 a.u. when Ns = 200. A similar convergence and accuracy was observed for the calculations with the central water molecule embedded in the charge distribution of liquid water. We have verified that by using the larger d-aug-cc-pVTZ basis set and 605 states in the SOS, the gas phase polarisability of the water molecule is 9.71 a.u., and the average value of the liquid phase is 9.36 a.u., which also means a 4% decrease relative to the gas phase. Therefore, our calculations indicate that the decrease of the water isotropic polarizability in the liquid relative to its gas phase value is not dependent on the basis set or SOS procedure.

12 11 10 2.6

μ /D

468

2.5 2.4

2

4

6

8

10

12

14

16

18

20

22

nw Fig. 1. Dependence of the dipole moment (l), ionisation potential (IP), and vertical electronic affinity (VEA) with the number of water molecules in the quantum system nw.

3.3. Ionisation potential and vertical electron affinity From a DE calculation, the MPW1PW91/aug-cc-pVDZ gas phase ionisation potential of water is 12.65 eV, which is in excellent agreement with experiment (12.62 eV) [33]. Average values for the IP and VEA from sequential MD/DFT calculations are reported in Table 2. The IP is dependent on the number of water molecules included in the DFT calculations. It decreases from 11.66 ± 0.07 eV (nw = 2) to 9.71 ± 0.06 eV (nw = 20), which is our best estimate. Experimental values for the first ionisation potential of liquid water are also reported in Table 2. They are in the 9.3 [34]–10.06 [35] eV range. Our best estimate for IP is in good agreement with experimental results [7,34] and it is 0.2 eV below the most recent experimental data [7]. It should be observed that the present estimate for the water first ionisation potential is based on DE calculations. This means, in contrast with predictions based on electron binding energies [8], that electronic density relaxation in the ionized species is taken into account and may lead to a stabilizing effect. Average values for vertical electron affinities are reported in Table 2. They exhibit a strong dependence on nw that is illustrated Table 2 Average values for the first ionisation potential (IP) and vertical electron affinity (VEA) from sequential molecular dynamics/DFT calculations nw

IP

VEA

1 2 3 5 6 8 10 12 13 15 20

12.82 ± 0.06 11.66 ± 0.06 11.10 ± 0.04 10.63 ± 0.04 10.51 ± 0.04 10.37 ± 0.04 10.24 ± 0.06 9.97 ± 0.03 10.02 ± 0.06 9.87 ± 0.06 9.71 ± 0.06

1.41 ± 0.04 1.05 ± 0.05 0.70 ± 0.05 0.44 ± 0.04 0.34 ± 0.04 0.13 ± 0.04 0.02 ± 0.04 0.07 ± 0.04 0.11 ± 0.04 0.20 ± 0.04 0.44 ± 0.05

Other values

9.3 ± 0.a;9.9b; 10.06c 10.75 ± 0.13e

1.2 ± 0.1d; 0.79 ± 0.08e; 0.12 ± 0.05f

DFT calculations were carried out at the MPW1PW91/aug-cc-pVDZ level. Values in eV. Results in italic are our best estimates. a Experimental value from Watanabe et al. [34]. b Experimental value from Winter and Faubel [7]. c Experimental value from Delahay [35]. d Experimental value for V0 from Grand et al. [3]. e Theoretical result from do Couto et al. [8]. f V0 from Coe et al. [6].

in Fig. 1, where the behaviour of VEA, IP, and l with the number of water molecules in the quantum system are compared. VEAs are negative for nw below 12. This means that in average, these aggregates (in the polarizing field of the surrounding water molecules) cannot stabilize, vertically, an excess electron. Our best estimate of VEA is 0.44 ± 0.05 eV for nw = 20. In our calculations it is assumed that the electron attachment is rapid with respect to the nuclear coordinates and that the anion is formed in the geometry of the neutral aggregate. VEA can be compared with minus the conduction band edge of water (V0) [6]. Values of V0 are also reported in Table 2. The present estimate of V0 is in good agreement with a recent prediction (0.79 ± 0.08 eV), which was based on Green’s function evaluation of electron binding energies, and supports the view [6] that V0 is significantly smaller than the typical literature value (1.2 ± 0.1 eV) [3]. 4. Conclusions Results from sequential MD/DFT calculations for the dipole moment, electronic polarisation, ionisation potential and vertical electron affinity of liquid water are reported. The dynamics was generated using the NCC polarisable model of liquid water. The calculation of the electronic properties was carried by using one hundred uncorrelated configurations from the MD simulation. Clusters of different sizes were embedded in the charge distribution of the NCC model. It has been verified that the water dipole moment is not dependent on the number of water molecules explicitly included in the quantum mechanical calculations. We are providing further evidence on the decrease of the isotropic polarisability of the water molecule in liquid water relative to its gas phase value. We predict that in comparison with the gas phase, a is 4% lower in liquid water. The ionisation potential and the vertical electron affinity are size dependent and cannot be estimated as a ‘local property’. Our best estimate for the water ionisation potential is in good agreement with experimental results. The vertical electron affinity was compared with minus the conduction band edge of water (V0) and it was concluded, that in agreement with recent works [6,9], it is smaller than the typical literature value [3]. Acknowledgment This work was partially supported by Fundação para a Ciência e a Tecnologia (FCT), Portugal (Grants No. POCI/MAT/55977/2004 and PTDC/QUI/68226/2006).

C. Millot, B.J. Costa Cabral / Chemical Physics Letters 460 (2008) 466–469

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]

S.S. Xantheas, Chem. Phys. Lett. 258 (2000) 225. C.A. Coulson, D. Eisenberg, Proc. R. Soc. London, A 291 (1966) 445. D. Grand, A. Bernas, E. Amouyal, Chem. Phys. 44 (1979) 73. P. Han, D.M. Bartels, J. Phys. Chem. 94 (1990) 5824. A. Bernas, C. Ferradini, J.-P. Jay-Gerin, Chem. Phys. 222 (1997) 151. J.V. Coe, A.D. Earhart, M.H. Cohen, G.J. Hoffman, H.W. Sarkas, K.H. Bowen, J. Chem. Phys. 107 (1997) 6023. B. Winter, M. Faubel, Chem. Rev. 106 (2006) 1176. P. Cabral do Couto, B.J.C. Cabral, S. Canuto, Chem. Phys. Lett. 429 (2006) 129. P. Cabral do Couto, B.J.C. Cabral, J. Chem. Phys. 126 (2007) 014509. J.V. Coe, et al., J. Chem. Phys. 92 (1990) 3980. C.G. Bailey, J. Kim, M.A. Johnson, J. Phys. Chem. 100 (1996) 16782. H. Abdoul-Carime, J.P. Schermann, C. Desfrançois, Few-Body Syst. 31 (2002) 183. L. Turi, A. Madaràsz, P.J. Rossky, J. Chem. Phys. 125 (2006) 014308. U. Niesar, G. Corongiu, E. Clementi, G.R. Kneller, D.K. Bhatacharya, J. Phys. Chem. 94 (1990) 7949. M. Neumann, J. Chem. Phys. 82 (1985) 5663. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. K. Coutinho, S. Canuto, Adv. Quantum Chem. 28 (1997) 89.

469

[18] T.S. Almeida, K. Coutinho, B.J.C. Cabral, S. Canuto, J. Chem. Phys. 128 (2008) 014506. [19] C.M. Breneman, K. Wiberg, J. Comput. Chem. 11 (1990) 361. [20] C. Adamo, V. Barone, J. Chem. Phys. 108 (1998) 664. [21] D.E. Woon, T.H. Dunning Jr., J. Chem. Phys. 98 (1993) 1358. [22] M.J. Frisch et al., GAUSSIAN-03, Rev. C.02 Gaussian Inc., Pittsburgh, PA, 2003. [23] J.-C. Soetens, M.T.C. Martins-Costa, C. Millot, Mol. Phys. 94 (1998) 577. [24] M. Sprik, J. Chem. Phys. 95 (1991) 6762. [25] L. Delle Site, A. Alavi, R.M. Lynden-Bell, Mol. Phys. 96 (1999) 1683. [26] K. Coutinho, R.C. Guedes, B.J.C. Cabral, S. Canuto, Chem. Phys. Lett. 369 (2003) 345. [27] S. Chalmet, M.F. Ruiz-López, J. Chem. Phys. 115 (2001) 5220. [28] S. Chalmet, D. Rinaldi, M.F. Ruiz-López, Int. J. Quantum Chem. 84 (2001) 1559. [29] P.L. Silvestrelli, M. Parrinello, J. Chem. Phys. 111 (1999) 3572. [30] Y. Tu, A. Laaksonen, Chem. Phys. Lett. 329 (2000) 329. [31] J. Kongsted, A. Osted, K.V. Mikkelsen, O. Christiansen, Chem. Phys. Lett. 364 (2002) 379. [32] A. Morita, J. Comput. Chem. 23 (2002) 1466. [33] O. Dutuit, A. Tabche-Fouhaile, I. Nenner, H. Frolich, P.M. Guyon, J. Chem. Phys. 83 (1985) 584. [34] T. Watanabe, H. Gerischer, J. Electroanal. Chem. Interfacial Electrochem. 122 (1981) 73. [35] P. Delahay, Acc. Chem. Res. 15 (1982) 40.