Generation and properties of bulk a-ZrO2 by molecular dynamics simulations with a reactive force field

    Generation and properties of bulk a - ZrO 2 by molecular dynamics simulations with a reactive force field S. Arash Sheikholeslam, Guangrui Maggie Xia, Cristian Grecu, Andr e´ Ivanov PII: DOI: Reference:

S0040-6090(15)01002-0 doi: 10.1016/j.tsf.2015.10.026 TSF 34716

To appear in:

Thin Solid Films

Received date: Revised date: Accepted date:

25 January 2015 8 October 2015 9 October 2015

Please cite this article as: S. Arash Sheikholeslam, Guangrui Maggie Xia, Cristian Grecu, Andr´e Ivanov, Generation and properties of bulk a - ZrO2 by molecular dynamics simulations with a reactive force field, Thin Solid Films (2015), doi: 10.1016/j.tsf.2015.10.026

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

PT

Generation and properties of bulk a-ZrO2 by molecular dynamics simulations with a reactive force field

SC

RI

S. Arash Sheikholeslama,∗, Guangrui(Maggie) Xiab , Cristian Grecua , Andr´e Ivanova a

MA

NU

Department of Electrical and Computer Engineering, the University of British Columbia, Vancouver, BC, V6T 1Z4 Canada b Department of Materials Engineering, the University of British Columbia, Vancouver, BC, V6T 1Z4 Canada

Abstract

AC CE P

TE

D

A classic molecular dynamic melting and cooling approach was used to generate bulk amorphous ZrO2 using a reactive force field. The density and the structure characteristics of the generated bulk ZrO2 were analyzed and compared with literature simulations and experimental data. The density dependence on the cooling rate in the simulations was established. The generated bulk ZrO2 can further be used for theoretical studies of amorphous ZrO2 properties. Keywords: amorphous zirconium dioxide, molecular dynamics, ReaxFF, annealing 1. Introduction

With the gate dielectric of metaloxidesemiconductor field-effect transistors approaching dimensions below 1 nm [1], there has been a surge in the research and applications of high-k dielectrics such as Hf and Zr oxides as they can meet both the scaling and coupling capacitance requirements. ZrO2 is considered as one of the leading candidates for alternative gate dielectrics [1]. Amorphous ZrO2 (a-ZrO2 ) has been in use in memory applications and in high electron mobility transistors [2, 3, 4]. Other applications of ZrO2 include the protection of nuclear fuel rods against corrosion and fuel cell ∗

Corresponding author Email address: [email protected] (S. Arash Sheikholeslam )

Preprint submitted to Thin Solid Films

October 12, 2015

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

electrolytes [5]. Understanding the chemical and structural characteristics such as density, energy bandgap, density of states, dielectric constant, structural defects, and diffusion behaviour is therefore very important. Only a handful of experimental studies of a-ZrO2 characteristics are available [6, 7]. Quantum chemistry methods such as density function theory (DFT) have been employed previously to simulate a small supercell of a-ZrO2 consisting of tens of atoms [8, 9, 10]. A Molecular Dynamics-Density Functional Theory (MD-DFT) simulation of ZrO2 was performed in [8], where cells of around 100 atoms were formed. Classical annealing followed by DFT annealing were performed and some geometrical characteristics of the resulting oxide have been analyzed [8]. In [10], a melt and quench approach has been used over a time window of 22.5ps. One limitation of this latter work is the fixed volume of the simulated cell, which means that expansion and contraction of the simulated cell could not be performed. Another general issue with quantum chemistry methods is that they are limited by the computational complexity compared to MD. For instance, hydrocarbon reactions performed by reactive force field (ReaxFF) molecular dynamics were shown to be about 10000 times faster than a quantum chemistry method [11]. The necessity to generate a-ZrO2 structures with a large number of atoms is amplified as one approaches problems such as doping, defects, and diffusion behaviour in bulk materials. In this work, classical molecular dynamics simulations were performed using Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software [12, 13] and ReaxFF [14, 11]. ReaxFF was introduced in [11] as a reactive force field for hydrocarbons which is capable of simulating bond forming and breaking. Classical annealing has been effectively used in forming SiO2 glasses both with Beest-Kramer-van Santen potentials [15] and ReaxFF force fields [16]. Due to its derivation from ab-initio simulations, ReaxFF enables forming a-ZrO2 with reasonable computational resources while staying within the realm of quantum mechanical accuracy. Our method employed a larger simulation cell of 1500 atoms (We have performed a similar work with a simulation cell consisting of 2000 atoms. No noticeable difference in the characteristics was observed). Classical annealing processes were used to generate realistic amorphous ZrO2 cells. The densities of the generated cells were compared with available simulations and experimental values. Structural properties such as bond lengths, bond angles, radial-distribution Function (RDF) were also calculated.

