An SCC-DFTB Repulsive Potential for Various ZnO Polymorphs and the ZnO-Water System.

We have developed an efficient scheme for the generation of accurate repulsive potentials for self-consistent charge density-functional-based tight-binding calculations, which involves energy-volume scans of bulk polymorphs with different coordination numbers. The scheme was used to generate an optimized parameter set for various ZnO polymorphs. The new potential was subsequently tested for ZnO bulk, surface, and nanowire systems as well as for water adsorption on the low-index wurtzite (101̅0) and (112̅0) surfaces. By comparison to results obtained at the density functional level of theory, we show that the newly generated repulsive potential is highly transferable and capable of capturing most of the relevant chemistry of ZnO and the ZnO/water interface.

Introduction

Zinc oxide is a wide band gap semiconductor (band gap 3.4 eV), used in several technologically important applications, such as heterogeneous catalysis,[1] gas sensors,[2] and microelectronic devices.[3] Many of these applications have in common that they rely on the specific electronic properties of ZnO and how they can be modified by dopants or by specific structural features. In particular, ZnO nanoparticles of various sizes can be grown to exhibit a large number of different shapes, such as wires, spheres, and helices.[2] ZnO nanoparticles can also be used in, for example, sunscreens to protect against UV irradiation, but the toxicity of the nanoparticles is debated[4,5] and may be related to the dissolution of Zn2+ ions.

The interaction of ZnO with water or OH-groups is of particular interest, because water in either the liquid or the gas phase tends to be present for most of the applications involving ZnO. At ambient conditions, ZnO exhibits the wurtzite crystal structure and principally exposes four different surfaces: the nonpolar ZnO(101̅0) and ZnO(112̅0) and the polar Zn-terminated ZnO(0001) and O-terminated ZnO(0001̅).[6] Adsorbed water adopts different structures for different surfaces. On ZnO(101̅0), experiments and calculations have shown that the most favorable adsorption structure for water is “half-dissociated”, where every other water molecule is dissociated.[7] This leads to a particularly stable hydrogen-bonded network on the surface. However, the half-dissociated adsorption structure coexists with a molecular adsorption structure on the surface,[8] and calculations suggest that the amount of dissociated water may increase upon addition of more water layers.[9] On ZnO(112̅0), the situation is to a large extent unknown, although it has been proposed by density-functional theory (DFT) calculations that water adsorbs either fully dissociated[10] or half-dissociated.[11] The clean polar Zn-terminated ZnO(0001) surface exhibits triangular pits[12] to stabilize the inherent polarity in the ZnO[0001] direction. On this surface, water adsorbs dissociatively and reduces the pit sizes.[13] Finally, the O-terminated ZnO(0001̅) surface is usually hydrogen covered under normal preparation conditions. Adsorption of water on this hydrogen-covered surface has been proposed to be molecular, while adsorption on the clean (non-hydrogen-covered) surface has been proposed to be dissociative.[14]

Recently, there has also been interest in crystal structures of ZnO other than wurtzite, such as the closely related cubic zincblende structure,[15] as well as high-pressure modifications such as the NaCl-type and CsCl-type structures[16−19] and the lesser-known “graphitic”[20−25] and low-density “body-centered tetragonal”[26] (BCT)[25,27−32] and cubane-type[33] structures. The graphitic and BCT structures have been proposed to form during thin-film growth of ZnO to quench the macroscopic dipole moment formed as wurtzite ZnO grows along the polar ZnO[0001]/ZnO[0001̅] directions. Each of these polymorphs have their own set of lattice parameters, formation energies, and electronic properties that may make them more suitable than the wurtzite structure for certain applications. Currently, all of the polymorphs except the CsCl, BCT, and cubane structures have been synthesized in experimental laboratories. However, a BCT-like structure has been shown to form dynamically at the ZnO(101̅0) surface.[32]

With this palette of ZnO structures and their individual characters, it is highly desirable to have access to a set of theoretical methods which would allow for the study of different polymorphs of ZnO and their interaction with water at experimentally relevant size and time scales, concerning both structural and electronic properties. Given the significant computational cost associated with very large-scale and simultaneously accurate theoretical calculations at the ab initio level today, other methods, which bridge the existing size and time gaps between experiments and theory, are needed.

Such methods are invariably parametrized or semiempirical, i.e., they rely on precalculated interaction parameters between different atomic species. The parameters are obtained by fitting them to the results of a set of reference calculations, usually performed at the ab initio level. Ideally, the method should contain only a small set of parameters that can be calculated from ab initio methods in a direct and easy manner. One such method is the density functional-based tight binding method with self-consistent charges[34] (SCC-DFTB), a method which has been successfully used in several applications, including biological systems[35] and inorganic materials.[36] SCC-DFTB allows for explicit calculations of the electronic structure (e.g., orbital population analysis and band structure calculations) of a system, which makes the approach suitable for studying, for example, charge transfer, a key ingredient in many chemical reactions.

There is an SCC-DFTB parameter set for the interactions between Zn, O, and H, called znorg-0-1.[37] It is based on the mio parameter set[34] for the O–O, O–H, and H–H interactions, while the Zn–O interaction was optimized to describe four-coordinated Zn and O atoms in the zincblende structure. There is also another mio-based Zn–O SCC-DFTB parameter set which was optimized for four-coordinated zincblende.[38] This parameter set is not publically available.

