Water adsorption on stepped ZnO surfaces from MD simulation

Surface Science 604 (2010) 741–752

Contents lists available at ScienceDirect

Surface Science journal homepage: www.elsevier.com/locate/susc

Water adsorption on stepped ZnO surfaces from MD simulation David Raymand a, Adri C.T. van Duin b, Daniel Spångberg a, William A. Goddard III c, Kersti Hermansson a,* a

Department of Materials Chemistry, The Ångström Laboratory, Box 538, 751 21 Uppsala, Sweden Department of Mechanical and Nuclear Engineering, Penn State University, University Park, PA 16802, USA c Materials and Process Simulation Center (MSC), California Institute of Technology, Pasadena, CA 91125, USA b

a r t i c l e

i n f o

Article history: Received 15 September 2009 Accepted for publication 11 December 2009 Available online 29 December 2009 Keywords: Zinc oxide Water Solid–gas interfaces Construction and use of effective interatomic interactions Molecular dynamics Density-functional calculations

a b s t r a c t This work presents a ReaxFF reactive force-field for use in molecular dynamics simulations of the ZnO– water system. The force-field parameters were fitted to a data-set of energies, geometries and charges derived from quantum-mechanical B3LYP calculations. The presented ReaxFF model provides a good fit to the QM reference data for the ZnO–water system that was present in the data-set. The force-field has been used to study how water is adsorbed, molecularly or dissociatively, at monolayer coverage on flat and stepped ZnO surfaces, at three different temperatures (10 K, 300 K, and 600 K). The stepped  0Þ-surface. Equilibsurfaces were created by introducing steps along the (0 0 0 1)-direction on the ð1 0 1  0Þ terraces, resulting in a half rium between molecular and dissociated water was observed on the ð1 0 1 dissociated, half molecular water monolayer. The equilibrium between dissociated and molecular water on the surface was found to be reached quickly (<10 ps). When water molecules desorb and the coverage  0Þ terraces, while steps remain largely hydroxfalls, the 1:1 water–hydroxyl ratio is maintained on ð1 0 1 ylated. The results show that structures that promote hydrogen bonding are favored and that the presence of steps promotes an increased level of hydroxylation in the water monolayers. Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction Zinc oxide (ZnO) is a material with a partly covalent and partly ionic character, which gives it many interesting properties. ZnO or ZnO-based compounds are used in heterogeneous catalysis to produce many chemicals, e.g. in the formation of methanol from CO and H2 [1]. Zinc oxide in its normal crystalline form has a hexagonal wurtzite-type structure, with a surface morphology dominated  0Þ; ð1 1 2  0Þ, (0 0 0 1) and ð0 0 0 1Þ  [2]. For ZnO by four faces: ð1 0 1  0Þ and ð1 1 2  0Þ faces together account for about powders, the ð1 0 1 80% of the total surface area (as determined from FTIR and HRTEM investigations [3]), and they therefore essentially determine the adsorption properties [3]. Water plays a key role in many important processes taking place on the ZnO surface and it is therefore important to understand how water interacts with the ZnO surface. For example, the reverse water–gas shift reaction, i.e. H2 + CO2 ? H2O + CO always has to be considered in the synthesis of methanol [4]. Most metal oxides chemisorb water molecules from the atmosphere when possible. The topic of the interaction between the  0Þ-surface and water has been investigated by Meyer ZnOð1 0 1 and co-workers in a set of publications [5–7]. Using both experiment (STM) and quantum-mechanical calculations (DFT), they pro* Corresponding author. Tel.: +46 18 471 3725; fax: +46 18 513 548. E-mail address: [email protected] (K. Hermansson). 0039-6028/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.susc.2009.12.012

posed that at monolayer coverage a water layer on the non-polar  0Þ-surface consists of alternating dissociated and molecZnOð1 0 1 ular water. The formed surface hydroxyl groups and adsorbed water molecules are ordered in a (2  1)-pattern. In their most recent, purely theoretical, publication on this topic [7], they concluded that the interaction between the water molecules themselves activates the dissociation reaction: isolated water molecules remain intact until a second water molecule donates a hydrogen bond, activating the first molecule for dissociation. This cooperative hydrogen-bonded network weakens the O–H bonds already activated by the Zn2+ ions on the surface and facilitates dissociation by further lowering the reaction barrier [7].  0Þ-surface, the MgO(1 0 0)-surface is known to Like the ZnOð1 0 1 be able to dissociate water molecules without promotion by a surface defect such as steps, kinks or vacancies. Giordano et al. [8] found, using DFT calculations, that in a water monolayer, some of the adsorbed water molecules dissociate due to the interactions within the ad-layer. In a subsequent DFT study, Odelius [9] found that mixed molecular and dissociative adsorption of water was the most favorable on MgO(1 0 0), in agreement with Ref. [8]. According to Odelius, the energy gain from the hydrogen bonding in a molecularly adsorbed monolayer is similar to the energy loss due to the accompanying non-optimal adsorption structure. On the other hand, in the case of mixed molecular and dissociative adsorption, the dissociated water is stabilized by hydrogen bond donation from an intact water molecule, and moreover, the hydroxyl group accepts a

742

D. Raymand et al. / Surface Science 604 (2010) 741–752

