Making amorphous ZnO: Theoretical predictions of its structure and stability

ZnO is a transparent semiconductor with optoelectronic, thermoelectric, and sensor applications, where using amorphous thin ﬁlms presents great advantages. However, growing amorphous (a) ﬁlms of pure ZnO proved challenging due to their rapid crystallization. We investigated the ability of bulk ZnO to form glass structures using well-tested interatomic potentials and a melt and quench procedure within isothermal-isobaric (NPT) ensemble. The geometries of some of the resulting structures were further optimized using density functional theory (DFT) calculations with the PBE functional. We demonstrate that cooling rates in melt and quench procedure equal or exceeding 100 Kps − 1 lead to formation of stable amorphous structures. However, ZnO samples tend to crystallize at lower cooling rates. This result does not depend on the size of the periodic cell used in the calculations for cells containing more than 324 atoms. Using simulation cells with up to 768000 atoms, we demonstrate that the expected average glass density is about 5.04 gcm − 3 and the coordination numbers of Zn and O atoms are around 3.9. We calculate radial distribution functions and characterize the structures of amorphous ZnO samples. Using both the activation-relaxation technique and simulated annealing, we show that the obtained amorphous structures have low propensity to crystallization.


I. INTRODUCTION
Crystalline zinc oxide (c-ZnO) is an important n-type wide band-gap oxide with a wide range of applications as a transparent conductor, varistor, gas sensor, piezoelectric transducer, and UV detector [1][2][3]. Recently, ZnO films have attracted a lot of attention due to their applications in, for example, the manufacture of heat mirrors used in gas stoves, as conducting coatings in aircraft glass for avoiding surface icing, and as thin film electrodes in amorphous silicon solar cells [4]. The attention has further been boosted by the development of efficient methods of film deposition using dip coating [4]. ZnO coatings are also produced using a range of other techniques, such as radio-frequency sputtering [5][6][7], pulse laser deposition [8], prefiring-final annealing method [9], and microemulsion method in AOT reverse micelles [10]. Thin ZnO films have been grown on different substrates, e.g., silica glass [8], soda-lime-silica glass [9], bare glass, and SiN/glass [7]. The thickness of these films range from 5 to 100 nm.
These, as well as recent advances in development of ZnO transparent thin film transistors, [7,11] have renewed interest in growing thin amorphous (a)-ZnO films. Several experimental studies have suggested that ZnO films grown on different substrates are indeed amorphous but did not provide compelling evidence to justify that. Where high-intensity x-ray diffraction or grazing incidence x-ray and neutron diffraction has been used, the films proved to be nanocrystalline [12,13]. The structure and electronic properties of polycrystalline thin a-ZnO films have, therefore, been studied in much greater * david.fonz.11@ucl.ac.uk detail and whether and how amorphous ZnO thin films can be stabilized still remains unclear.
It is often suggested that rapid cooling of the melt can freeze disorder even in a non-glass-former [14]. There is a question, however, whether the structures of the glass produced from the melt and those of noncrystalline thin films grown on a substrate using different deposition methods are indeed equivalent. Recent nuclear magnetic resonance (NMR) experiments on alumina and hafnia films grown using different techniques demonstrated that both the Al coordination populations and the extent of the "topological" disorder (distributions of bond distances and angles, etc. that contribute to NMR line widths and shapes) are similar in the glass and the film [15,16]. This led to speculations that the amorphous film reaches a supercooled liquidlike metastable equilibrium during deposition [15,17,18]. Recent simulations [19] suggest that the glassy materials prepared by physical vapor deposition (PVD) have very similar structures to those prepared by gradual cooling from the liquid phase, and that PVD glasses correspond to liquid-cooled glasses prepared at extremely slow cooling rate (CR). These observations provide support to using a melt/quench procedure and 3D bulk samples, which is a common practice for simulating glasses, to model the structure inside thin oxide films grown on substrates.
Here we use melt/quench procedure to investigate whether ZnO can produce a bulk glass state. These results will be useful for growing and understanding the structures of thin a-ZnO films. Since experimental structural information on a-ZnO is almost nonexistent, theoretical modeling can provide valuable insights into its possible density and structure. The success of such approach has been demonstrated by the recent advances in Sb phase change memory [20] where molecular dynamics (MD) melt/quench simulations have provided insights into how successful amorphization attempts of Sb can be.
A similar approach has been used to predict amorphous structures of other non-glass-forming oxides, such as HfO 2 and Al 2 O 3 [21,22]. The structure of bulk a-ZnO has also been studied using the melt/quench procedure and mainly using ab initio molecular dynamics (AIMD) methods [23][24][25][26][27][28]. Such calculations can provide an accurate description of interatomic interactions and electronic structures but use relatively small periodic cells and are several orders of magnitude more expensive than those using classical interatomic potentials (IPs) for describing interactions between atoms. The accuracy of AIMD calculations is often compromised by reducing the kinetic energy cutoff, using only the gamma point, NVT ensembles, and large time steps and, in some cases, quenching the structures too quickly. The reported a-ZnO properties are often based on a single atomic configuration and lack statistics required to describe distributions of structural characteristics.
More extensive studies can be carried out using classical IPs. Recently, Lin et al. [26] and Roy et al. [25] have modelled a-ZnO structures using IPs developed by Binks and Grimes [29], and Wang et al. [30], respectively. In the first study, simulated-annealing basin-hopping algorithm (with an initial density of 4.6 gcm −3 ) was employed. A very low a-ZnO density of 4.47 gcm −3 (a reduction by more than 21% with respect to the c-ZnO phase) is reported with an average coordination number (CN) of 3.71. In the second study [25], the lattice parameters and atomic positions of a 128-atom unit cell were extracted from a 4096-atom melt/quench simulation. These simulations report a much higher CN of 3.80 and mean O-Zn-O bond angle of 109.02 • (with the standard deviation of 25.7 • ). However, one structure is not sufficient to represent the distribution of a-ZnO structural properties and, as we show below, a 128-atom cell is too small to converge the structural properties of a-ZnO. Furthermore, the IPs used [30] overestimate both the experimental expansion coefficient by a factor of about three and the DFT surface energies for the two nonpolar ZnO surfaces-(1010) and (1120)-by about 50%, which may affect the a-ZnO topology. Therefore, there is scope for more extensive simulations of a-ZnO structures.
In this paper, we used IPs [31] and melt/quench simulations to create ensembles of amorphous structures and then DFT calculations to refine their structural characteristics. This allows us to investigate distributions of structural properties as a function of CR and cell size (see, for example, Ref. [22]). We then used the activation-relaxation technique (ART) [32][33][34][35][36] to explore the energy landscape of 20 independent a-ZnO structures. We found the a-ZnO structures to be stable with respect to local structural rearrangements. In only one structure out of 20, such rearrangements lead to significant structural changes, forming well-defined crystal motifs. These results suggest that stable a-ZnO samples could be produced under fast cooling from melt.
The paper is organized as follows. In the next section, we provide a detailed description of the computational methods and the results of testing two different IPs. We then continue by giving a detailed description of the atomic structure of a-ZnO. We discuss the short-and medium-range structural characteristics of a-ZnO across several samples and analyze different modeling parameters, including cell sizes and CRs.