The znorg-0-1 potential has been used to study intrinsic defects in ZnO nanowires[39] and adsorption of various molecules on ZnO surfaces.[11,40,41] However, we have found that this potential displays overbinding at long Zn–O distances which tend to be present in structures with high coordination numbers. This has implications for the comparison of different structural phases of ZnO. For example, experimentally, the cubic NaCl-type structure of ZnO is less stable than the wurtzite structure, i.e., at ambient conditions the wurtzite structure is the most stable. However, in this paper we find that the znorg-0-1 potential gives the opposite result. The overbinding at long Zn–O distances is also problematic when comparing adsorption configurations for water molecules on ZnO surfaces in which the water oxygen atom coordinates one or two Zn atoms on the surface.

In an effort to overcome these problems, we have here generated a new SCC-DFTB repulsive potential which accurately handles all relevant phases of ZnO, for different bond lengths and coordination numbers. We achieved this by, in contrast to previous works, including structures of different coordination numbers directly in the parametrization. In this procedure, we have adopted ideas that has previously been used for the parametrization of a reactive force field for zinc oxide.[42] We will show that this new parametrization also improves the chemical description of ZnO surfaces, nanowires, and especially of the ZnO/water interface, for which we obtain good agreement with higher-level theoretical calculations (DFT) and to available experimental data.

In this work, we propose a method for efficient optimization of the repulsive potential within the SCC-DFTB framework. Specifically, based on a set of reference DFT calculations, we optimize the Zn–O repulsive potential to accurately describe the low-energy polymorphs of ZnO under various strains. We then use this optimized parameter set, which we will call znopt, to study ZnO nanowires and the ZnO/water interface. In the paper, we repeatedly demonstrate the accuracy of the calculated SCC-DFTB results with respect to calculations at the GGA-PBE[43] density functional level, which constitutes our framework of reference.

Methods

SCC-DFTB

The SCC-DFTB method has been described in detail elsewhere.[44,45] Here we will just briefly present the points necessary to appreciate the results of this paper.

SCC-DFTB is an approximative quantum-chemical method where the total energy is expressed as a second-order expansion of the DFT energy with respect to charge density fluctuations. The total energy is usually expressed as the sum of three terms: the band-structure term EBS, the second-order term Esecond, and the repulsive term Erep. The band-structure and repulsive terms are also found in “normal” (non-SCC) DFTB, while the second-order term is unique to the SCC-DFTB method.

The band-structure term corresponds to the sum of the energies of all occupied electronic eigenstates of the Hamiltonian. The eigenstates are expressed in a minimal basis of pseudoatomic orbitals generated from an all-electron calculation for each atomic species. The (distance-dependent) overlap matrix elements and the matrix elements representing the Hamiltonian in this basis are precalculated once and stored in Slater–Koster tables.

The second-order term contains the energy contributions due to charge density fluctuations in the system. These fluctuations are approximated by atomic charges derived from Mulliken population analysis.[46] The energy associated with the interaction of these atomic Mulliken charges is derived using the approximation of spherically symmetric Gaussian-shaped atomic charge densities, leading to an analytical expression that determines the on-site and interatomic contributions to the second-order term. These contributions depend on the individual atomic hardnesses, which are usually expressed in terms of the Hubbard U parameter.

The band-structure term and the second-order term constitute the electronic part of the total energy. The total energy of the system also includes the repulsive term, which contains the ion–ion repulsion (hence the name) as well as exchange-correlation contributions and “double-counting” corrections. In practice, it is usually approximated by a sum of pairwise atomic interactions:where Vrep(r) is the potential between atoms I and J at a distance r. These pairwise potentials are obtained by fitting them either to experimental data, or to the results of theoretical calculations of higher accuracy, e.g., ab initio calculations.

Scheme for Fitting of the Repulsive Potential

To make our newly generated repulsive potential compatible with the well-established mio parameter set,[34] and to some extent with the znorg-0-1(37) parameter set (which also involves interactions between Zn and the elements C, S, N, and P, that are not treated in the present work), we have kept all the Hubbard parameters for the elements and Slater–Koster tables from the mio and znorg-0-1 parameter sets. The only part that we have changed is the Zn–O repulsive potential from the znorg-0-1 set (i.e., all the other repulsive potentials are also kept the same as in znorg-0-1).

The znorg-0-1 repulsive potential was only trained to reproduce results for the cubic zincblende polymorph of ZnO.[37] In the present work we show that the repulsive energy term obtained with the znorg-0-1 potential is too small at longer distances (above 2 Å) which are present in ZnO structures with high Zn coordination numbers. Thus, using the znorg-0-1 potential, the NaCl-type polymorph of ZnO, which is six-coordinated and therefore exhibits relatively long Zn–O bonds, is considerably more stable than the experimentally found wurtzite structure, which is four-coordinated and exhibits relatively short Zn–O bonds. We therefore set out to generate a new repulsive potential (znopt), with an increased repulsion at long distances, to obtain the correct stability order of the wurtzite and NaCl-type polymorphs of ZnO.

To achieve this, we performed volume scans for ZnO in the wurtzite structure and the NaCl-type structure (see Figure 1) using density functional theory (details are given in section 2.3). The necessity of including structures with different coordination numbers and different bond lengths in the parametrization has been pointed out before,[45,47] but this seems to not yet have been the practice for SCC-DFTB parametrizations of solid-state materials.

Figure 1

Unit cells for the considered bulk structures of ZnO: wurtzite, NaCl-type, zincblende, CsCl-type, graphitic, body-centered tetragonal (BCT), and cubane. Oxygen is depicted in red and zinc in gray. The coordination numbers for the Zn and O atoms are given in parentheses.