2

ACCEPTED MANUSCRIPT

PT

2. Validation of ReaxFF for ZrO2 systems

AC CE P

TE

D

MA

NU

SC

RI

As explained earlier, classical MD was used to form our a-ZrO2 cells. The force field that we used in this work was developed first to study BaZrO3 proton conductors [17]. This force field contains all of the elements needed for this work. To validate the use of this force field, we calculated the cohesive energy as a function of cell volume for five different ZrO2 crystal structures: Monoclinic, Cubic, Tetragonal, OrthorhombicI, and OrthorhombicII [18, 19, 20, 21, 22]. We made supercells, each consisting of 6 × 6 × 6 unit cells and computed the cohesive energy while changing the cell volume. To be comparable with [23], we scaled our results accordingly (Figure 1). From Figure 1, the energy ordering and values of the five crystal structures EM ono < EOrthoI < ET etra < ECube < EOrthoII are consistent with the DFT calculations given in [23].

Figure 1: Unit cell binding energy versus lattice constant for five ZrO2 structures.

3. Methods

We started with a supercell of monoclinic ZrO2 crystal consisting of 1500 atoms [18]. The dimensions of this cell were X = 25.75˚ A, Y = 25.72˚ A and 3 Z = 26.58˚ A. The initial density was 5.68 g/cm . An annealing process was used to form the a-ZrO2 cells. There are two variable parameters in our simulations: the heating temperature and the cooling rate. Since the simulations were performed over time periods considerably shorter than those that would be associated with real experiments, the heating temperatures should be higher than the a-ZrO2 melting point. In this work, 4000K, 6000K and, 8000K were chosen as heating temperatures. Three different cooling rates of 1015 K/s, 1014 K/s and, 1013 K/s were initially 3

ACCEPTED MANUSCRIPT

MA

NU

SC

RI

PT

used in this work to obtain the basic characteristics. This means a total of nine systems were simulated. After the preliminary results were obtained, other cooling rates were used to establish a quantitative relation between the cooling rate and the final density. Initially, randomization of the initial monoclinic crystal [18] was performed, using an NVT ensemble, where NVT denotes a system with constant number of particles, volume, and temperature. Temperature is regulated via a Nose-Hoover thermostat [24], [25]. After the initial stage, an NPT heating and cooling was performed. NPT is similar to NVT but with constant pressure (P) and variable volume. Note that the pressure in all of our NPT simulations was fixed at 1 atm. An NPT ensemble will allow the system to expand and contract (similar to most of the experimental processes) but one cannot trivially control the final density. The process was performed through isothermal annealing of the system and then linearly cooling it down to the desired temperature. The annealing procedure was as follows:

AC CE P

TE

D

1. A monoclinic crystal simulation cell was initiated at 4000K and linearly cooled down to 300K using an NVT ensemble. The cooling rate at this step does not matter as our main goal in this stage is to randomize the system and remove any dependence on the initial structure. An NVT ensemble was used so that the initial density of the cell remains fixed. Time steps of 0.2f s were used for this part. Time steps of 0.1f s were used for the rest of the annealing process. A total of nine simulation cells were initialized in this manner. 2. Each simulation cell was again heated up to the heating temperature, using an NPT ensemble for 1ps. 3. The simulation cell was kept at the heating temperature using an NPT ensemble till the system was totally melted (i.e, the coordination numbers distribution changed to a more uniform distribution). In this simulation the samples were heated for t > 40ps. 4. The simulation cell was linearly cooled down back to 300K using an NPT ensemble. 5. The simulation cell was relaxed at 300K using an NPT ensemble for 1ps. Figure 2 depicts the annealing process as explained above. To better visualize the system during melting and cooling, Figure 3 shows the evolution of Zr and O coordination number distribution during the melting and cooling processes. We chose the 8000K data to generate these results. A comparison 4

SC

RI

PT

ACCEPTED MANUSCRIPT

NU

Figure 2: Annealing process for three anealing temperatures and cooling rates.

MA

between Figure 3(a) and 3(b) show that Zr and O atoms rearranged during the melting. The arrangement reached a steady-state as shown in Figure 3(b) and 3(c). 4. Results