II. COMPUTATIONAL METHODS
Over 500 a-ZnO structures were created using IPs and a melt/quench method, from which about 200 structures were fully optimized using DFT calculations. Details of these calculations are discussed below.

A. Generation of amorphous structures
Amorphous ZnO structures were generated using MD simulations and the LAMMPS code [37]. Two different sets of IPs were used for comparison: a simple force field that uses a Buckingham potential [29] and a more complex Born-Mayertype IP [31].
The latter has been split into different regions (with polynomial potentials used to spline the regions together), each of which deal with a specific coordination sphere of the opposed ions. This force field has been used to study the physical properties of crystalline ZnO, such as dielectric response, elastic, piezoelectric constants, and vibrational spectra, stability of surfaces, the phase diagram, and energetics of point defects in Refs. [38][39][40][41]. On average, the calculated ZnO static and high-frequency dielectric constants, elastic constants, and vibrational spectra were reported to be within ca. 2.94, 3.70, and 1.1% from the experimental data, respectively (see Table I in Ref. [31]).
To reduce the computational cost of MD simulations, a rigid ion model (removal of the spring constants between cores and shells in O species) of the IP reported by Whitmore, Sokol, and Catlow (WSC) [31] and Binks and Grimes (BG) [29] was used, with ionic charges of +2 and −2 for Zn and O, respectively. We note that the computationally less demanding rigid ion model of WSC's IP has been used to explore the energy landscape of ZnO clusters [42]. A similar two-step process is widely used where IPs are employed to filter low-energy candidates for a subsequent DFT refinement [41,43,44]. The main concern during the generation of amorphous structures is to obtain the correct topology. Inaccurate bond lengths and bond angles can be fixed during a full DFT optimization step. As one of the tests, we show below that rigid ion WSC IPs reproduce the experimental thermal expansion coefficients for ZnO with a good agreement with experiment.
To our knowledge, the density of a-ZnO is unknown. Therefore, a-ZnO structures were generated in a NPT ensemble with a Nosé-Hoover thermostat and barostat, using the following procedure: (i) equilibration of the system for 50 ps at 300 K; (ii) temperature was linearly increased to 5000 K for 50 ps; (iii) the system was equilibrated for 500 ps at 5000 K; (iv) structures were cooled down to 300 K with a CR of 100 Kps −1 ; (v) the system was equilibrated for 50 ps at 300 K. Cubic periodic cells were used in all cases. A range of CRs were tested, ranging from 0.75 to 800 Kps −1 , as discussed below.
Structural parameters of c-ZnO are similar for both WSC and BG's IPs. But, since the amorphous structures generated are essentially frozen melts, their density depends on whether the IPs used reproduce the experimental thermal expansion correctly. This was tested (see  [45] and yellow [46] lines. Blue lines represent the fit of the calculated lattice parameters with respect to temperature. A 5×5×3 supercell (300 atoms) of the wurtzite ZnO hexagonal lattice was used with the lattice parameters allowed to vary independently under zero external pressure. Both simulations (temperature increases/decreases) for each pair of IPs are similar, here we only display the runs where the temperature is increased from 50 to 2000 K; (a) WSC's IP; (b) BG's IP.
was by decreasing the temperature from 2000 to 50 K with the same rate. Similar behavior was obtained for the two separate runs. A 5×5×3 supercell (300 atoms) of the wurtzite ZnO hexagonal lattice was used with the lattice parameters, allowed to vary independently under zero external pressure. It can be seen that WSC IPs [31] show a good agreement with experimental data. The thermal expansion along the lattice parameter c obtained with the second IP [29] is, however, less satisfactory: above 1500 K, there is a very wide distribution of the c parameter with its average increasing constantly. This results in a much wider distribution of a-ZnO densities: The density standard deviation over 100 a-ZnO structures (produced with MD melt and quench with a CR of 100 Kps −1 and a 324-atom cell) is ca. 60% larger than for the WSC's potential. As we show below, the atomic structures generated by a rigid ion version of WSC IP [31] are in a good agreement with DFT calculations of a-ZnO structures. Therefore, these IPs are used in all further calculations discussed in this paper.