Atomic positions and lattice parameters were optimized with DFT, giving a Zn–O bond length of 2.17 Å in the NaCl-type structure and a (shortest) bond length of 2.01 Å in the wurtzite structure. The optimized DFT geometry was then scaled isotropically in a volume range of 73% to 133% in steps of approximately 6%. Single-point energy calculations were performed at each volume, using both DFT and SCC-DFTB. The ideal SCC-DFTB repulsive energy, i.e., the repulsive energy which would result in perfect agreement between DFT and SCC-DFTB and which we denote Ẽrep, is the difference between the DFT total energy and the electronic part of the SCC-DFTB energy for each structure S, calculated with respect to the DFT-optimized wurtzite structure Sref:where EDFT is the total DFT energy and Eel is the electronic part of the SCC-DFTB energy. The ideal value for Ẽrep(Sref) is the repulsive energy value that gives the desired atomization energy of the solid. It is, however, in principle possible to choose any arbitrary value for this quantity. The absolute values of the atomization energies in the SCC-DFTB calculations will then not match the DFT values, although the relative atomization energies will match. By sacrificing some of the accuracy of the atomization energy, it may be possible to improve other aspects such as surface relaxation[37] or water adsorption energies. This has been the approach used in this work, i.e., the absolute values of the atomization energies were varied until all other properties considered in this work were reproduced satisfactorily by the SCC-DFTB calculations. The error introduced in the absolute atomization energies can thus effectively be regarded as an error in the description of the isolated atoms, but because we are primarily interested in the chemistry of ZnO with a coordination number near 4, we believe this to be a sound approach.

The goal of the parametrization was thus to find Vrep of eq 1 so that the difference between Ẽrep and Erep was minimized. We chose to express Vrep using the four-range Buckingham potential, which divides Vrep into four different functional forms depending on the interatomic distance:This expression allows for some flexibility without giving rise to unphysical variations or kinks in the potential, as the first and second derivatives are continuous in the entire distance range. In the parametrization, we varied the five limits rmin to rmax manually and added a tapering function to let the repulsive potential smoothly approach zero at r = 3.0 Å, which is smaller than the second nearest neighbor distance in any of the considered polymorphs of ZnO. We subsequently optimized the coefficients to minimize the difference between Ẽrep and Erep. The fitting was done with the GULP[48] software.

Electronic Structure Methods

The SCC-DFTB calculations were performed with the DFTB+ package[49] together with both the original znorg-0-1 repulsive potential by Moreira et al.[37] and with our new improved repulsive potential znopt. An empirical dispersion correction, similar to the Grimme dispersion correction[50] for DFT, was not used, because we validated our SCC-DFTB repulsive potential against results obtained with DFT without such a dispersion correction.

The DFT reference calculations were performed using the generalized gradient approximation exchange-correlation functional PBE[43] in an implementation involving a plane-wave basis set with energy cutoff 500 eV and projector augmented wave (PAW) type pseudopotentials.[51,52] We explicitly treated one, six, and twelve valence electrons for H, O, and Zn, respectively. The PBE DFT functional was used as a reference because it was the functional used to derive the SCC-DFTB atomic basis for Zn.[37] The PBE calculations were performed using the VASP[53−55] package.

Periodic boundary conditions in three dimensions were employed throughout. Geometry optimizations were performed using the conjugate gradient algorithm until all forces on the atoms were smaller than 5 meV/Å. Convergence tests with respect to the k-point sampling were made for each system under study (for the bulk calculations, we typically used a Γ-centered 7 × 7 × 7 k-point grid and correspondingly smaller grids for the surface and nanowire calculations).

Systems Studied

ZnO Bulk Polymorphs

In this work, we have calculated properties such as lattice parameters, bulk moduli, and band gaps for the CsCl-type, NaCl-type, wurtzite, zincblende, graphitic, BCT, and cubane structures of ZnO. The unit cells for each of the seven polymorphs are shown in Figure 1. The hexagonal wurtzite structure is characterized not only by its lattice parameters a and c but also by the internal parameter u which corresponds to the fractional displacements of the Zn and O sublattices along the c direction. We have chosen to model the BCT structure as an “ideal” BCT structure, i.e., with a = b so that the structure is tetragonal. The BCT structure then possesses two internal parameters, one of which we call u and which corresponds to the fractional displacements of the Zn and O sublattices along the a or b direction, and the other which corresponds to the Zn–O–Zn angle α depicted in Figure 1. The internal parameters u and α are defined in a similar way for the cubane polymorph.

In the parametrization procedure, only the wurtzite and NaCl-type polymorphs were included, while the other polymorphs were used for testing purposes. Thus, to ensure the transferability of our generated parameter set, we validated the parameters by calculating energy–volume scans for the high-density CsCl, intermediate-density graphitic and zincblende, and low-density BCT and cubane polymorphs. The optimized structures were strained from 73% to 133% in steps of approximately 6%. For each volume under consideration, the c/a ratio in the wurtzite, BCT, and graphitic structures was optimized, and all atoms were allowed to fully relax.

The bulk modulus B0 was obtained by fitting the energy–volume curves to a Murnaghan type equation of state. The atomization energies Eat of the bulk polymorphs were calculated aswhere n is the number of ZnO formula units in the crystal unit cell, and EZn-atom, EO-atom, and EZnO-crystal correspond to the energies of an isolated Zn atom, an isolated (spin-polarized) O atom, and the ZnO crystal, respectively.

Clean ZnO Surfaces