hydrogen bond from the surface hydroxyl group formed from H and an oxygen ion from the MgO-lattice [9]. It is reasonable to assume  0Þ, where the surface that the same conditions apply to ZnOð1 0 1 structure offers similar hydrogen-bonding conditions.1 Subsequent to the studies of Meyer et al. [5], Cooke et al. [10]  0Þ and ZnOð1 1 2  0Þ, studied a monolayer of water on ZnOð1 0 1 using DFT calculations and force-field calculations. They studied ZnO surfaces under different external conditions by considering  0Þ-surface, they the effect of non-stoichiometry. For the ZnOð1 1 2 found that fully dissociative adsorption is favored. The reason reported is that this adsorption structure offers a stabilizing hydrogen-bond network, while molecular adsorption offers little or no hydrogen bonding. The force-field presented in the current study was applied in molecular dynamics simulations in Ref. [11], which reproduced  0Þ presented by Meyer et al. the adsorption pattern on ZnOð1 0 1 and predicted that the pattern would remain even if the temperature increases. Apart from these studies, the theoretical literature concerning the atomic-level structure of water/ZnO interfaces is quite scarce. For a more general review of theoretical studies of ZnO we refer to Ref. [12]. The experimental literature of the water/ZnO system is more abundant and we refer to Ref. [4] for a comprehensive review. However, it is difficult to identify whether water is adsorbed molecularly or dissociatively using experimental techniques. The difficulties lie in the similarity between the vibrational spectra of H2O and OH bound to the surface. Determining whether the adsorption is molecular or dissociative is important, since the dissociation products (OH, H and O) are chemically very different from water [13]. The type of adsorption may depend on the level of coverage, which adds to the complexity. For the Cr2O3(0 0 0 1)surface, for example, it has been found, using DFT calculations, that the adsorption was molecular at low water coverage, mixed molecular and dissociative at medium coverage, and fully dissociative at high coverage. However, at full coverage the adsorption was found to revert to a mixed molecular and dissociative type [14]. In this work the reactive force-field model (ReaxFF) [15,16] has been used to model interactions between ZnO surfaces and water molecules. The force-field approach allows us to perform molecular dynamics (MD) simulations for models consisting of thousands of individual atoms during hundreds of picoseconds or more, rather than up to a hundred atoms and tens of picoseconds which is typical for a so-called ab inito-MD simulation. The ReaxFF forcefield is able to simulate the breaking and formation of bonds during dynamics. The force-field parameters within the ReaxFF framework are developed using results of quantum-mechanical (QM) calculations. The ReaxFF framework was initially developed for hydrocarbons [15,16], and has since then been applied to many such systems. It has also been employed in the study of several oxide systems, MoOx, VxOy and BixOy [17,18] oxides and Si/SiO2 [16] and Al/Al2O3 [19] oxide interfaces. We have previously published parameters for ZnO [20]. In the current study we extend the ZnO force-field to include the interaction with water (see Section 2.2). The layout of the paper is as follows. Section 2 describes the methods used, including details about the ReaxFF model the force-field optimization procedure and the quantum-mechanical calculations used to generate the data-set against which the force-field is trained. In Section 3, we assess the quality of the new force-field model and apply it in molecular dynamics simulations to study the effect of temperature and surface steps on the 1 In the discussion in the present paper no distinction will be made between favorable H  O interactions and hydrogen bonding and we will therefore not discuss whether (or not) it is appropriate to label H  O interactions between two negative OH species as H-bonds.

reactivity of the ZnO surface towards water dissociation, especially with respect to the character of the adsorption, molecular or  0Þ-terraces. dissociative, around steps on surfaces with ZnOð1 0 1 Geometry optimizations were also performed to find the dependence of the water adsorption energy on the presence of a surface step.

2. Method 2.1. The ReaxFF reactive force-field method In the ReaxFF method, the forces are derived from a general energy expression, see Eq. (1)

Esystem ¼ Ebond þ Eover þ Eunder þ Elp þ Eval þ EvdWaals þ ECoulomb

ð1Þ

The partial contributions in Eq. (1) include bond energies (Ebond), energy contributions to penalize over-coordination and (optionally) stabilize under-coordination of atoms (Eover and Eunder), lone-pair energies (Elp), valence angle energies (Eval) and terms to handle non-bonded Coulomb (ECoulomb) and van der Waals (EvdWaals) interaction energies. All terms except the last two include a bond-order dependence and as such depend on the local environment of each atom. The bond-orders are given by a general relationship between bond-order and inter-atomic distance. The bond-order expression includes a set of parameters that are adjusted to fit, e.g., the volume–energy relations (equations of state) for crystals with four, six and eight coordination. The bond-order is then calculated directly from the instantaneous inter-atomic distances, which are continuously updated during dynamics. The Coulomb energy (ECoulomb) of the system is calculated using a geometry dependent charge distribution determined using the electronegativity equalization method (EEM) [21], in which individual atomic charges vary in time. This feature allows ReaxFF to describe charge transfer in chemical reactions. All other non-bonded interactions (short-range Pauli repulsion and long-range dispersion) are included in the van der Waals term (EvdWaals).The non bond-order dependent terms (ECoulomb and EvdWaals) are screened by a taper function and shielded to avoid excessive repulsion at short distances. A detailed description of the individual terms can be found in Refs. [15,16].

2.2. Fitting the ReaxFF parameters for zinc oxide–water systems In the current study, we have re-optimized the force-field parameters for ZnO from our previous study [20], to incorporate recently modified H and O parameters for bulk water systems [22]. The remaining force-field parameters were re-optimized against a training set to ensure the quality of the ZnO description using the new H and O parameters. Included in the training were the same data points as in our previous pure ZnO study, as well as additional data regarding dissociative and molecular water/ ZnO interactions. The training set is a collection of data points (energies, geometries, charges) derived from quantum-mechanical (QM) calculations. The data points were chosen to give adequate descriptions of the volume–energy relations for a number of zinc metal and zinc oxide phases, a number of low-index ZnO surfaces and gas-phase zinc hydroxide clusters. These data points were taken from Ref. [20]. Additional data points were added to the training set (i.e. new QM calculations were performed). These were the water adsorption structures published by Meyer et al. [7]: nine adsorption  0Þ-surface geometries for single water molecules on the ZnOð1 0 1 and dissociation profiles for the molecular (1  1)-M and dissoci-

D. Raymand et al. / Surface Science 604 (2010) 741–752

ated (1  1)-D monolayers and the (2  1)-HD half dissociated monolayer on the same surface. In the study of Meyer et al. [7], the PBE-functional [23] was used. In the present study, the structures from Ref. [7] were reoptimized using the B3LYP-functional [24,25] to remain consistent with previous data for zinc oxide in the ReaxFF training set. The same ZnO basis-set as for the other data points (involving ZnO) in the training-set was used (see Section 2.4.1 for computational details). The re-optimized parameters are the bond parameters for Zn–Zn, Zn–O and the valence angle parameters for Zn–Zn–O, O–Zn–O, O–O–Zn and H–O–Zn. A successive one-parameter search technique was used to minimize the sum of squares error function

Error ¼

2 n  X ðxi;qm  xi;ReaxFF Þ i

ri

ð2Þ