B. Exploration of the potential energy landscape
Thermal activation may lead to relaxation of amorphous structures decreasing the degree of structural disorder [47][48][49]. This effect can be studied using simulated MD annealing. However, structural relaxation often takes place on timescales that are orders of magnitude longer than afforded by MD simulations even using IPs.
Here, we have employed the ART [33,35,50,51] to sample the energy landscape of a-ZnO samples. This open-ended search method finds saddle points and adjacent minima around a relaxed configuration and it is capable of reaching energy landscape regions that involve the jump over high energy barriers (a task difficult to perform using more standard methods such as MD). ART has been used to study energy landscapes of amorphous Si [36,52] and SiO 2 [34].
Implementation of this method consists of three key steps: leaving the harmonic well, convergence to a saddle point, and relaxation into a new minima. The first step is carried out by picking a random atom. Then, this atom and its neighbors (within a 3.5 Å cutoff) are moved in random directions. Next, the system is pushed to the adjacent saddle point and converged to a new local minimum. The new minimum is accepted only if the energy difference between the two minima is below an energy threshold that is based on a Boltzmann weight. A Metropolis temperature of 0.25 eV has been used as a criterion to accept/reject new events (minima).
In the present paper, this technique has been used to search for new structures close to our initial amorphous structure. Structural parameters are also analyzed at different points of the energy landscape as well as the energetic cost for such changes. The potential energy landscape is described by WSC's IP [31]. Twenty independent ART runs were performed (enumerated from run1 to run20) using periodic cells containing 324 atoms. The initial configurations were taken from our MD melt/quench simulations at a CR of 100 Kps −1 . Ten-thousand trial events were generated for each run, with a success rate of ca. 82% (ca. 164 000 accepted events out of 20 000). New minima and saddle points for each event were stored and analyzed. Lattice parameters of every accepted event were locally optimized and structural parameters calculated, including density, average CN, average Zn-O-Zn angle bond, and standard deviation of the bond angle distribution.