In the calculation of structural properties of the clean ZnO surfaces, we used twenty-layer thick ZnO(101̅0) and ZnO(112̅0) slabs with a vacuum gap of at least 15 Å between neighboring systems. The systems were constructed using the bulk lattice parameters native to each method (i.e., the PBE-optimized bulk wurtzite lattice was used to construct the surfaces for the PBE calculations, etc.). Side and top views of these systems are shown in Figure 2. We characterized the surface structure through the angle θ that the ZnO “dimers” in the top layer along the polar ZnO[0001]/ZnO[0001̅] directions make to the corresponding layers in the bulk (indicated in Figure 2). Additionally, we calculated the surface energy aswhere Eslab is the total energy of the slab supercell, n is the number of ZnO formula units in the cell, EZnO-bulk is the energy of a formula unit of ZnO in the wurtzite bulk, and A is the total surface area of the slab supercell (twice the area of one face). The surface energies were converged within 0.01 J/m2 for the eight-layer thick slabs. For this reason, eight-layer thick slabs were used for the adsorption of water (see the following sections).

Figure 2

Top and side views of the wurtzite ZnO(101̅0) and ZnO(112̅0) structures geometry-optimized using the SCC-DFTB parametrization znopt. In the top views, two layers are shown with the top layer atoms larger and brighter than the second-layer atoms. In the side views, six layers are shown for ZnO(101̅0) and four layers for ZnO(112̅0), and the relaxation angle θ is indicated.

No calculations were performed for the polar ZnO(0001) and ZnO(0001̅) surfaces, because these typically exhibit significant surface reconstructions and possibly surface metallization.[56,57] Because of the complexity involved in modeling these surfaces, they lie outside the scope of the current article.

ZnO Nanowires

We generated wires in the hexagonal wurtzite structure and in the NaCl-type structure, see Figure 3, and calculated formation energies as a function of wire diameter d for diameters in the range 9 to 44 Å using SCC-DFTB, and for diameters in the range 9 to 27 Å using PBE. The largest nanowire studied using PBE contained 200 formula units, while the largest nanowire studied using SCC-DFTB contained 588 formula units. A vacuum gap of at least 15 Å was introduced between neighboring wires. The hexagonal wires of the wurtzite structure were terminated by {101̅0} facets and were periodic in the [0001̅] direction. The square wires of NaCl-type structure were terminated by {100} facets and were periodic in the [001] direction. The formation energies of the nanowires per ZnO were calculated aswhere n is the number of ZnO formula units in the supercell, Ewire is the total energy of the wire supercell, and EZnO-bulk is the energy per ZnO formula unit in wurtzite bulk ZnO.

Figure 3

View along the [0001] direction of a wurtzite ZnO nanowire and along the [001] direction of a NaCl-type ZnO nanowire. The diameter d is indicated by the black arrow.

Water Adsorption on ZnO Surfaces

Single water molecules were adsorbed onto the ZnO(101̅0) and ZnO(112̅0) surfaces in various adsorption configurations, similar to the calculations performed by Meyer et al.[7] and grosse Holthaus et al.[11] The calculations were performed in a 3 × 2 supercell for ZnO(101̅0) (ca. 10 Å × 10 Å) and in a 2 × 2 supercell for ZnO(112̅0) (ca. 11 Å × 10 Å). Additionally, full-monolayer calculations in a 2 × 2 cell for both surfaces were performed in fully dissociated, half-dissociated, and molecular configurations. The water molecules were adsorbed onto only one side of the eight-layer thick slab. The adsorption energy per water molecule was calculated aswhere N is the number of water molecules per supercell, Etot is the energy per supercell of the system with water adsorbed onto the ZnO slab, Eslab is the energy per supercell of the clean optimized ZnO slab, and EH is the energy of an optimized isolated water molecule. Thus, the more negative the adsorption energy, the more stable the system.

Results and Discussion

The Optimization of the SCC-DFTB Repulsive Potential

Our optimized znopt repulsive potential is compared to the original znorg-0-1 repulsive potential in Figure 4. At distances longer than 2 Å, the znopt potential is more repulsive than the znorg-0-1 potential. Figure 5 shows energy–volume curves calculated using znorg-0-1, znopt, and the reference PBE method for the wurtzite and NaCl-type bulk structures of ZnO. These were the polymorphs that were explicitly fitted against in the znopt parametrization procedure. By inspection of Figure 5, it is clear that we have addressed the problem of the too stable NaCl-type structure that plagued the znorg-0-1 repulsive potential. In the following sections, we will show that we improve the descriptions of many other types of systems as well.

Figure 4

The original znorg-0-1 repulsive potential[37] (dashed line) and the new optimized repulsive potential znopt of this work (solid line) as a function of Zn–O distance r.

Figure 5

Volume scans calculated with three different methods for the bulk polymorphs of ZnO included in the parametrization of the znopt potential (the NaCl-type and wurtzite polymorphs): PBE, the znopt potential, and the original znorg-0-1 potential.[37].

ZnO Bulk Polymorphs

The calculated optimized lattice parameters, band gaps, and bulk moduli of the various bulk phases considered (CsCl-type, NaCl-type, wurtzite, zincblende, graphitic, BCT, and cubane) of ZnO are given in Table 1. The energy–volume curves calculated with both PBE and the znopt SCC-DFTB repulsive potential for all considered phases are shown in Figures 5 and 6. We find that the agreement between the znopt and the reference PBE results in general is very good, both for those structures that formed part of the parametrization and for the others. With both methods, the most stable structure is the wurtzite structure (as is the case in experiment), followed closely by the zincblende structure. At high pressures, a phase transformation into the NaCl-type structure can be expected while at highly negative pressures, the BCT and cubane structures are the most stable.