AC CE P

TE

D

Slices of the generated bulk a-ZrO2 after the above steps are shown in Figure 4a-4c. As expected, a faster cooling rate generated lower density cells. The major structural characteristics of the resulting systems are given in Table 1. The resulting RDF’s are plotted in Fig. 5 for 8000K heating Cooling rate Heating temperature Density(g/cm3 ) Zr-O bond length (˚ A) A) RDFO−O first peak (˚ A) RDFZr−O first peak (˚ RDFZr−Zr first peak (˚ A)

4000 4.84 2.189 2.756 2.143 3.368

1015 K/s 6000 4.83 2.182 2.808 2.033 3.389

8000 4.84 2.185 2.676 1.985 3.367

4000 5.07 2.192 2.817 2.090 3.544

1014 K/s 6000 5.07 2.194 2.672 2.119 3.409

8000 5.09 2.190 2.818 2.091 3.545

4000 5.13 2.192 2.824 2.095 3.370

1013 K/s 6000 5.18 2.202 2.620 2.232 3.396

8000 5.12 2.193 2.850 2.155 3.545

Table 1: Comparison of the basic structural properties of results for various temperatures and heating temperatures.

temperatures. The amorphous nature of the results is implied by the RDF plots. From the data we obtained there seem not to be a clear relation between heating temperatures, cooling rates, and the location of RDF peaks. Our results for Zr − O bond lengths are depicted in Fig. 6. The bond lengths slightly increase as the cooling rates decrease and they are independent of heating temperature (Table 1 and Fig. 6). We also plotted the O − Zr − O and Zr − O − Zr angles. Unlike silica, the bond angles are not scattered in a Gaussian form around a peak (Fig. 7 and Fig. 8). They however seem to be 5

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Figure 3: Coordination number distribution of Zr and O during the melting and cooling process: (a) T = 8000 K, t = 0.2 ps; (b) T = 8000 K, t = 16.7 ps; (c) T = 8000 K, t = 66.7 ps; (d) T = 6000 K, t = 267.2 ps; (e) T= 4000 K, t = 459.7 ps; (f) T = 8000 K, t = 652.2 ps; (g) T=300 K

6

K

and (c) T=8000 1014 K/s.

K

and

MA

NU

(a) T=8000 K and 2.57× (b) T=8000 1015 K/s. 1015 K/s.

SC

RI

PT

ACCEPTED MANUSCRIPT

K

and

TE

D

(d) T=8000 1013 K/s.

AC CE P

Figure 4: Sliced section of the final samples for various cooling rates. The depth of these slices is 4.5˚ A. Images were produced with OVITO [26] visualization software.

Figure 5: The RDF of the final system for the three cooling rates and 8000 K heating temperature.

scattered in a two-Gaussian function. Therefore, for the probability density of the bond angles we have: P r(θ) = αe

(

−(x−θ1 ) 2 ) σ1

+ βe

(

−(x−θ2 ) 2 ) σ2

(1)

This means that the angle data is scattered around two values (θ1 and θ2 ) instead of one and σ1 and σ2 are the standard deviations for those values. Ta7

SC

RI

PT

ACCEPTED MANUSCRIPT

NU

Figure 6: The bonds length distributions for three different cooling rates and heating temperatures.

AC CE P

TE

D

MA

ble 2 shows the parameters of our two-Gaussian fit for three different heating temperatures at 1013 K/s cooling rate. It can be observed that for O −Zr −O both θ1 and θ2 slightly decrease by increasing the heating temperature while for Zr − O − Zr θ1 and θ2 slightly increase.

Figure 7: The O −Zr −O angle distribution for the three cooling rates and 8000 K heating temperature. A two-Gaussian is fit on 1013 K/s cooling rate data.

Heating temperature Analyzed angle α E1 σ1 β E2 σ2

4000K OZrO ZrOZr 3.371 5.592 76.47 100.6 12.93 9.272 1.472 1.974 119.9 127.8 41.65 25.62

6000K OZrO ZrOZr 3.303 5.68 75.67 100.9 12.53 10.46 1.48 1.774 118.7 130.9 42.93 23.9

8000K OZrO ZrOZr 3.341 6.029 75.58 101.3 12.24 9.852 1.502 1.885 118.2 131.3 42.84 22.57

Table 2: Parameters of the two-Gaussian fit of the angle data.

8

SC

RI

PT

ACCEPTED MANUSCRIPT

NU

Figure 8: The Zr − O − Zr angle distribution for the three cooling rates and 8000 K heating temperature. A two-Gaussian is fit on 1013 K/s cooling rate data.

AC CE P

TE

D

MA

One reason for this two-Gaussian shape of the bond angle distribution (in contrast with silicon dioxide which is a Gaussian) is the zirconium-oxygen multiple bonds; i.e; unlike the silicon dioxide where the coordination numbers are fixed at four for silicon and two for oxygen, here there are various coordinations for each atom (Table 3). This causes some of the bonds to be bent more than others. One can observe that there is no significant relation between the final density and the heating temperature. This is expected to be true as long as the heating temperature is well above the melting point of the initial crystal. The glass transition temperature for different cooling rates can be obtained as shown in Figure 9, which shows the volume of simulation cell during the cooling process for an heating temperature of 6000K and cooling rates of 1015 K/S and 1014 K/S. The glass transition temperature (and the density)

Figure 9: cell volume versus temperature during the cooling from 6000K heating temperature with two different cooling rates.

of the annealed simulation cell is a function of the cooling rate. In fact, a 9

ACCEPTED MANUSCRIPT

1 ) + log10 (5 × 1016 ))−5.22 + 5.149 (2) cooling rate