C. Periodic ab initio calculations
All the electronic structure calculations were performed using the Gaussian and plane-wave method implemented in the CP2K code [53] employing the generalized gradient approximation (GGA) Perdew-Burke-Ernzerhof (PBE) functional.
A kinetic energy cutoff of 600 Ry (8163 eV) was sufficient to converge the bulk lattice energy (four atoms unit cell) to less than 1 meV.
The DZVP-MOLOPT-SR-GTH basis sets [54] were employed in all atomic species with the GTH pseudopotentials [55]. When an atomic relaxation was performed, the forces on all ions were converged to less than 0.02 eVÅ −1 . A 300-atom unit cell was used for the crystalline phase. The lattice constants for crystalline ZnO in wurtzite structure obtained both using IPs and DFT are in a good agreement with experiment [56][57][58] and theory [38,39]--see Table I. Structural parameters are reproduced within a 2% and 0.5% 014202-3 for PBE and IP, respectively. All crystal structures in this paper were generated using the VESTA package [59]. All the figures in this work have been produced using GNUPLOT [60].

III. RESULTS AND DISCUSSION
The melt/quench method using classical IPs [31] in the NPT ensemble was employed to create more than 500 amorphous structures across eight different cell sizes ranging from 96 to 768 000 atoms. The lattice parameters and atomic structure of about 200 systems, containing 96 and 324 atoms, were fully optimized using DFT and the PBE functional. In modeling amorphous systems, small simulation cells exhibit spurious long range order due to periodic boundary conditions. This can be counteracted by using larger cells including up to hundreds of thousands of atoms. Although feasible when using classical forcefields, these are too computationally expensive for DFT calculations of the electronic structure and structural properties. Importantly, the optimal CR and cell size discussed below is practical only for melt/quench simulations using IPs. To obtain good statistics and fully DFT optimize the geometry of amorphous structures, a compromise cell size is 324 atoms. As we show below, using much larger cells does not change significantly the a-ZnO structural parameters obtained.
The melt/quench procedure used here consists of five steps: equilibration at room temperature, ramping to the melting point, melting the material, cooling down, and equilibration at room temperature. Below, we provide a detailed description of the generation of the a-ZnO structures.