Table 1

Calculated Lattice Parameters a and c, Internal Parameter u, Atomization Energy Eat, Bulk Modulus B0, and Band Gap Egap for Various ZnO Polymorphs at Different Levels of Theorya

	PBE	znopt	znorg-0–1
Wurtzite (P63mc, no. 186)
a [Å]	3.29	3.21	3.29
c [Å]	5.31	5.25	5.38
u	0.379	0.374	0.375
Eat [eV]	7.37	9.88	9.92
B0 [GPa]	129	161	161
Egap [eV]	0.73	4.33	3.78
NaCl (Fm3̅m, no. 225)
a [Å]	4.34	4.34	4.36
Eat [eV]	7.07	9.56	10.06
B0 [GPa]	166	227	302
Egap* [eV]	0.77	2.32	2.23
Zincblende (F4̅3m, no. 216)
a [Å]	4.62	4.54	4.65
Eat [eV]	7.35	9.88	9.92
B0 [GPa]	129	162	160
Egap [eV]	0.65	4.27	3.73
Graphitic (P63/mmc, no. 194)
a [Å]	3.46	3.33	3.53
c [Å]	4.58	4.87	4.49
Eat [eV]	7.23	9.75	9.87
B0 [GPa]	99	178	195
Egap [eV]	0.95	4.20	3.35
CsCl (Pm3̅m, no. 221)
a [Å]	2.69	2.71	2.64
Eat [eV]	5.93	8.86	9.34
B0 [GPa]	160	264	385
Egap* [eV]	0	0.48	1.19
(Ideal) BCT (P42/mnm, no. 136)
a [Å]	5.63	5.51	5.64
c [Å]	3.29	3.25	3.33
u	0.363	0.367	0.365
α [deg]	90.1	91.9	92.0
Eat [eV]	7.32	9.73	9.77
B0 [GPa]	105	142	131
Egap [eV]	0.75	4.52	3.94
Cubane (I4̅3m, no. 217)
a [Å]	6.29	6.19	6.34
u	0.326	0.329	0.327
α [deg.]	91.0	91.9	91.7
Eat [eV]	7.15	9.36	9.41
B0 [GPa]	97	111	104
Egap [eV]	1.32	5.33	4.61

The band gaps for the NaCl and CsCl type structures are indirect (denoted with an asterisk). The angle α for the BCT and cubane structures is defined in Figure 1. The space groups and space group numbers for the various structures are also given.

Figure 6

Volume scans for five different bulk polymorphs of ZnO not included in the znopt parametrization (CsCl-type, graphitic, zincblende, BCT, and cubane), calculated with and PBE and znopt.

There are, however, some minor discrepancies between the PBE and znopt results. For example, the equilibrium volume of the wurtzite structure is somewhat smaller with znopt compared to PBE, although the znopt-optimized lattice parameters are closer to the experimental values (a = 3.250 Å, c = 5.207 Å, u = 0.3825).[58] Additionally, the band gap obtained with znopt is 4.33 eV while the PBE band gap is 0.73 eV, so the znopt band gap is closer to the experimental value of 3.4 eV.[3] PBE greatly underestimates the band gap because of electron self-interaction, while SCC-DFTB overestimates the band gap because of the use of a minimal basis set. Here, we should point out that the original znorg-0-1 potential in fact gives a somewhat smaller band gap (3.78 eV). This is a result of the fact that the Zn–O bonds are longer in the znorg-0-1-optimized structure compared to the znopt-optimized structure.

The underestimation of the band gap by PBE reaches an extreme for the CsCl-type polymorph of ZnO, where a metallic solution is obtained (as also noted by Wróbel and Piechota).[18] With znopt, a small indirect band gap of 0.48 eV remains. Uddin and Scuseria[17] performed hybrid density functional calculations on this polymorph and also noted a decrease of the band gap compared to the wurtzite structure, albeit only by 1.2 eV (from 2.9 to 1.7 eV). For the NaCl-type structure, we obtain an indirect band gap of 2.32 eV with znopt, which is similar to the experimental value of 2.45 ± 0.15 eV.[59]

A more detailed view of the electronic structure of the various polymorphs is given in Figure 7. In Figure 7a, the electronic DOS obtained by znopt and PBE is shown for the wurtzite polymorph of ZnO. The energies have been aligned with respect to the O2s energy, which is the most “corelike” state that is explicitly treated in the calculations. It is clear that the larger band gap obtained with znopt mainly comes from a shift of the unoccupied states, i.e., the conduction band. The occupied valence band, on the other hand, is well-described compared to PBE. In Figure 7b, schematic views of the DOS for all seven polymorphs are shown. The agreement between the znopt and PBE-calculated valence bands for the four-coordinated wurtzite, zincblende, BCT, and cubane-type structures, as well as the five-coordinated graphitic structure, is very good. The positions of the valence bands relative to the O2s states with PBE and znopt match very well, and the znopt-calculated valence bands are only slightly more narrow than the corresponding PBE-calculated valence bands. For the higher-coordinated NaCl-type and CsCl-type polymorphs, the agreement between znopt and PBE is worse, with the valence bands being much too narrow in the znopt calculations. Thus, the valence bands of these high-coordinated phases are not very well described, and the znopt potential is consequently not very suitable to study these phases, even if the formation energies relative to the wurtzite phase are qualitatively correct.