where xi;qm is a QM value, xi;ReaxFF is a ReaxFF calculated value, and ri is a weight assigned to data point i. In total 49 parameters were fitted to a training set containing 380 data points. The result of the reoptimization and a comparison between the reference (QM) data and the model data (ReaxFF) is presented in Section 3.1. 2.3. System descriptions 2.3.1. The systems for the quantum-mechanical calculations  0Þ-surface systems were modeled as slabs, periodic All ZnOð1 0 1 in two dimensions (x and y), with free surfaces in the third dimension (+z and z). The periodic surface cells of the slabs were constructed from the optimized wurtzite bulk cell axes (a = 3.28 Å and c = 5.28 Å) and were then kept fixed in all calculations. In all cases, six layers thick (10 Å) ZnO-slabs were used. The initial guesses for all the H2O/ZnO geometries were taken from the study of Meyer et al. [7]. Ref. [20] contains the system descriptions for the Zn, ZnO and Zn(OH)x data points from our earlier study. Adsorption of single water molecules We performed geometry optimization calculations for single water molecules adsorbed in  nine different configurations, on the ZnO(1010)-surface (i.e. water was adsorbed on one side of the slab). Figs. 1 and 2 show comparisons between the optimized adsorption geometries and energies using B3LYP and the final ReaxFF model, respectively. The structures are presented and labeled in the same way as in Ref. [7]: top and side views of the geometries are labeled as a through i. Calculations for reference states were also performed: a clean  0Þ-surface and an isolated water molecule. The surfaces ZnOð1 0 1 were represented by a Zn36O36-slab (each layer containing six Zn and six O), with a surface area of 9.83  10.55 Å2. Monolayer dissociation paths The following calculations for two different monolayer dissociation paths were performed: (1) a monolayer of molecular water ((1  1)-M) dissociating to a fully hydroxylated surface ((1  1)-D) and (2) a monolayer of molecular water where first every other water molecule dissociates ((2  1)HD) and, second, the remaining water molecules dissociate to form a fully hydroxylated surface. Both paths share common starting and finishing structures, but the second passes through the ordered pattern of alternating dissociated and molecular water molecules found by Meyer et al. as the most favorable structure. Fig. 3 shows the calculated dissociation profiles for B3LYP and ReaxFF along with optimized geometries for the (1  1)-M, (2  1)-HD and (1  1)-D monolayers. The profiles are discussed in Section 3. The dissociation profiles were constructed by linear interpolation of the coordinates of all atoms in nine increments between the three monolayer adsorptions structures. Here, the surfaces were represented by a Zn12O12-slab (each layer containing two Zn and two O) with a surface area of 6.55  5.28 Å2, corresponding to two surface cells.

743

2.3.2. The systems for geometry optimization of stepped surfaces using ReaxFF The influence of a surface step on the adsorption energy of a water monolayer was studied by constructing seven terraced ZnO–wurtzite surface systems with increasing step concentration (on one side of the slab). The steps were parallel to the (0 0 0 1) 0Þ-facets; the surfaces were direction and the terraces were ð1 0 1  0Þ; ð6 1 7  0Þ; in order of increasing step concentration: ð1 0 1  0Þ; ð4 1 5  0Þ; ð2 1 3  0Þ and ð1 1 2  0Þ. For the clean sur 0Þ; ð3 1 4 ð5 1 6 faces (i.e. no water), the surface ions are three-coordinated (including the ions around the step) compared to the four-coordination in the wurtzite structured bulk. To represent the surfaces, 3D-periodic slab-models were used, where the c-axis of the simulation box was elongated to 100 Å making the slabs non-interacting. The slabs were geometry optimized starting from three different geometries: (1) no water, (2) one water molecule per adsorption site (molecular adsorption) and (3) one dissociated water per adsorption site (dissociative adsorption). The adsorption site is here defined as one zinc and surface oxygen pair on the slab surface. We also performed geometry optimizations for two reference  0Þ-surface with cases: (1) a half dissociated monolayer on the ð1 0 1 the (2  2)-HD pattern and (2) an isolated water molecule. Fig. 4 shows the optimized geometries of the water covered (molecular or dissociated) stepped surfaces along with their corresponding adsorption energies plotted as a function of step concentration. For a closer view of the water geometry on the terraces and at the steps, see Figs. 3 and 5, respectively. A detailed description of the systems, including the box dimensions, is given in Table 1. 2.3.3. The systems for MD simulations using ReaxFF The geometry optimized structures were taken as initial geometries (see Section 2.3.2). 2.4. Computational details 2.4.1. Quantum-mechanical calculations The QM-calculations that provided the surface data used to parametrize ReaxFF were performed using the CRYSTAL06 program [26]. All calculations for the ZnO surfaces were performed using 2-dimensionally periodic systems (as mentioned), Gaussian type basis-sets, and the B3LYP hybrid DFT functional as implemented in the CRYSTAL06 program [26]. A basis-set optimized by Jaffe and Hess [27] was used for the zinc and oxygen atoms of ZnO. This basis-set has been applied with good results in several earlier QM-studies of ZnO polymorphs and surfaces [28–30]. The 6-311+G* basis-set [31–33] was used for the H and O atoms of the water molecules. This basis-set for water was chosen to remain consistent with previous data used to parametrize ReaxFF. The Brillouin zone was in all cases sampled with a 4  4 k-point mesh of the Monkhorst-Pack type. During the geometry optimization of all the slabs, the bottom three layer ion positions were fixed at their corresponding bulk positions (following the study of Meyer et al.) The top three layers including the water molecule (when present) were allowed to fully relax. All structures were optimized using the built-in geometry optimizer of CRYSTAL06 with default parameters. The force-field training set energies were computed by taking the energy difference between the adsorption structure and the reference structures (a clean slab and an isolated water molecule). The dissociation profiles were constructed by optimizing the end point geometries (i.e. (1  1)-M, (2  1)-HD, (1  1)-D) and linearly interpolating the coordinates between in nine increments (as mentioned above). These increments (27 in total) were not geometry-optimized, and therefore only give an estimation of the dissociation profile. The energies entered in the force-field training-set were computed by using the molecularly adsorbed mono-

744

D. Raymand et al. / Surface Science 604 (2010) 741–752

a) B3LYP

b) B3LYP

c) B3LYP

d) B3LYP

e) B3LYP

a) ReaxFF

b) ReaxFF

c) ReaxFF

d) ReaxFF

e) ReaxFF

h) B3LYP

i) B3LYP

f) B3LYP

f) ReaxFF

g) B3LYP

g) ReaxFF

h) ReaxFF

i) ReaxFF

 0Þ-surface using B3LYP and ReaxFF. The side view shows the top four layers (out of Fig. 1. The adsorption geometries for a single water molecule optimized on the ZnOð1 0 1 six), and the top view shows the top two layers. Zn, O and H are represented by small gray, large red and small white spheres respectively. To guide the eye, the O in the adlayer are blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