A. Melt cooling rate and glass formation
Assuming that the melt has been equilibrated (in this paper, this was achieved by melting the material for 500 ps), the CR and the size of the simulation cell control the topology of final structure. In previous simulations of other oxides, e.g., SiO 2 [61], Al 2 O 3 [21], HfO 2 [22], CRs range from 0.75 to 8 Kps −1 and cell sizes contained from 300 to 600 atoms. In ZnO, such low CRs, however, frequently generate highly ordered structures. As mentioned above, ZnO shows a strong trend to form crystalline phases, particularly the wurtzite structure. Other thermodynamically stable zincblende and rocksalt structures (+0.021 eV and +0.179eV per ZnO pair, when compared to ZnO wurtzite [38]) are rarely seen in experiment. Since growing a-ZnO films is so difficult, establishing a CR required for achieving stable amorphous phase is a critical question.
To study the effect of CRs on the structures produced, we have tested 14 different CRs spanning two orders of magnitude (from 0.75 to 800 Kps −1 ) and eight simulation cell sizes containing from 96 up to 768 000 atoms. The effect of CR was analyzed mainly using a 324-atom unit cell. A much bigger 6144-atom unit cell has also been used for tests with similar results. As a measure of how energetically distant the amorphous structure is from the lowest energy crystal structure, we use the so-called excess energy ( E)-the energy difference per ZnO unit of the amorphous structure and the wurtzite structure. This energy is stored in the amorphous network and its value with respect to crystal in some cases, such as a-Si, has been measured experimentally using calorimetry [62,63]. Densities of amorphous films are usually lower than those for crystals. In this paper, therefore, the effect of the CR was monitored by comparing densities, Zn and O CNs and E for obtained structures. Figure 2 shows the density distribution and the average density with respect to the CR. Twenty different runs per CR were used. We observe that CRs smaller than 100 Kps −1 tend to produce (semi)crystalline structures, with densities close to the wurtzite phase (see Fig. 2) containing interstitial ions. Whereas for CRs exceeding 100 Kps −1 the obtained structures are amorphous with no long-range order. Such fast CRs have been recently used to form monatomic metallic glasses using ultrafast laser quenching [65].
To investigate whether experimentally more accessible CRs could produce amorphous structures, we have analyzed the structures produced using CRs of 1, 10, and 50 Kps −1 using a cubic 324-atom cell and 20 structures per CR. At CR = 1 Kps −1 , 100% of obtained structures are crystalline. The number of crystallized structures drops to 60% and 5% for faster CRs of 10 and 50 Kps −1 , respectively. We note that constrains imposed by the cubic cell prohibit the creation of wurtzite or zincblende polymorphs at CR = 1 Kps −1 . Therefore, interstitial O ions are formed, as can be seen in Fig. 2. However, tests on 320-atom (a multiple of 8, the number of atoms in ZnO zincblende) cubic cells produced crystalline structures with no defects.
By increasing the CR by one order of magnitude (10 Kps −1 ), we observe a cell size dependence (discussed in more detail below); bigger cells are needed to converge the structural parameters. For 324 atoms, more than half of the samples are crystalline with several point defects, whereas bigger cell sizes are amorphous but less disordered than those produced using 100 Kps −1 .
At CR = 50 Kps −1 , the density is on average reduced by ≈10.38% to 5.09 gcm −3 with respect to the crystalline phase. In comparison, densities of other amorphous oxides, such as SiO 2 , Al 2 O 3 , HfO 2 [21,22,61], obtained using a similar melt/quench procedure and NPT ensemble are by 12−15% smaller than the corresponding crystalline values. With respect to CRs equal to 100 Kps −1 , these 50 Kps −1 structures have a smaller and higher percentage of 3-and 4-coordinated ions, respectively: the percentage of 3-, 4-, and 5-coordinated ions is approximately 11, 86, 3%, respectively (compared with 12, 85, 3%, Table II). Moreover, another sign of a more ordered structure is the average E, which is 8% lower than in 100 Kps It turns out that CR exceeding 100 Kps −1 produce highly unstable and thermally sensitive structures. These structures usually have small energy barriers to undergo considerable structural changes. We have analyzed the structure and energetics of 20 324-atom structures produced at 400 Kps −1 and compared them with those of a CR of 100 Kps −1 . Even though the average density difference between CR = 400 Kps −1 and 100 Kps −1 is within less than 1%, these structures have a 8% larger E (0.54 eV per ZnO pair) than those for 100 Kps −1 (see Table II). There is also a discrepancy in distribution of CNs: the percentage of 3-, 4-, and 5-coordinated ions is 15, 82, 3%, respectively (compared with 12, 85, 3%, Table II ). At the higher CR, there is also a slightly broader bond angle distribution than that produced at CR = 100 Kps −1 (see Fig. 3).
The higher instability for the 400 Kps −1 configurations was confirmed by simulated MD annealing of five structures at 600 K for 50 ps. On average, the density of annealed 400 Kps −1 structures increases by 1.81%, while it is only 0.90% for five 100 Kps −1 structures. However, the main difference is the change of the atomic structure. After annealing, the 400 Kps −1 structures change substantially with a high degree of crystal motifs, whereas the morphology of the 100 Kps −1 samples remains the same.
These results demonstrate that low CR produce crystalline structures, whereas high CR do not allow the system to relax and the obtained structures keep many highly unstable motifs frozen from the melt. We find that a CR of 100 Kps −1 is fast enough to prevent the system to attain a crystalline structure, but at the same time, is slow enough to allow    highly unstable motifs from the melt to relax. All amorphous structures discussed below were, therefore, created with the CR of 100 Kps −1 . This CR has also been used in some of the previous simulations of a-ZnO [23,27].