Figure 7

(a) Electronic DOS for ZnO in the wurtzite structure, calculated with the znopt potential and PBE. The energies are aligned with respect to the average O2s energy. A Gaussian broadening of width 0.2 eV has been applied. In the znopt-panel, the major projections onto the atomic orbitals are given. The (unoccupied) conduction band is lightly shaded. (b) Schematic electronic DOS for seven different polymorphs of ZnO. The valence band (VB) is heavily shaded and the conduction band (CB) is lightly shaded. In each panel, the upper band structure was obtained with znopt, and the lower band structure was obtained with PBE.

The znopt-calculated bulk moduli are in general too high compared to PBE. This can be seen in Figures 5 and 6, where the znopt-calculated energy–volume curves are too steep compared to PBE, although good agreement with experiment is obtained (B0 = 142.6 GPa for wurtzite and B0 = 202.5 GPa for the NaCl-type structure).[60]

We note that the atomization energies of all the polymorphs are quite overestimated, as is often the case in SCC-DFTB calculations.[37,61] The reason for this is that by imposing the correct atomization energy (which we found can be achieved by increasing the Zn–O repulsive potential by roughly 0.65 eV near the Zn–O equilibrium distance in wurtzite, cf., eq 2), the adsorption of water molecules (see coming sections) on the ZnO surfaces becomes far too weak. Our reported znopt potential thus corresponds to a trade-off between correct bulk ZnO and water adsorption descriptions.

Finally, we comment on the PBE-calculated results for the BCT polymorph. Our calculated equilibrium volume (26.1 Å3/ZnO) and relative stability (0.05 eV/ZnO less stable than the wurtzite structure) agree very well with results of Morgan[29] and Zwijnenburg et al.[30] However, the results are considerably different from the equilibrium volume (33 Å3/ZnO) and relative stability (1.15 eV/ZnO less stable than wurtzite) reported by Zhang et al.,[33] who performed similar PBE calculations and suggested that the cubane polymorph is the most stable low-density polymorph of ZnO. In contrast to the results of Zhang et al.,[33] our calculations show that the BCT structure in fact is more stable than the cubane structure.

ZnO Surfaces

The calculated surface energies of the ZnO(101̅0) and ZnO(112̅0) together with the corresponding relaxation angles are given in Table 2. Although our znopt-calculated surface energies are still too high compared to PBE, we improve on both the surface energies and relaxation angles compared to znorg-0-1. SCC-DFTB gives higher surface energies than PBE because of the higher atomization energy (see Table 1), so that it is more unfavorable to break ZnO bonds.

Table 2

Calculated Surface Energies Esurf and Relaxation Angles θ (defined in Figure 2) for the ZnO(101̅0) and ZnO(112̅0) Surfaces of the ZnO Wurtzite Structure

	PBE	znopt	znorg-0-1
Esurf(101̅0) [J/m2]	0.86	1.22	1.31
Esurf(112̅0) [J/m2]	0.89	1.29	1.38
θ(101̅0) [deg]	10.4	11.1	13.1
θ(112̅0) [deg]	11.9	11.9	14.1

ZnO Nanowires

The calculated formation energies as a function of the nanowire diameter are given in Figure 8 for the reference PBE method, as well as the znopt and znorg-0-1 potentials. The PBE and znopt results are comparable, with the NaCl-type nanowire being consistently less favorable than the wurtzite-type nanowire (experimentally, ZnO nanowires indeed exhibit the wurtzite structure).[62] In fact, the difference in formation energies between the two different phases appear to be almost constant over the diameter range considered, with the difference for PBE amounting to 0.24 eV and the difference for znopt amounting to 0.35 eV, i.e., slightly larger than the PBE difference.

Figure 8

Formation energies ΔEf for ZnO nanowires in the wurtzite and NaCl-type structures as a function of wire diameter, calculated with PBE, the znopt potential, and the znorg-0-1 potential.

Conversely, the znorg-0-1 potential favors the NaCl-type nanowire for all wire diameters larger than roughly 10 Å. The results for the nanowires thus follow closely the results of the bulk calculations, where the NaCl-type structure was more stable than the wurtzite structure for the znorg-0-1 potential.

Both the formation energies for the NaCl-type and wurtzite-type nanowires are higher with znopt than with PBE (for example, the formation energy of the wurtzite nanowire of diameter 23 Å is 0.22 eV/ZnO with PBE but 0.31 eV/ZnO with znopt). This is consistent with the higher surface energies obtained with znopt (see previous section).

We additionally calculated band gaps of the nanowires as a function of wire diameter. In agreement with a previous SCC-DFTB study,[37] we find that the band gap decreases with increasing size until it converges to a value slightly below the one for the optimized bulk. This is due to surface states in the nanowire, and similar effects are observed for the surface slabs.

The ZnO/Water Interface

The structures of all studied water/ZnO systems are shown in Figure 9, and the corresponding adsorption energies are given in Table 3 for PBE, znopt, and znorg-0-1. The adsorption energies are also given in a graphical diagram in Figure 10. In general, the znorg-0-1 potential overestimates the interaction between the water molecules and the ZnO surface (resulting in much too negative, i.e., too stable, adsorption energies), compared to both PBE and znopt. In ref (11), even more negative adsorption energies for the ZnO(112̅0)/water interface obtained with the znorg-0-1 potential were reported, but those numbers contained errors.[63]

Figure 9