layer as reference and calculating the differences with the other structures along the dissociation profile. The energies were corrected for basis-set superposition errors using the counterpoise method [34]. Thus, all structures in Fig. 1

except geometries e, f and i, were divided into two fragments: the ZnO-slab and the water molecule. Structures e, f and i were divided into three fragments: the ZnO-slab, H and OH. Mulliken charges were used to evaluate the OH character, radical or ionic. In all cases

745

D. Raymand et al. / Surface Science 604 (2010) 741–752

PBE B3LYP

-5

ReaxFF-

kcal/mol

-10 -15 -20

a

b

c

d

e

f

g

h

i

 0Þ-surface using B3LYP and ReaxFF. Adsorption energies from Ref. [7] calculated with Fig. 2. The adsorption energies for a single water molecule optimized on the ZnOð1 0 1 PBE are given for reference.

(a) 5

_ ZnO(1010)-(1x1)-M

B3LYP

Energy/(kcal/mol)

4 3 2 1 0

_ ZnO(1010)-(2x1)-HD

_ ZnO(1010)-(1x1)-D

-1 -2

arbitrary reaction coordinate

_ ZnO(1010)-(1x1)-M

(b)

ReaxFF

5

Energy/(kcal/mol)

4 3 2 1 0

_ ZnO(1010)-(2x1)-HD

_ ZnO(1010)-(1x1)-D

-1 -2

arbitrary reaction coordinate  0Þ using (a) B3LYP and (b) ReaxFF along with optimized geometries for the fully molecular Fig. 3. The calculated dissociation profiles for a water monolayer on ZnOð1 0 1 (1  1)-M, half dissociated (2  1)-HD and fully dissociated (1  1)-D monolayers. Zn, O and H are represented by small gray, large red and small white spheres, respectively. To guide the eye, the O in the ad-layer are blue and the periodic images of the simulation box show in flat gray. The two traces correspond to a direct transition from molecular to dissociated adsorption (green curve), and a transition via the half-dissociated structure (red curve). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

mostly OH-character was found. For the BSSE calculations of the monolayer dissociation paths, the structures were divided into two

fragments, the ZnO-slab and the water monolayer. The monolayer was in turn divided into two water molecules.

746

D. Raymand et al. / Surface Science 604 (2010) 741–752

_

_ (6170)

_ (5160)

_ (4150)

_ (3140)

_ (2130)

_ (1120)

-26

Eads (kcal/mol)

Dissociative

Molecular

Dissociative

Molecular

Dissociative

Molecular

Dissociative

Molecular

(1010)

Molecular

-27 -28 Dissociative

-29

Mixed molecular/ dissociative on (1010) 0

0.02

0.04

0.06

0.08 0.10

0.12

0.14

0.16

Steps/Å Fig. 4. The optimized structures of the (molecularly or dissociated) water-covered stepped surfaces, along with their corresponding adsorption energies (in kcal per mole of adsorbed water molecules – dissociated or not) plotted as a function of step concentration. All stepped surfaces are seen to prefer dissociation. All surfaces are viewed along the (0001) direction. The adsorption energy of the mixed molecular-dissociated (2  1)-HD monolayer is included to show the effect of favorable hydrogen bonding. Zn, O (of ZnO), O (from the water molecules), and H are represented by small gray, large red, large blue, and small white spheres, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Our earlier study, Ref. [20], contains the computational details about the Zn, ZnO and Zn(OH)x data points. 2.4.2. Geometry optimizations using ReaxFF The geometries were optimized using the built-in conjugate gradient minimizer of the ReaxFF software package with a convergence criterion of 0.01 kcal/mol between two consecutive

optimization steps. During the geometry optimization the volume of the simulation box was kept fixed. The energies reported in Fig. 4 were calculated by taking the energy difference between the water covered surface and the sum of the clean surface and corresponding number of isolated water molecules (i.e. the number of adsorption sites per simulation box).

747

D. Raymand et al. / Surface Science 604 (2010) 741–752

Molecular