B. Atomic structure of a-ZnO
Using the selected CR, we continue with the analysis of structural parameters for 324-atom cells produced by a melt/quench procedure. The geometries of 100 a-ZnO structures were fully optimized using DFT with the PBE functional. CNs, radial, and bond angle distribution functions were calculated using the R.I.N.G.S. code [66] with a cutoff of 2.45 Å. The radial cutoff was chosen as the back of the first peak of the radial distribution function graph (see Fig. 4).
The average radial distribution and partial pair distribution functions, from 100 amorphous structures, are shown in Fig. 4. The RDF profiles of a-ZnO structures produced using IPs and after DFT optimization are very similar. There is only an adjustment of the Zn-O bond length within the cell (Fig. 5) but the structures and their pair distribution functions (Fig. 4) and bond angle distribution (Fig. 3) remain the same. We note that the first coordination shell for the cation is very narrow, preserving the ZnO bulk character. This characteristic short-range behavior has been observed in other amorphous oxides, see, e.g., Refs. [28,67,68].
The maxima of RDF peaks of the second-shell Zn-Zn and O-O distances match closely those for c-ZnO. The distributions are, however, much wider than in the case of the first shell. The profiles are in good agreement with the AIMD results by Pandey et al. [23] The reduction of density with respect to wurtzite ZnO calculated using IP and DFT is 11.27 and 8.3%, respectively.  Fig. 3), with Gaussian-like distribution, at ca. 89 • and 112 • . The second peak resembles the one in the wurtzite structure, whereas squarelike structures with right angles are present in rock salt ZnO, one of the ZnO polymorphs seen experimentally at high pressures. In rock salt, however, ions are 6-coordinated, a characteristic that is not seen in a-ZnO.
In general, most short-range configurations are well defined, 4-coordinated, and with an environment similar to that in c-ZnO. Medium-range structures are, however, highly affected by disorder with a wide distribution of their structural properties. The medium-range order concerns the spacial distribution of the tetrahedra and the connections between them, e.g., corner-sharing, edge-sharing, or face-sharing. Considering the Coulomb interaction, these connections have different energetic implications. For c-ZnO, corner-sharing of tetrahedra is the most energetically favorable configuration with a cation-cation distance of ca. 3.2 Å (ca. 60% longer than the anion-cation bond distance), whereas for the perfect edgesharing polyhedra is ca. 2.2 Å. Face-sharing configurations have the shortest Zn-Zn distances and are the most energetically expensive. In the case of the edge-sharing tetrahedra, however, the strong Zn-Zn Coulomb repulsion increases the 2.2 Å distance, as observed in our a-ZnO structures-the average Zn-Zn distance for edge-sharing tetrahedra is 2.87 Å (see Fig. 6). This distance matches the Zn-Zn and O-O shoulder of the main peak in Fig. 4. The relative height of this shoulder and the main peak suggests that the amount of edgesharing tetrahedra is approximately half of the corner-sharing tetrahedra. Face-sharing polyhedra are highly unstable and are not seen in our a-ZnO structures.
The average excess energy per ZnO pair relative to the crystalline phase calculated using IP and DFT is 0.5 eV and 0.32 eV, correspondingly (see Table II). These values are similar to those reported using AIMD by Walsh et al. [67] (0.352-0.376 eV) for more complex ZnO systems-In 2 O 3 (ZnO), InGaZnO 4 , InAlZnO 4 -and by Pandey et al. [23] for a-ZnO (0.02 eV per atom difference for two different models).