Top views of various water arrangements on the ZnO(101̅0) and ZnO(112̅0) surfaces, as calculated with the znopt potential, are shown. The coloring scheme follows that of Figure 2, with the addition of blue water oxygen atoms and white hydrogen atoms.

Table 3

Calculated Adsorption Energies (in eV per water molecule) for Water on the ZnO(101̅0) and ZnO(112̅0) Surfacesa

configuration	PBE	znopt	znorg-0-1
Single Water Molecule on ZnO(101̅0)
M-1	–0.98	–1.06	–1.26
M-2	–0.60	–0.69	–0.88
D-1	–0.89	–0.90	–1.31
D-2	–0.66	–0.62	–0.72
Full Water Monolayer on ZnO(101̅0)
molecular	–1.07	–1.09	–1.27
half-dissociated	–1.17	–1.14	–1.29
dissociated	–0.96	–0.89	–0.99
Single Water Molecule on ZnO(112̅0)
M-1	–1.07	–1.21	–1.39
M-2	–0.92	–1.02	–1.21
D-1	+0.34	+0.63	+0.60
Full Water Monolayer on ZnO(112̅0)
molecular-1	–0.93	–1.17	–1.34
molecular-2	–0.92	–1.08	–1.25
half-dissociated	–1.06	–1.04	–1.15
dissociated	–0.99	–0.87	–0.95

The structures are displayed in Figure 9.

Figure 10

Adsorption energies of all the ZnO/water structures considered in this work, except the D-1 configuration for ZnO(112̅0), which has been removed for clarity. The structures are shown in Figure 9, and the values can also be found in Table 3.

The agreement between PBE and the newly generated znopt potential is in general very good. One of the reasons our znopt potential gives weaker adsorption than the znorg-0-1 potential is the greater value of the Zn–O repulsive potential at the water-Zn distances (2.0–2.1 Å). In the following, we will discuss in detail the results obtained for single-molecule and full-monolayer adsorption of water on the ZnO(101̅0) and ZnO(112̅0) surfaces.

Single Water Molecule Adsorption on ZnO(101̅0)

Meyer et al.[7] also performed DFT calculations of isolated water molecules on ZnO(101̅0) using PBE. We found that only four of the nine adsorption structures reported by Meyer et al.[7] correspond to true local minima (at least, using our implementation of DFT, the other five structures represented shallow regions on the potential energy surface). The four structures are the d, e, f, and i structures in Figure 1 of ref (7), and we choose to call them M-2, D-2, D-1, and M-1, respectively (see Figure 9), where the M denotes molecular adsorption and the D denotes dissociated adsorption, and the species labeled “1” is more stable than the species labeled “2”.

In the molecular M-1 configuration (Figure 9a), the water oxygen atom coordinates one Zn atom, which adopts an almost bulklike tetrahedral coordination of four oxygen atoms. One of the water hydrogen atoms forms a hydrogen bond to an oxygen atom on the surface. In the dissociated D-1 configuration (Figure 9c), the water (hydroxyl) oxygen is in a bridgelike position between two neighboring Zn atoms, and the dissociated H atom coordinates one of the surface oxygens. The M-2 configuration (Figure 9b) corresponds to molecular adsorption almost on top of a Zn atom, and the D-2 configuration (Figure 9d) is the corresponding dissociated configuration.

For these four configurations at the PBE level, we obtain the same relative stability and similar adsorption energies as those obtained by Meyer et al.[7] The M-1 configuration is the most stable with an adsorption energy of −0.98 eV, followed by the D-1 configuration (Eads = −0.89 eV). Considerably less stable are the M-2 (Eads = −0.60 eV) and D-2 (Eads = −0.66 eV) configurations.

The old znorg-0-1 potential gives much too strong adsorption for these four adsorption structures, particularly for the most stable M-1 and D-1 configurations where the adsorption energies are 0.3–0.4 eV more negative than the PBE values. In contrast, our optimized znopt potential is in excellent agreement with the PBE results, and also predicts the M-1 configuration to be more stable than the D-1 configuration, as is the case in PBE but not with znorg-0-1. The improvement is displayed in a clear manner in Figure 10a. With znorg-0-1, the D-1 configuration becomes the most stable because there is an increased coordination of the water O atom to the ZnO surface, so the relative stability of the M-1 (low-coordinated) and D-1 (high-coordinated) adsorption configurations of water on ZnO(101̅0) exactly follow the relative stability trend for the bulk phases of wurtzite (low-coordinated) and NaCl-type (high-coordinated), irrespective of the method used. This could explain why the correct stability order of M-1 and D-1 is obtained with znopt but not with znorg-0-1.

Full Water Monolayer Adsorption on ZnO(101̅0)

Increasing the coverage of water to full monolayer coverage (1 ML, one water molecule per surface Zn atom) on the ZnO(101̅0) surface, we distinguish between three different cases: 1 × 1 molecular adsorption (Figure 9e), 2 × 1 half-dissociated adsorption (Figure 9f), and 1 × 1 dissociated adsorption (Figure 9g). In the half-dissociated case, every other water molecule along the [12̅10] direction of ZnO is dissociated. In the PBE calculations, the half-dissociated adsorption is the most stable (Eads = −1.17 eV), in agreement with experiment[8] and previous DFT calculations.[7,64] The molecular configuration is less stable (Eads = −1.07 eV) and the dissociated configuration even more so (Eads = −0.96 eV).