Dissociative

 0)-surface. The energy difference per water between the molecular and Fig. 5. Detailed side and top view of the adsorption structure of water at step sites on the ZnO(1 1 2 dissociated cases is 1.9 kcal/mol. Zn, O (of ZnO), O (from the water molecules), and H are represented by small gray, large red, large blue, and small white spheres, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 1 Slab models used with ReaxFF to study stepped surfaces. Values are given per repeating unit (simulation box).

a b

Surface

Number of adsorption sites per simulation boxa

Number of ZnO formula units per simulation box

Surface cell area (Å2)

Slab thicknessb (Å)

Number of steps per simulation box

Step concentration (step/Å)

 0Þ ð1 0 1  0Þ ð6 1 7  0Þ ð5 1 6  0Þ ð4 1 5  0Þ ð3 1 4  0Þ ð2 1 3  0Þ ð1 1 2

48 56 48 40 64 48 64

288 312 280 240 320 320 320

19.9  21.6 21.6  21.7 21.4  18.5 21.5  15.2 21.5  24.0 21.4  17.5 21.4  22.9

14.9 15.3 15.2 17.7 11.6 19.6 14.6

0 1 1 1 2 2 4

0.000 0.046 0.054 0.066 0.083 0.114 0.174

An adsorption site corresponds to one surface Zn2+-ion and one surface O2-ion. Water is only adsorbed on one side of the slab. The slab thickness is calculated as the difference in z-coordinate between topmost and bottommost ions of the optimized slabs.

2.4.3. Molecular dynamics simulations All simulations were performed using the ReaxFF software package. The simulations reported in this paper were carried out with constant number of atoms, constant volume and constant temperature (NVT). The velocity Verlet [35] integrator was used with a timestep of 0.25 fs. The Berendsen thermostat [36] with a temperature damping constant of 0.1 ps was used. For each surface, the simulations were run for 100 ps at three different temperatures (10 K, 300 K and 600 K), yielding a total of 42 simulations (seven surfaces with either a molecularly or dissociatively adsorbed monolayer as initial guess). 3. Results and discussion 3.1. Parametrization of the ReaxFF force-field The re-optimization of the ZnO ReaxFF force-field resulted only in a small change of the original ZnO parameters. A comparison between data from QM calculations, experiment and the current and previous ReaxFF models is given in Tables 2–4. The properties tested in Ref. [20] remain essentially the same. The

four ZnO polymorphs are seen to be reproduced in the correct relative order of stability by the ReaxFF model. Moreover, the model gives a good description of elastic properties of wurtzite structured ZnO. The relative order of the surface energies is reproduced for the five fitted ZnO-facets. As seen in Tables 2 and 3, the new ReaxFF model, which includes H2O reactions, retains the accuracy for describing both the ZnO-only and Zn-only systems obtained in Ref. [20]. For the ZnO surfaces, the new ReaxFF parameters improve the accuracy over the previous ReaxFF force-field, as seen in Table 4. 3.1.1. Adsorption energies for single water molecules The present ReaxFF description reproduces the QM results adequately for the adsorption of single water molecules, both energies and geometries (Figs. 1 and 2). The worst agreement is for cases d and e, where the ReaxFF leads to less stability than the QM. This may be because ReaxFF finds these to be geometrically similar to geometries a, b and c. In optimizing the Zn–O parameters there is a competition between reproducing these structures and the (1  1)-M-adsorption structure; the result is a compromise.

748

D. Raymand et al. / Surface Science 604 (2010) 741–752

Table 2 QM results and ReaxFF results (at 0 K) compared to experimental data (at room temperature) for the cell axes, heat of formation ðDf HÞ, bulk moduli and elastic constants of the four polymorphs (B1–B4) of ZnO. Property

B3LYP (Ref. [20])

ReaxFF (this work)

ReaxFF (Ref. [20])

Experiment

3.29 5.30 91.2 144 222.9 116.3 103.5 212.8 57.1

3.25 [37] 5.21 [37] 83.3 [38] 141, 143 [27,39] 209.7 [40] 121.1 [40] 105.1 [40] 210.9 [40] 42.47 [40]

130 4.62 0.97

Polymorph does not occur in nature

Wurtzite(B4) P63 mc a (Å) c (Å) Df HB4 (kcal/mol) Bulk modulus (GPa) c11 (GPa) c12 (GPa) c13 (GPa) c33 (GPa) c44 (GPa)

3.28 5.28

3.28 5.28 84.4 139 207.9 113.3 81.4 287.5 52.0

Bulk modulus (GPa) a (Å) Df HB1 —Df HB4 ðkcal=molÞ

162 4.60 0.50

130 4.60 1.19

Bulk modulus (GPa) a (Å) Df HB3 —Df HB4 ðkcal=molÞ

202 4.30 8.57

227 4.26 8.15

283 4.44 7.72

203, 228 [39,41] 4.27 [41]

Bulk modulus (GPa) a (Å) Df HB2 —Df HB4 ðkcal=molÞ

183 2.68 37.89

 Caesium chloride (B2) Pm3m 327 2.63 43.24

407 2.64 37.58

Polymorph does not occur in nature

136

 Zincblende(B1) F 43m

 Rocksalt(B3) Fm3m

Table 3 QM results and ReaxFF results (at 0 K) compared to experimental data (at room temperature) for the cell axes, cohesive energy ðEcoh Þ, bulk moduli and elastic constants of Znmetal. Property

PBE (Ref. [20])

ReaxFF (this work)

ReaxFF (Ref. [20])

Experiment

2.73a 4.46a 31.6 87.7 193 41 42 188c 71c

2.67 [42] 4.95 [42] 31.1 [43] 64.5–75.1b 163b 31b 48b 60c 39c

87.7 3.86 0.00a

Polymorph does not occur in nature

73.5 3.06 2.71

Polymorph does not occur in nature

30.2 2.71 6.62

Polymorph does not occur in nature

hcp P63/mmc a (Å) c (Å) Ecoh (hcp) (kcal/mol) Bulk modulus (GPa) c11 ðGPaÞ c12 ðGPaÞ c13 ðGPaÞ c33 ðGPaÞ c44 ðGPa

2.63 5.06 66.7

2.74a 4.49a 31.0 110 224 56 53 253c 74c  fcc F 43m

Bulk modulus (GPa) a (Å) Ecoh (fcc)–Ecoh (hcp) (kcal/mol)

81.8 3.86 1.31

78 3.87 0.00a

Bulk modulus (GPa) a (Å) Ecoh (bcc)–Ecoh (hcp) (kcal/mol)

84.6 3.06 2.71

89 3.11 2.15

Bulk modulus (GPa) a (Å) Ecoh (sc)–Ecoh (hcp) (kcal/mol)

64.2 2.71 6.00

30 2.77 6.11

 bcc Im3m

 sc Pm3m

a Valence angle terms for pure Zn were not included; therefore ReaxFF will preserve the ideal c/a ratio ((8/3)1/2) of a close-packed crystal. Also the volume–energy relationship for the fcc and hcp phases will be identical. b Experimental values taken from the review article by Ledbetter [44]. c These hcp-specific stresses are not correctly reproduced for the reason stated in a.

3.1.2. Monolayer dissociation paths The ReaxFF model is able to reproduce the correct order of stability between the three adsorption patterns ((1  1)-M, (1  1)-D and (2  1)-HD) and the barriers are similar (and small), see Fig. 3. Fully dissociative adsorption ((1  1)-D) is found to be least favorable because the energy loss from the dissociation is not compensated by any additional hydrogen bonds. Compared to the B3LYP reference data, the (1  1)-M-adsorption structure is more stable, although the difference is small

(1 kcal/mol). The method of linearly interpolating the optimized end structures, i.e. not optimizing the nine increments between the optimized end structures, gives an upper estimate of the dissociation barrier height. The actual barrier during dynamic simulation is likely to be lower. The final force-field parameters are presented in Tables 5–10. The values of the bond- and valence angle parameters indicate that the ReaxFF model describes the ZnO lattice to have a substantial covalent character, including relatively high angular force constants for the O–Zn–O and Zn–O–Zn angles. For the H–Zn

749

D. Raymand et al. / Surface Science 604 (2010) 741–752 Table 4 QM results versus ReaxFF model results for surface energies of the five ZnO surfaces present in the ReaxFF training set. B3LYP (Ref. [20]) ReaxFF (this work) ReaxFF (Ref. [20]) Surface energy Surface energy Surface energy (J m2) (J m2) (J m2)  0Þ Wurtzite ð1 0 1  0Þ Wurtzite ð1 1 2 Zincblende (1 0 0) Rocksalt (1 0 0) Rocksalt (1 1 0)

1.32 1.39 1.30 2.00 1.16

1.31 1.45 1.37 1.68 1.28

0.96 1.06 1.00 1.10 0.73

Table 5 The final ReaxFF atom parameters for Zn. For O and H, parameters from Ref. [22] was used. Definitions of the individual ReaxFF parameters in this table and Tables 6–10 can be found in Refs. [15,16]. r0

H O Zn

H O Zn

pov/un

0.8930 19.4571 1.2450 3.5500 1.8862 3.0614 van der Waals parameters

Table 8 The final ReaxFF valence angle parameters for H–O–Zn, O–Zn–O, Zn–O–Zn, O–Zn–Zn and O–O–Zn. Parameters for H–H–H, O–O–O, H–O–O, H–O–H, O–H–O and H–H–O are taken from Ref. [22]. Valence angle

H0;0 (deg)

ka (kcal/ mol)

kb (1/ radian)2

pv,1

pv,2

H–H–H O–O–O H–O–O H–O–H O–H–O H–H–O H–O–Zn O–Zn–O Zn–O–Zn O–Zn–Zn O–O–Zn

0.0000 80.7324 75.6935 85.8000 0.0000 0.0000 77.5446 10.4094 37.8790 16.9624 60.2000

27.9213 30.4554 50.0000 9.8453 15.0000 8.5744 9.9016 38.9915 32.3525 30.3241 20.0000

5.8635 0.9953 2.0000 2.2720 2.8900 3.0000 2.3157 0.7072 0.2657 0.2697 0.5000

0.0000 1.6310 1.0000 2.8635 0.0000 0.0000 0.4543 2.0000 0.4403 2.0000 1.0000

1.0400 1.0783 1.1680 1.5800 2.8774 1.0421 2.3770 2.6162 1.1000 3.0708 2.0000

H0;0 is (180° – equilibrium valence angle).

Coulomb parameters

g (eV)

v (eV)

c (Å1)

9.6093 8.3122 5.7915

3.7248 8.5000 2.0219

0.8203 1.0898 0.4828

rvdW (Å)

a

cvdW (Å1)

 (kcal/mol)

1.3550 2.3890 1.9200

8.2230 9.7300 11.5134

33.2894 13.8449 18.3776

0.0930 0.1000 0.2998

Table 6 The final ReaxFF van der Waals and bond radius parameters for the H–Zn bond and the O–Zn bond. Parameters from Ref. [22] was used for the H–O bond. 1

Bond

rr ðÅÞ

rvdW (Å)

 (kcal/mol)

cvdW (Å )

H–O H–Zn O–Zn

0.9215 0.1000 1.9804

1.2885 1.8227 2.1414

0.0283 0.0987 0.2744

10.9190 12.0654 9.7703

Table 7 The final ReaxFF bond energy and bond-order parameters for the Zn–Zn and Zn–O bonds. The H–Zn is defined as a zero covalent interaction. The H–H, H–O and O–O parameters from Ref. [22] were used. Bond

Dre ðkcal=molÞ

pbe,1

pbe,2

pbo,1

pbo,2

pcov

H–H H–O O–O H–Zn O–Zn Zn–Zn

153.3936 160.0000 142.2858 0.0000 159.9755 38.4643

0.4600 0.5725 0.2605 0.0000 0.4548 0.6944

6.2500 1.1150 0.3451 0.5000 1.3099 0.5059

0.0970 0.0920 0.1225 0.2000 0.4787 0.0814

6.0552 4.2790 5.5000 8.0000 4.6717 6.0333

0.7300 0.5626 0.6051 0.5000 0.0375 0.2129

interaction, no bond parameters were derived, i.e. it was only considered as non-bonded interaction, which explains the low sigmabond radius (0.1 Å) and dissociation energy (0.00 kcal/mol). In principle parameters for H–Zn bonds could be included into the ReaxFF description, but since these interactions do not play a major part in the ZnO/water interaction they were excluded.

3.2. H2O structure on stepped ZnO surfaces One monolayer of H2O was adsorbed on the seven different ZnO faces (see Table 1 for step concentration) and the geometries were optimized. As mentioned above, in all cases, the optimization was started from two geometries; fully molecular and fully dissociated.

Table 9 The ReaxFF torsion angle parameters taken from Ref. [22]. Torsion angle

V1 (kcal/ mol)

V2 (kcal/ mol)

V3 (kcal/ mol)

ptor,1

pcot,1

H–O–O–H H–O–O–O O–O–O–O

2.2500 0.4723 2.5000

6.2288 12.4144 25.0000

1.0000 1.0000 1.0000

2.6189 2.5000 2.5000

1.0000 1.0000 1.0000

Table 10 The ReaxFF hydrogen bond parameters taken from Ref. [22]. Bond

rhb (Å)

phb,1

phb,3

phb,3

O–H–O

2.1200

3.5800

1.4500

19.5000

The optimized geometries and adsorption energies (Eads per mole of water molecules) are shown in Fig. 4 as a function of step concentration. For reference, the value of the mixed molecular/disso 0Þ is included ciative (2  1)-HD adsorption pattern on ZnOð1 0 1 in the graph to show the effect of a more favorable hydrogen-bond pattern. Both molecular and dissociative adsorption show an almost linear dependence on step concentration. Molecular adsorption is seen to be more stable than dissociative on the flat  0Þ-surface. In all other cases dissociative adsorption is ZnOð1 0 1 more favorable. Our findings are consistent with the conclusions  0Þ-surface that fully dissociaof Cooke et al. [10] for the ZnOð1 1 2 tive adsorption is favored because of a stabilizing hydrogen-bond network that is not present in the fully molecular case. The balance between molecular/dissociative adsorption on the various ZnO surfaces is governed by hydrogen bonding. This balance is delicate since the hydrogen-bonded network is disturbed by the presence  0Þ), of a step (even in the case were the terraces are largest ð6 1 7 and the energy difference between molecular and dissociative adsorption increases. MD simulations were performed for all surfaces with T = 10 K, T = 300 K and T = 600 K. The results are presented in Fig. 6, which shows the coverage of the surface adsorbates (H2O or OH) as a function of time, starting from a fully molecular water monolayer as well as from a fully dissociated monolayer. The figures also give the total coverage (H2O plus OH), where any desorbed water molecules are excluded. Common for all our surface simulations is that equilibrium is not reached with respect to molecular/dissociative adsorption at T = 10 K. For the simulations at T = 300 K and T = 600 K, equilibrium is reached regardless of the initial adsorption configuration – molecular or dissociative. For both T = 300 K

D. Raymand et al. / Surface Science 604 (2010) 741–752

T=10K

T=600K

_ (6170)

100 80 60 40 20 0

(b)

Time(ps)

Time(ps)

(c)

0 20 40 60 80 0 20 40 60 80 0 20 40 60 80 100

100 80 60 40 20 100 80 60 40 20 0

(d)

0 20 40 60 80 0 20 40 60 80 0 20 40 60 80 100

coverage(%)

100 80 60 40 20 0

T=600K

_ (2130)

Dissociative

(e)

100 80 60 40 20

T=300K

Molecular

Dissociative

100 80 60 40 20 0

T=10K Molecular

coverage(%)

Time(ps) T=600K

100 80 _ 60 (3140) 40 20

(f)

0 20 40 60 80 0 20 40 60 80 0 20 40 60 80 100

0 20 40 60 80 0 20 40 60 80 0 20 40 60 80 100

Time(ps)

Time(ps)

T=10K

100 80 60 40 20

(1120)

100 80 60 40 20 0

(g)

T=300K

T=600K Molecular

_

Dissociative

coverage(%)

T=300K

T=600K

_ (4150)

Time(ps) T=10K

T=300K

Dissociative

Dissociative

100 80 60 40 20 0

coverage(%)

_ (5160)

T=10K Molecular

coverage(%)

T=600K

Molecular

100 80 60 40 20

T=300K

T=600K

0 20 40 60 80 0 20 40 60 80 0 20 40 60 80 100

0 20 40 60 80 0 20 40 60 80 0 20 40 60 80 100

T=10K

T=300K

Dissociative

(a)

100 80 60 40 20

Molecular

_ (1010)

Dissociative

100 80 60 40 20 0

T=300K

Molecular

coverage(%)

T=10K

100 80 60 40 20

coverage(%)

750

0 20 40 60 80 0 20 40 60 80 0 20 40 60 80 100

Time(ps) Fig. 6. The time variation of the coverage of the surface adsorbates (H2O or OH) during the MD simulations, starting from a fully molecular water monolayer and a fully dissociated monolayer at three different temperatures (10 K, 300 K and 600 K). The figures also show the total coverage (H2O plus OH), where any desorbed water molecules are excluded. The color scheme for the graph is: H2O (green), OH (red) and total (blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

and T = 600 K, the equilibrium is seen to be reached quickly (in less than 20 ps). At T = 600 K, water desorbs from the surface and the coverage falls. The T = 10 K simulations show very few features (with respect to dissociation and re-association): only three cases

started from a molecular configuration display some features (see Fig. 6b,e and f). In these cases water is able to dissociate close to a step. Thus, since the system is able to cross the barrier for dissociation already at T = 10 K, we conclude that it is very low close

751

D. Raymand et al. / Surface Science 604 (2010) 741–752

(a)

T=300K

(b)

80 60 40 dissociative start molecular start

20 0

T=600K

100

0

0.04

0.08 0.12 steps/Å

0.16

OH-coverage (%)

OH-coverage (%)

100

80 60 40 dissociative start molecular start

20 0

0

0.04

0.08 0.12 steps/Å

0.16

Fig. 7. The fraction of OH-coverage plotted with respect to step concentration for the different stepped ZnO surfaces, at (a) 300 K and (b) 600 K. Red squares represents OHcoverage starting from fully dissociative start, green circles a fully molecular start. The values plotted are the average fraction of OH-coverage during the last 12.5 ps (1000 MD snapshots) from each simulation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the T = 600 K simulation; the hydroxylation of the step is larger than the hydroxylation of the terraces. The MD results are consistent with the findings from the geometry optimizations: steps favor hydroxylation because the oxide structure can accommodate the formation of a more favorable Hbond network and a stronger adsorption. On the terraces, mixed adsorption is favored. 4. Conclusions  Fig. 8. ZnOð6170Þ-surface (steps/Å = 0.046). Snapshot from the end of the T = 600 K simulation. The step shows an increased hydroxylation. Zn, O (of ZnO), O (from the water molecules), and H are represented by small gray, large red, large blue, and small white spheres, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

to a step. We observe no occurrence of the reverse reaction, i.e. reassociation of two hydroxyl groups (Olattice–H and Znlattice–OH) to molecular water.  0Þ-surface at T = 300 K, the results agree with For the ZnOð1 0 1 the predictions of Meyer et al. [7], i.e. an ordered (2  1)-HD pattern is observed, where half of the water molecules are dissociated (see Fig. 6a). At T = 600 K, water molecules evaporate and the coverage falls, but the 1:1 water–hydroxyl ratio is maintained. During the 100 ps of simulation, the coverage falls to 60%, i.e. a 30% hydroxyl and 30% water coverage. The reason for this, as stated before, is that it is possible to form a more favorable hydrogen-bonded network in a mixed molecular/dissociated monolayer. The effect is seen to be greater than 1 kcal/mol.  0Þ; ð4 1 5  0Þ; ð5 1 6  0Þ; On the stepped surfaces, ð6 1 7  0Þ; ð2 1 3  0Þ and ð1 1 2  0Þ, the level of hydroxylation increases ð3 1 4 with increasing step density, as seen in Fig. 6. At both T = 300 K and T = 600 K, the increase is visible, but it is more pronounced at T = 600 K, where water molecules have started to desorb, thus shifting the water–hydroxyl ratio towards hydroxylation on the surface, since only molecular water is able to evaporate. This is perhaps most apparent in Fig. 6g, where the level of hydroxylation  0Þincreases, in spite of a decrease in total coverage. On the ð1 0 1 terraces, the 1:1 water:hydroxyl equilibrium is maintained during desorption, while the steps remain mostly hydroxylated. To quantify this, the level of hydroxylation with respect to step concentration was plotted (taken as the average fraction of OH-coverage during the last 12.5 ps from each simulation), see Fig. 7. The plot confirms that the number of hydroxylated surface sites is mainly a function of step concentration since the level of hydroxylation shows an almost linear dependence on step concentration. Fig. 8  0Þ-surface at the end of shows a snapshot of a step on the ð6 1 7

A good fit of the ReaxFF model to the QM reference data for the ZnO–water system was achieved. Monolayer structures that promote hydrogen bonding are favored. Steps on the surface promote an increased hydroxylation, while on terraces a dynamic equilibrium between water and hydroxyl groups is present. An increased temperature leads to water desorption on the terraces leaving hydroxylated surface steps. The presence of hydroxyl groups around steps is likely an important effect to take into account in the context of catalytic processes on ZnO surfaces. The current ZnO/H2O force-field combined with ReaxFF’s existing capability to model hydrocarbons is an important step in our long term goal to directly model metal/metal oxide interface catalysis. Acknowledgments This work was supported by the Swedish Research Council (VR) and by NSF Grant #DMR 0427177. Computer time was provided by SNIC. The computations were performed on UPPMAX under Project s00707-54. Many useful discussions with Dr. Tomas Edvinsson and Dr. Pavlin Mitev are acknowledged. References [1] G.A. Somorjai, Introduction to Surface Chemistry and Catalysis, John Wiley & Sons, Inc, New York, USA, 1994. [2] V. Henrich, A. Cox, The Surface Science of Metal Oxides, Cambridge University Press, Cambridge, UK, 1994. [3] D. Scarano, G. Spoto, S. Bordiga, A. Zecchina, C. Lamberti, J. Am. Chem. Soc. 276 (1992) 281. [4] C. Wöll, Prog. Surf. Sci. 82 (2007) 55. [5] B. Meyer, D. Marx, O. Dulub, U. Diebold, D. Langenberg, C. Wöll, Angew. Chem. Int. Ed. 43 (2004) 6642. [6] O. Dulub, B. Meyer, U. Diebold, Phys. Rev. Lett. 95 (2005) 136101. [7] B. Meyer, H. Rabaa, D. Marx, Phys. Chem. Chem. Phys. 8 (2006) 1513. [8] L. Giordano, J. Goniakowski, J. Suzanne, Phys. Rev. Lett. 81 (1998) 1271. [9] M. Odelius, Phys. Rev. Lett. 82 (1999) 3919. [10] D.J. Cooke, A. Marmier, S. Parker, J. Phys. Chem. B 110 (2006) 7985. [11] D. Raymand, T. Edvinsson, D. Spångberg, A. van Duin, K. Hermansson, Proc. SPIE 7044 (2008) 70440E. [12] C. Catlow, S. French, A. Sokol, A. Al-Sunaidi, S. Woodley, J. Computat. Chem. 29 (2008) 2234. [13] M. Henderson, Surf. Sci. Rep. 46 (2002) 1.

752

D. Raymand et al. / Surface Science 604 (2010) 741–752

[14] D. Costa, K. Sharkas, M.M. Islam, P. Marcus, Surf. Sci. 603 (2009) 2484. [15] A.C.T. van Duin, S. Dasgupta, F. Lorant, W.A. Goddard III, J. Phys. Chem. A 105 (2001) 9396. [16] A.C.T. van Duin, A. Strachan, S. Stewman, Q. Zhang, X. Xu, W.A. Goddard III, J. Phys. Chem. A 107 (2003) 3803. [17] W.A. Goddard III, A.C.T. van Duin, K. Chenoweth, M. Cheng, S. Pudar, J. Oxgaard, B. Merinov, Y.H. Jang, P. Persson, Top. Catal. 38 (2006) 1. [18] K. Chenoweth, A.C.T. van Duin, P. Persson, M.-J.J.O. Cheng, W.A. Goddard III, J. Phys. Chem. C 112 (37) (2008) 14645. [19] Q. Zhang, T. Cagin, A.C.T. van Duin, W.A. Goddard III, Y. Qi, G. Hector, Phys. Rev. B 69 (2004) 045423. [20] D. Raymand, M. Baudin, A. van Duin, K. Hermansson, Surf. Sci. 602 (2008) 1020. [21] G.O.A. Janssens, B.G. Baekelandt, H. Toufar, W.J. Martier, R.A. Shoonheydt, J. Phys. Chem. B 99 (1995) 3251. [22] A.C.T. van Duin, V. Bryantsev, W.A. Goddard III, in preparation. [23] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [24] A.D. Becke, J. Chem. Phys. 98 (1993) 5968. [25] C. Lee, W. Yang, R. Parr, Phys. Rev. B 37 (1988) 785. [26] R. Dovesi, V. Saunders, C. Roetti, R. Orlando, C. Zicovich-Wilson, F. Pascale, B. Civalleri, K. Doll, N. Harrison, I. Bush, P. D’Arco, M. Llunell, CRYSTAL06 User’s Manual, 2006. [27] J.E. Jaffe, A.C. Hess, Phys. Rev. B 48 (1993) 7903.

[28] J. Jaffe, A. Hess, J. Chem. Phys. 104 (1996) 3348. [29] A. Wander, N. Harrison, Surf. Sci. 457 (2000) 342. [30] Y. Noel, C.M. Zicovich-Wilson, B. Civalleri, P. D’Arco, R. Dovesi, Phys. Rev. B 65 (2001) 014111. [31] R. Krishnan, J.S. Binkley, R. Seeger, J.A. Pople, J. Chem. Phys. 72 (1980) 650. [32] T. Clark, J. Chandrasekhar, G.W. Spitznagel, P.V.R. Schleyer, J. Comp. Chem. 4 (1983) 294. [33] K. Raghavachari, G.W. Trucks, J. Chem. Phys. 91 (1989) 1062. [34] S. Boys, F. Bernardi, Mol. Phys. 19 (4) (1970) 553. [35] W. Swope, H. Andersen, P. Berens, K. Wilson, J. Chem. Phys. 76 (1982) 637. [36] J.C. Berendsen, J. Chem. Phys 81 (1984) 3684. [37] J. Albertsson, S.C. Abrahams, Å. Kvick, Acta Crystallogr. Sec. B 45 (1) (1989) 34. [38] W.L. Masterton, E.J. Slowinski, C.L. Stanitski, Chemical Principles, CBS, New York, USA, 1987. [39] S. Desgreniers, Phys. Rev. B 58 (1998) 14102. [40] T. Bateman, J. Appl. Phys. 33 (1962) 3309. [41] W.P.H. Karzel, M. Kofferlein, W. Schiessl, M. Steiner, U. Hiller, G.M. Kalvius, D.W. Mitchell, T.P. Das, P. Blaha, K. Schwarz, M.P. Pasternak, Phys. Rev. B 58 (1998) 14102. [42] J. Brown, J. Inst. Metals 83 (1955) 49. [43] G. Aylward, T. Findlay, S.I. Chemical Data, John Wiley & Sons, Ltd., Milton, Australia, 1987. [44] H. Ledbetter, J. Phys. Chem. Ref. Data 6 (1977) 1181.