C. Density of a-ZnO and effect of simulation cell size
The results of the previous section show that DFT geometry optimization does not change significantly a-ZnO structures obtained using IPs. Therefore, we use IPs to model bigger structures, where DFT calculations are unfeasible. Figure 7 shows the normalized density distribution of the a-ZnO structures using NPT with different cell sizes and a CR of 100 Kps −1 . First, we note that large fluctuations in small cells create a very broad distribution of densities. The dispersion in density decreases as the volume increases: For example, a 324-atom cell has a standard deviation of 0.0494 gcm −3 , while for a 96 000-atom cell it is only 0.0030 gcm −3 . One can see that for the largest cells, the density converges to about 5.04 gcm −3 .  Table II. It is interesting to compare density distributions, shown in Fig. 7, with a different approach where 324 000-and 768 000atom cells are divided into a 3D grid of smaller supercells with equal volumes and density is calculated for each cell. We note that the number of atoms in each cell may deviate from stoichiometry due to the rigid grid. Nevertheless, the density distribution reflects the local density fluctuations in different regions of the large cell. Semitransparent lines added in Fig. 7 for the 324, 768, 1500, and 6144 profiles were obtained as follows. The 324 000-atom cell was geometrically divided into 1000 equivalent cells. The remaining three distributions were obtained by dividing the 768 000-atom cell into 1000, 512, and 125 equivalent cells. As a result of variations in the number of atoms in each cell and local density fluctuations, we observe distributions of densities. However, due to the large number of cells, these distributions are much smother than the ones obtained for smaller ensembles of cells above and exhibit Gaussian-like behavior for all added distributions. The average a-ZnO density obtained in the 768 000 atoms cell is 5.04 gcm −3 . We also note that for all different 324 000-and 768 000-atom samples, the RDF and bond angle distribution profiles are very similar. The density distribution within those samples is also very similar.
Pair distribution functions are show in Fig. 8(a). As noted above, a periodic cell size of 96 atoms is not enough to reproduce the average properties of the a-ZnO system. There are peaks after 3 Å, which correspond to (semi)crystalline samples. Table II displays the average CN for the Zn-O pairs as a function of cell size. The average IP CN for a volume bigger than a 32-atom cell changes no more than 0.5% and the average density difference fluctuates within 0.01 gcm −3 . Moreover, values obtained using IPs show an excellent agreement with DFT calculations (numbers in parenthesis). We note that the average CN for O and Zn after DFT optimization for the smallest cell size (96 atoms) is the same as reported by Pandey et al. [23], who used an melt/quench AIMD procedure  Table II for both figures. and a cubic cell containing 128 atoms, a CR of 100 Kps −1 , and a PBE functional. Moreover, the agreement among average properties for a 324-and a 768 000-atom cell size is also good.
Finally, we checked how the obtained structures depend on the cell size for different CRs. Three CRs (1, 10 and 100 Kps −1 ) were used for cell sizes containing 324, 768, 1500, and 6144 atoms [see Table III and Fig. 8(b)]. At 100 Kps −1 , all structures are amorphous and the average structural properties have converged using a 324-atom cell. At 10 Kps −1 , there is a cell-size dependence: more than half of the structures show high degrees of crystallinity using cells with 324 atoms. Bigger cells produce amorphous structures, but less disorder than those obtained using 100 Kps −1 . At 1 Kps −1 , every structure is highly ordered independent of the cell size used.