MA

Density = −4.229(log10 (

NU

SC

RI

PT

slower cooling rate results in a denser final structure. Densities of various simulations are given in Table 1. All of the obtained densities are within the accepted theoretical density range for amorphous zirconium dioxide [9, 10]. In order to further analyze the relation between the cooling rate and final structure density, one more simulation data point was added at 2 × 1015 K/s cooling rate. Figure 10 shows the density of annealed simulation cell versus cooling rates for the three different heating temperatures. A curve is fitted on 8000 K data using a simple power function of the form f (x) = axb + c and we found the parameters using nonlinear least square method (where x is the shifted logarithm of the inverse rate). The result is the following function:

AC CE P

TE

D

Based on Equation 2, as the cooling rate approaches zero, the density approaches 5.149 g/cm3 which is not very different than the densities we have obtained using 1013 K/s cooling rate. We can also choose the density of the final system by varying the cooling rate using Equation 2.

Figure 10: This plot and the fitted curve shows the relation between the cooling rate and the final density of the cell.

4.1. Comparison with experimental a-ZrO2 density To demonstrate the use of Equation 2, we reproduced the density of an experimental study of a-ZrO2 [7], which used electron beam deposition to form a thin zirconia film sample of 4.2g/cm3 .This density is well bellow the theoretically accepted numbers and is an unstable form of a-ZrO2 according to [10]. However, we tried to reproduce this system using Equation 2. Substituting 4.2 g/cm3 for our density, this yields a cooling rate 10

ACCEPTED MANUSCRIPT

SC

RI

PT

of 1015.3633 K/s ≈ 2.3083 × 1015 K/s. This means that the cooling must be performed over 33333 steps of 0.1f s. We adjusted the number of steps to 30000 (≈ 2.57 × 1015 K/s) and performed an 8000 K heating process using this cooling rate. We achieved a density of 4.18 g/cm3 . We refer to the cell with this density as the low density cell.

MA

NU

4.2. Coordination analysis and comparisons with DFT simulations In this section, we compare the results from the two previous works on smaller simulation cells using DFT [8, 10]. ReaxFF parameters are optimized to achieve geometry and energy characteristics computed by plain wave quantum chemistry methods. Therefore, DFT calculations are good benchmarks of quality for our simulations. O and Zr coordination distributions are given in Table 3 which contains the coordination distribution results for four different cooling rates including the low density (4.18 g/cm3 ) sample. It also contains the data from [8, 10]. The main observations are:

AC CE P

TE

D

1. Table 3 shows that a slower cooling rate results in a higher seven fold Zr coordination (Zr(7)) as well as higher O(3) and O(4). This is somewhat the case with Zr(8) with only one exception. Crystalline ZrO2 (monoclinic and tetragonal phase) has Zr(7), Zr(8), O(3) and O(4). This is consistent with the observed trend. 2. Our coordination results are closely consistent with the DFT data given in Table 3. Note that the 1015 K/s results are closer to [8]. Simulation size Simulation type Cooling rate O(2)(%) O(3)(%) O(4)(%) O(5)(%) Zr(4)(%) Zr(5)(%) Zr(6)(%) Zr(7)(%) Zr(8)(%) Zr(9)(%)

1500 Atoms MD ReaxFF 2.57 × 1015 K/s 1015 K/s 1014 K/s 20.4 7.2 1.1 63.1 69.8 73.8 15.9 22.4 24.7 0.6 0.6 0.4 5.4 0.2 0 26.6 15 5.2 41.6 43.4 48.8 22.8 34.8 38.2 3 6.4 7.6 0.6 0.2 0.2

1013 K/s 0.6 74.2 24.9 0.3 0.2 4.2 48.6 39.8 7 0.2

96 Atoms DFT Ref.[10] DFT Ref.[8] 3.22 6.45 69.35 64.52 29.03 30.65 1.61 1.61 0 0 6.25 18.75 37.5 34.38 50 40.63 6.25 6.25 0 0

Table 3: Comparison of our obtained coordinations with other DFT simulations. The cut-off radius was chosen to be 3˚ A.

11

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

There is one difference between this work and both of the DFT results: there is a higher distribution of Zr(6) than Zr(7). To explain this, one needs to compare the coordination distribution in Figure (a) (right after the initial NVT annealing) with the the data given in Table 3 (after the final cooling). The ratios of Zr(6)/Zr(7) in Figure (a) are more similar to those predicted by DFT. We hypothesize that this may be related to the constant volume of the simulation cell in the DFT method. By allowing the simulation cell to change its volume in an NPT ensemble, we are able to observe the coordination distribution rearrangement (Figure (b)-(g)).

12

ACCEPTED MANUSCRIPT

PT

5. Conclusions and future works

AC CE P

TE

D

MA

NU

SC

RI

Classical Molecular Dynamics was used to develop supercells of a-ZrO2 . Structural analysis of samples heated at different temperatures and cooling rates were carried out. The heating temperature does not affect the final result as long as it is well beyond the melting point. Therefore, the annealing process can be performed with temperatures as low as 4000K (The computation time for annealing at low temperature is lower). A quantitative relation between the final density and cooling rate was established. This relation was further used to construct cells of lower densities. With our method, we were able to achieve densities very similar to experimental thin films fabricated by electron beam deposition. The supercells obtained by our method can further be used for quantum chemistry calculations. Based on the required accuracy and limitations in place one can further minimize the cell energy using DFT or just use it in its current form. In this work no further quantum relaxation of the supercells was performed. In terms of computation hardware and time requirements, each simulation in this work used a cluster of 36 processors with clock rates between 2.53.06 GHz and 12 × 2 GB memory for runs of 48 hours CPU time. This is much faster than typical quantum chemistry simulations (DFT calculations on systems of sizes in the range studied in this work are impractical due to computational complexity). In a separate work we analyze the hydrogen diffusion coefficients in a-ZrO2 with different densities and to investigate whether or not a relation between the density and diffusion coefficient can be established. This will help assess the reliability characteristics of a-ZrO2 based transistors. 6. Acknowledgement Our work made use of the infrastructure and resources of WestGrid (www.westgrid.ca) and Compute Canada Calcul Canada (www.computecanada.ca). The research is partly funded by Natural Sciences and Engineering Research Council of Canada (NSERC) and Cisco Systems, Inc. [1] ITRS, International Technology Roadmap for Semiconductors, Executive Summary (2011). URL http://www.itrs.net/Links/2011ITRS/2011Chapters/ 2011ExecSum.pdf 13

ACCEPTED MANUSCRIPT

RI

PT

[2] A. Salan, H. Grampeix, J. Buckley, C. Mannequin, C. Valle, P. Gonon, S. Jeannot, C. Gaumer, M. Gros-Jean, V. Jousseaume, Investigation of HfO2 and ZrO2 for Resistive Random Access Memory applications, Thin Solid Films 525 (2012) 20 – 27.

SC

[3] D. Panda, T.-Y. Tseng, Growth, dielectric properties, and memory device applications of ZrO2 thin films, Thin Solid Films 531 (2013) 1–20.

MA

NU

[4] G. Ye, H. Wang, S. Arulkumaran, G. I. Ng, R. Hofstetter, Y. Li, M. J. Anand, K. S. Ang, Y. K. T. Maung, S. C. Foo, Atomic layer deposition of zro2 as gate dielectrics for algan/gan metal-insulator-semiconductor high electron mobility transistors on silicon, Applied Physics Letters 103 (14) (2013) 142109.

D

[5] M. Youssef, B. Yildiz, Hydrogen defects in tetragonal ZrO2 studied using density functional theory., Physical chemistry chemical physics : PCCP 16 (4) (2014) 1354–65.

AC CE P

TE

[6] J. Livage, K. Doi, C. Mazires, Nature and Thermal Evolution of Amorphous Hvdrated Zirconium Oxide, Journal of the American Ceramic Society 51 (6) (1968) 349–353. [7] M. Winterer, Reverse Monte Carlo analysis of extended x-ray absorption fine structure spectra of monoclinic and amorphous zirconia, Journal of Applied Physics 88 (10) (2000) 5635–5644. [8] E. A. Chagarov, A. C. Kummel, Generation of Realistic Amorphous Al2O3 And ZrO2 Samples By Hybrid Classical and First-Principle Molecular Dynamics Simulations, ECS Transactions 16 (10) (2008) 773– 785. [9] D. Ceresoli, D. Vanderbilt, Structural and dielectric properties of amorphous ZrO2 and HfO2 , Phys. Rev. B 74 (2006) 125108. [10] X. Zhao, D. Ceresoli, D. Vanderbilt, Structural, electronic, and dielectric properties of amorphous ZrO2 from ab initio molecular dynamics, Phys. Rev. B 71 (2005) 085107. [11] A. C. T. Van Duin, S. Dasgupta, F. Lorant, W. A. Goddard, ReaxFF: A reactive force field for hydrocarbons, Journal of Physical Chemistry A 105 (41) (2001) 9396–9409. 14

ACCEPTED MANUSCRIPT

RI

[13] S. N. Laboratories, LAMMPS Website (2015). URL http://lammps.sandia.gov

PT

[12] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of Computational Physics 117 (1) (1995) 1 – 19.

NU

SC

[14] H. M. Aktulga, J. C. Fogarty, S. A. Pandit, A. Y. Grama, Parallel reactive molecular dynamics: Numerical methods and algorithmic techniques, Parallel Computing 38 (4-5) (2012) 245–259.

MA

[15] K. Vollmayr, W. Kob, K. Binder, Cooling-rate effects in amorphous silica: A computer-simulation study, Phys. Rev. B 54 (1996) 15808– 15827.

D

[16] J. C. Fogarty, H. M. Aktulga, A. Y. Grama, A. C. T. van Duin, S. A. Pandit, A reactive molecular dynamics simulation of the silica-water interface, The Journal of Chemical Physics 132 (17) (2010) 174704.

AC CE P

TE

[17] A. C. T. Van Duin, B. V. Merinov, S. S. Han, C. O. Dorso, W. A. Goddard, ReaxFF reactive force field for the Y-doped BaZrO3 proton conductor with applications to diffusion rates for multigranular systems, Journal of Physical Chemistry A 112 (45) (2008) 11414–11422. [18] C. J. Howard, R. J. Hill, B. E. Reichert, Structures of the ZrO2 Polymorphs at Room Temperature by High-Resolution Neutron Powder Diffraction, Acta Crystallographica Section B 44 (1988) 116–120. [19] R. W. G. Wyckoff, Crystal Structures, Vol. 1, Interscience Publishers, 1963. [20] G. Teufer, The crystal structure of tetragonal ZrO2 , Acta Crystallographica 15 (11) (1962) 1187. [21] E. H. Kisi, C. J. Howard, R. J. Hill, Crystal structure of orthorhombic zirconia in partially stabilized zirconia, Journal of the American Ceramic Society 72 (9) (1989) 1757–1760. [22] L.-G. Liu, New high pressure phases of ZrO2 and HfO2, Journal of Physics and Chemistry of Solids 41 (4) (1980) 331 – 334.

15

ACCEPTED MANUSCRIPT

PT

[23] G. Jomard, T. Petit, A. Pasturel, L. Magaud, G. Kresse, J. Hafner, Firstprinciples calculations to describe zirconia pseudopolymorphs, Physical Review B 59 (6) (1999) 4044–4052.

SC

RI

[24] S. Nose, A unified formulation of the constant temperature molecular dynamics methods, The Journal of Chemical Physics 81 (1) (1984) 511– 519.

NU

[25] W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Physical Review A 31 (3) (1985) 1695–1697.

AC CE P

TE

D

MA

[26] A. Stukowski, K. Albe, Extracting dislocations and non-dislocation crystal defects from atomistic simulation data, Modelling and Simulation in Materials Science and Engineering 18 (8) (2010) 085001.

16