Excellent agreement between the PBE and znopt calculations is obtained also for the full coverage case of water adsorption on ZnO(101̅0) (see Table 3 and Figure 10b). For the original znorg-0-1 potential, the most favorable adsorption energy in the full coverage case (the half-dissociated configuration) is −1.29 eV per water molecule, which is less stable than the isolated water molecule in the D-1 (Eads = −1.31 eV) configuration, i.e., there is an effective repulsion between the water molecules. However, it is known both from experiment[8] and higher-level theoretical calculations[7,64] that the 2 × 1 half-dissociated network stabilizes the adsorbed water. The performance of the znorg-0-1 potential is thus not very good in this case. The znopt potential performs much better, with adsorption energies that are within 0.07 eV of the PBE-calculated values.

Single Water Molecule Adsorption on ZnO(112̅0)

On the ZnO(112̅0) surface, we find two stable adsorption sites that are both of the molecular kind (M-1 and M-2). These correspond to configurations B and A in ref (11), respectively. In the M-1 configuration (Figure 9h), the water molecule forms two hydrogen bonds to the surface, while in the M-2 configuration (Figure 9i), only one hydrogen bond is formed. Here, the znopt potential overestimates the binding strength by about 0.1 eV compared to the PBE values (Table 3 and Figure 10c). These values considerably improve upon the values obtained with znorg-0-1, that are about 0.3 eV too negative compared to PBE. Nevertheless, the relative stabilities of the two adsorption sites are the same with all three methods, and the differences in adsorption energies between M-1 and M-2 are all about 0.2 eV.

For completeness, we have calculated the adsorption energy of a dissociated water molecule on this surface (D-1, see Figure 9j). In this case, we were unable to stabilize the dissociated hydrogen atom near the remaining hydroxyl fragment. The only local minimum found had the dissociated hydrogen atom being displaced almost a full surface unit cell along the [0001̅] direction away from the hydroxyl group. Although being a local minimum, the adsorption energy obtained with PBE, znopt, and znorg-0-1 is positive, which implies that the structure is not stable with respect to desorption. The adsorption energies obtained with either SCC-DFTB repulsive potential is about 0.3 eV higher than the PBE value, but because the structure is so unstable, it is unlikely to appear during for example molecular dynamics simulations.

Full Water Monolayer Adsorption on ZnO(112̅0)

For the ZnO(112̅0) surface, similar to the case of the ZnO(101̅0) surface, the adsorbed water monolayer can form molecular (Figure 9k and 9l), half-dissociated (Figure 9m), or dissociated (Figure 9n) structures. PBE predicts that the half-dissociated configuration is the most stable. This result is consistent with the DFT calculations of große Holthaus et al.[11] but at odds with the results of Cooke et al.,[10] who found that the fully dissociated configuration is 0.06 eV per water more stable than the half-dissociated configuration. This may be related to the choice of supercell, as Cooke et al.[10] used a 1 × 1 supercell, while in the current work (as well as in ref (11)) a 2 × 2 supercell was used, which gives room for two nonidentical molecularly adsorbed water molecules alongside the dissociated water molecules. In any case, the energy differences between the different configurations are small, and coexistence of various configurations under normal conditions can be expected, as noted previously.[10,11]

Table 3 and Figure 10d show that the agreement between the PBE and SCC-DFTB calculations is not very good in that the PBE result favors half-dissociated adsorption while the SCC-DFTB results strongly favor molecular adsorption. Even so, the znopt potential significantly improves on the znorg-0-1 results, where the largest error compared to PBE is 0.4 eV for the molecular-1 configuration, while our largest error with the znopt potential is 0.2 eV (also for the molecular-1 configuration). However, for a study where the ZnO(112̅0)/water interface is in focus (not this one), the SCC-DFTB parameters may need to be improved. This may involve changing the O–H repulsive potential or the Slater–Koster tables. In the present work, we limited ourselves to only changing the Zn–O repulsive potential, and already this modification resulted in major improvements, as we have seen.

große Holthaus et al.,[11] who used the original znorg-0-1 potential, stated that the half-dissociated configuration was the most stable one at full monolayer coverage, but they did not actually report the adsorption energy of a fully molecular adsorption layer. große Holthaus et al.[11] subsequently performed MD simulations with additional layers of liquid water on top of the surface and found that the dissociated or half-dissociated water layers nearest to the ZnO surface began to convert into molecular layers. The authors attributed this effect to the formation of a developed hydrogen-bonded network above the surface, although a more likely explanation, in our opinion, is that the molecular monolayer actually is the most stable for both the znorg-0-1 and the znopt potentials, as we have shown here.

Conclusions

We have developed an efficient scheme to generate SCC-DFTB repulsive potentials. The scheme is based on energy–volume scans of bulk polymorphs with different coordination numbers and was applied here to ZnO, although it can easily be generalized to other materials. The key to a successful and transferable parametrization appears to be the inclusion of structures with different coordination numbers in the parametrization procedure, something which to our knowledge has not been done before for SCC-DFTB parametrizations of solid-state materials.

The previously reported znorg-0-1(37) potential gives the incorrect stability order of the wurtzite and NaCl phases (Figure 5). Our newly generated znopt potential, on the other hand, gives the correct stability order for the NaCl-type and wurtzite structures and additionally performs very well for other possible bulk structures of ZnO (CsCl-type, zincblende, graphitic, body-centered tetragonal (BCT), and cubane). The improved description of chemical properties extends into lower-dimensional properties such as surface relaxation, surface energies, and nanowire formation energies. Finally, we found that the znopt potential greatly improves on the description of the ZnO/water interface at both low and high coverage, particularly for the ZnO(101̅0) surface where excellent agreement is obtained compared to reference DFT data.