D. Exploration of the a-ZnO Potential Energy Landscape
The results presented above suggest that bulk a-ZnO structures could be produced from melt using appropriate CRs. Difficulties in producing such structures experimentally suggest that they may be unstable and prone to crystallization. To explore the energy landscape of a-ZnO, we have employed the event-based structural-evolution code ART. This approach has the advantage of exploring energy landscape areas far from the initial local minima, which is difficult to access for MD techniques due to the limited timescales. By perturbing a-ZnO structures by moving atoms into different local minima, we explore their stability. Moreover, we report the structural parameters of a-ZnO at different points of the energy landscape and the energetic cost for such structural changes.
ART has been used to find local transition states and analyze their energetic barriers in amorphous Si and SiO 2 [34,36,52]. Each successful event corresponds to a structure changing its initial atomic configuration for another metastable structure. These are strong perturbations of the structure and some of these moves require activation energies of several eV. As a result of many such moves, the structure is effectively reshuffled. Figure 9 shows the E profile for each ART accepted event of 20 independent runs (run1 to run20) with different initial structures. Density, CN, and standard deviation of the Zn-O-Zn bond angle distribution are also shown. The 20 ART runs can be grouped in three categories when comparing changes in their structure and energetics: no significant changes (18 runs), small changes (run1), and significant changes (run6). Both run1 and run6 are labeled in Fig. 9. There is strong correlation with a Pearson correlation coefficient [69] of over 0.95, for the run 6, between energies and density, CN and standard deviation of the Zn-O-Zn bond angle distribution. There is a weak or nonexistent correlation among these properties for the other runs.
Since the structural distortion to find a new event is random, we have picked two ART runs with contrasting behavior: run6 and run15. The initial configurations of both runs are again evaluated five more times, each with the ART approach following the same steps as the ones described in the Methodology section: a total of six ART runs per initial configuration. This helps us to identify whether a structure is stable/unstable or it is an artifact of the random process. Run6 and run15 were selected because they show significant and insignificant structural evolution, respectively.
From run6 calculations, we note that every ART simulation results in a highly crystalline ZnO structure. The potential energy profiles with respect to an accepted event are, however, different. Minor structural and energetic differences were found among the six well-relaxed ZnO structures. On average, E ≈ 0.10 eV per ZnO pair was calculated for the final structures.
Run15 behaves slightly differently, though. At the end of the ART simulation, five out of six runs share very similar structural properties with total energies being within 7 meV per ZnO pair and an average E of ≈0.47 eV per ZnO pair. One of the ART runs, however, decreases its E to ≈0.34 eV. This decrease in energy is accompanied by creation of some crystal-like motifs, which are not that common in the other five configurations. We note that this structure still remains amorphous with a relatively high E. It is very unlikely for run15 initial configuration to reach the energy landscape area where structures have an E per ZnO pair of ≈0.34 eV. The reason for this is that only one out of six ART simulations run6 with only a few hundred ART accepted events. It is clear, therefore, that run6 has some structural features that are not present in run15 (or any other run), which helps to nucleate crystallization, even though their initial average properties are similar.
Following the latter observations, we have analyzed the characteristics of initial 20 structures (run1 to run20): E, RDF, pair radial distribution, CN, average Zn-O-Zn bond angle, Zn-O-Zn bond angle standard deviation, and X-ray diffraction. There is no clear energetic or structural difference among the 20 initial structures. We note, however, that there is small difference in the O-O pair distribution function profile of run6 when compared to the others. In every case but in run6, there are two well-defined peaks at ca. 2.8 Å and ca. 3.3 Å, which correspond to edge-sharing and corner-sharing tetrahedra, respectively. In run6, however, instead of a welldefined peak at 2.8 Å, there is a shoulder that overlaps with the main 3.3 Å peak. This suggests that the number of edgesharing tetrahedra in run6 is less and this may be helping the crystallization.
To gain a further insight into this phenomena, we carry out a simulated MD annealing process at 400 K for 25 ns on run6 and run15. Other annealing temperatures (600, 800, and 1000 K) were also tested with 400 K being the lowest temperature needed to crystallize run6. There is increase of more than 7% in density for run6, whereas there is no significant changes for run15. In contrast to run15, run6 shows a high degree of crystallinity. We conclude that amorphous structures with average structural and energetic properties have different medium-range motifs that trigger crystallization.
Finally, we note that structural relaxation is a process that takes place on timescales that are orders of magnitude larger than MD time steps. In the case of run6, up to 10 7 MD steps (nanoseconds) at 400 K are required to obtain the (semi)crystalline structure. It took, however, only hundred of ART events to achieve a similar effect. The ART moves represent the most likely physical trajectories followed during structural relaxation. Moreover, this method is much less affected by the computational cost increase caused by high activation energy barriers than MD techniques [34].

IV. SUMMARY AND CONCLUSIONS
To summarize, the results of our calculations suggest that stable a-ZnO structures should be possible to produce by cooling of ZnO from melt. The stability of amorphous ZnO against crystallization at room temperature strongly depends on the rate at which it was cooled from the melt and the rates around 100 Kps −1 provide stable structures. We have shown that reliable predictions can be obtained using IPs and large simulation cells, where DFT calculations are currently unfeasible. The predicted DFT density of such structures is about 4.97 gcm −3 and excess energy about 0.3 eV per ZnO unit. ART and MD annealing analyses demonstrate that such amorphous structures are stable and do not crystallize easily. Unlike in the case of glass-formers, such as SiO 2 [70] or silicates [71], there is no convergence of the amorphous structure properties with decreasing the CR. Just the opposite-at CR lower than 100 Kps −1 , our samples tend to crystallize forming defective crystalline structures, in agreement with experimental observations. Such high CRs can be achieved by using ultra-fast laser quenching [65] and have been used to form monatomic metallic glasses, which are very difficult to form otherwise. Strain at film-substrate interface as well as impurities and dopants, such as Sn, Ge, and others [72,73], improve stability and performance of amorphous ZnO based films.
Our simulations provide strong support to the notion that, by quenching ZnO rapidly enough, one can succeed in creating a glass of pure ZnO with sufficient stability against crystallization even around room temperature. The fact that creating such systems is becoming increasingly possible has been demonstrated by the recent advances in Sb phase change memory [20]. Further studies will help to identify optimal conditions to forming stable a-ZnO films under deposition and anneal conditions.