Investigating Grain Boundary Structures and Energetics of Rutile with

Subscriber access provided by SUNY DOWNSTATE

Article

Investigating Grain Boundary Structures and Energetics of Rutile with Reactive Molecular Dynamics Patrick Jacob Shamberger, Jennifer L. Wohlwend, Ajit K Roy, and Andrey A Voevodin J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b02695 • Publication Date (Web): 03 May 2016 Downloaded from http://pubs.acs.org on May 4, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Investigating Grain Boundary Structures and Energetics of Rutile with Reactive Molecular Dynamics 1

Patrick J. Shamberger* Jennifer L. Wohlwend 2 Ajit K. Roy 4 Andrey A. Voevodin 2,3

1

Department of Materials Science and Engineering; Dwight Look College of Engineering; Texas A&M University; College Station, TX, 77843, USA 2

Nanoelectronic Materials Branch; Materials and Manufacturing Directorate; Air Force Research Laboratory; Wright Patterson AFB, OH 45459, USA 3

Universal Technology Corporation; Dayton, OH 45459, USA

3

Department of Materials Science and Engineering; University of North Texas; Denton, TX 76203, USA Corresponding Author: Patrick J. Shamberger e-mail: [email protected] phone: (979) 458-1086 Abstract Determining quantitative grain boundary (GB) energies as a function of microscopic orientation parameters is essential in order to understand the population of boundaries present in polycrystalline ceramics and films, and the physical properties that result from these boundaries. Here, we investigate the use of two reactive potentials, COMB3 and ReaxFF, to predict free surface and grain boundary structures and energies in the TiO2 rutile system, and compare these results against previously reported ab initio surface and interfacial energies. We demonstrate reactive MD potentials to be generally capable of reproducing key features anticipated for GB structures and energetics, including relative GB and surface energy, charge distributions and potential for different polar and non-polar terminations, and energy cusps at low-energy interfaces (e.g., coherent twin boundaries, coherent site lattice boundaries). This work establishes the foundation for further use of reactive MD to simulate libraries of oxide GBs and dynamic processes occurring along those GBs.

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile 1 Introduction Grain boundary interfacial energies, γGB, vary as a function of the macroscopic and microscopic parameters which describe the orientation and structure of that interface. Generally, this dependence is broken down to five macroscopic parameters (three degrees of freedom describing the crystallographic rotation between neighboring grains; two degrees of freedom describing the interface normal) and three microscopic parameters (describing translational displacement of one crystal with respect to the other).1-2 Determining quantitative γGB as a function of these orientation parameters is essential in order to understand the population of grain boundaries (GBs) present in polycrystalline ceramics and films. The resulting GB population plays a critical role in dictating mechanical properties (ductility, toughness), dielectric properties (charge trapping, leakage current), and ionic diffusion rates in polycrystalline materials. Interest in oxide interfaces and surfaces has existed for a number of decades, but has recently gained increased awareness due to the discovery of unexpected phenomena, such as two-dimensional sheets of conductive charges, existing at polar oxide interfaces.3-5 Despite their importance, γGB of oxide systems as a function of all eight microscopic parameters remain relatively unexplored for most oxide systems. Experimental approaches to quantify γGB are often challenging, and are limited by the kinetics of low-temperature equilibration.6-9 Thus, validated simulation techniques to investigate the structure and energetics of oxide interfaces are of critical interest to the community. First-principles density functional theory (DFT) has been utilized extensively to predict energies and configurations of oxide surfaces and interfaces.10-12 However, DFT is extremely limited in the number of atoms that can be simulated simultaneously, thereby limiting which particular interface orientations are accessible. Furthermore, due to the computational costs, most DFT studies are limited to only a few specific boundaries, and are therefore not appropriate for developing an overall distribution of the energetics of a specific oxide system. Classical atomistic models with empirical potentials have demonstrated utility in predicting γGB in a number of metallic systems, and, in contrast with DFT, allow for more rapid simulation of systems with greater numbers of atoms.13-15 However, oxide surfaces and interfaces are inherently more complicated than their metal counterparts, due to the potential existence of localized changes in electron density distribution that are expected to occur near interfaces due to over- or under-coordination of metal cations,10, 16 which correspond with a 2 ACS Paragon Plus Environment

Page 3 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile change in the degree of bond covalency near the interface.10, 17 Empirical potentials can be broadly subdivided into reactive and non-reactive potentials. Nonreactive potentials are generally held to be reliable only relatively near equilibrium; reactive potentials are anticipated to better recreate interface energetics, as they can account for over- and under-coordinated configurations, as well as bond breaking and reforming using variable bond-order formalisms, both of which are likely to occur at free surfaces and interfaces.18-19 Furthermore, many reactive potentials incorporate charge re-equilibration, which is necessary to more accurately describe mixed valence states expected to occur near interfaces in metal oxides.20-22 An excellent review discussing the details of reactive potentials has been published recently.18 Here, we investigate the use of reactive empirical potentials in molecular dynamics (MD) simulations to investigate the structure and energy of surfaces and GBs in oxide systems. This validation study is a critical component of utilizing MD to investigate dynamic processes occurring along GBs. We focus on one well-investigated model oxide system, TiO2 (rutile), and utilize two different well-documented reactive MD potentials for investigating grain boundary structure and energetics: the third generation charge optimized many body (COMB3) potential,23 and the reactive force field (ReaxFF) potential.24-25 These two methods utilize slightly different approaches in calculating variable bond-order, as well as different approaches to training an empirical potential set; for more detailed discussion, see 18. Thus, comparison of results derived from the two potentials adds greater insight into the use of reactive MD techniques to simulate GBs, in contrast to results which may be specific to a given potential. We first investigate different free surfaces, including low Miller-index and select members of {h0l} and {hk0} surfaces, to validate the use of empirical reactive MD potentials. We report structures, surface energies (γsurf), as well as charge distributions and internal electrostatic potentials, for different surface orientations and terminations from the two potentials of interest and compare these results against previously reported results. We next investigate a subset of symmetric tilt GBs in the rutile system, where neighboring grains are related by either a mirror reflection across the GB, or a two-fold rotation around an axis lying in the plane of the GB, including the effect of rigid-body translation in the plane of the GB. Again, we compare the output of MD simulations against previous simulation results for GBs, principally those obtained by DFT. We report grain boundary energy, γGB, as a function of initial rigid-body translation parallel to the grain boundary, and discuss the accessibility of different local energy minima. We 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile particularly discuss results for low-energy GBs consisting of rotation and reflection twins or pseudo-twins. We conclude by investigating γGB in symmetric tilt GBs across a wide range of rotation angles between original surface orientations, in order to demonstrate the utility of this method in rapidly scanning γGB space.

2 Methods 2.1 Simulation System An orthorhombic computational cell (Fig. 1) is defined with three-dimensional periodic boundary conditions. In the case of a single crystal (reference state), the crystal lattice is continuous throughout the cell. For free surface simulations, only half of the cell is occupied by atoms while the remaining volume is a vacuum. Thus, two free surfaces exist in the system, on opposite sides of the slab; these surfaces may be equivalent or, in the case of polar terminations (outermost layers with net electrostatic dipoles), they may be dissimilar. In this latter case, reported surface energies represent the average of the surface energies of the two free surfaces. For grain boundary simulations, two equi-dimensional crystals with different orientations exist in the cell, separated by a GB. Thus, two GBs exist in the system – one in the middle of the cell, and one on the upper and lower bounds of the cell (Fig. 1). Again, these GBs are not strictly equivalent in all cases, in which case reported γGB represents the average of the two. The crystal lattice in the upper and lower crystals are related to each other by some crystallographic transformation; in this study, that transformation may take the form of either a reflection across the GB (Fig. 1), or a two-fold rotation around an axis that lies in the plane of the GB, followed by a rigid-body translation of one grain with respect to another. The dimensions of the computational cell are calculated to meet periodic geometric constraints imposed by the tetragonal crystal symmetry of rutile and by the selected interfacial plane (see supplementary information). The surface area of each slab is given by = ∙ ∙ , where x and z are the minimum periodic repeating length, and , are integers which scale the overall surface area. Initial slab areas are between 125 and 425 Å2. Repeated calculations on different super-cells with slab areas up to 3000 Å2 perturbed γsurf by less than 2 %, demonstrating that initial slab areas were sufficiently large. Slabs were defined with interface separation distances, = ∙ , where y is the minimum periodic repeating length, and is some integer. Thicknesses of slabs, , were always greater than 2.5 nm, 4 ACS Paragon Plus Environment

Page 5 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile which is a comparable distance to that used in previous studies.13-15 Repeated calculations on different super-cells with thicknesses between 2.5 nm and 5.0 nm perturbed γsurf by less than 2 %, demonstrating that initial slab areas were sufficiently thick to avoid interaction between neighboring interfaces. Given these computational volumes, GB simulation cells contained between 972 and 2880 atoms; free surface cells contained half this number. Computational cells always contained a stoichiometric 1:2 ratio of titanium and oxygen atoms. Reported γsurf are calculated as the area-normalized difference between the potential energy for the system with two free surfaces (Eslab) and the potential energy of the single crystal with the equal number of atoms · relaxed from experimentally determined initial atomic positions using the same process: = − · ⁄2 ,

(1)

The surface area, 2Asurf, accounts for both top and bottom free surfaces in the computational cell. Similarly, GB energies are calculated as the area-normalized difference between the potential energy for the entire bicrystal system with two GBs and of the single crystal:

!

= "#$ − · %⁄2

!.

(2)

Systems are relaxed to some metastable configuration, about which energy oscillates slightly due to thermal fluctuations. Reported γsurf or γGB for a particular interface are calculated from minimum potential energies observed after relaxation. The principle source of numerical uncertainty in eq. 1 and 2 derives from uncertainty on the estimate of ETiO2; (±1.5 meV/TiO2) this value is estimated by repeated relaxation of single crystals of different sizes. Reported uncertainty in tables and figures is propagated from uncertainty for ETiO2. 2.2 Surface Orientations Relaxed atomic configurations and γsurf for a particular free surface depend on both the crystallographic indices of the investigated plane, as well as the particular termination selected. This work includes: 1) investigation of the four principle low Miller-index crystallographic planes in the rutile structure – (100), (001), (110), (101), and select representatives of 2) {h0l} families of planes, and 3) {hk0} families of planes. For low Miller-index crystallographic planes, all possible planar terminations (i.e., those which can be generated simply by removing all atoms that fall on one side of a plane that defines a free surface) are evaluated. For both {h0l} and {hk0} families, all non-polar terminations, as well as all polar planar terminations are 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile evaluated. In addition, those surface reconstructions which are composed of terraced low Millerindex surfaces are always considered. These reconstructions always represent non-polar terminations, but are not necessarily planar, as defined above.

2.3 Grain Boundary Orientations We investigated symmetric tilt GBs based on a subset of those free surfaces investigated in the previous section. GB orientations are related by either a reflection or a 2-fold rotation across {h0l} and {hk0} families of planes with the different terminations considered above, coupled with a rigid-body lattice translation of one crystal (referred to here as refl+trans or rot+trans, respectively). In the general case, rigid-body translations are required to avoid relatively high-energy conformations where atoms exist in close proximity to each other on either side of the crystallographic interface; this situation is even more significant in the case of ionic materials, where these proximal atoms may have like charge. In certain cases, refl+trans or rot+trans may result in coherent rotation or reflection twins, or pseudo-twins (where the cation sub-lattice shares a certain symmetry across the GB, but the anion lattice does not). In order to test these different configurations, one lattice is displaced by rigid-body translations relative to the other in two orthogonal directions parallel to the GB (& , & ; Fig. 1). This has the effect of sampling different initial GB configurations in an effort to find low-energy configurations, and to measure the accessibility of these configurations (i.e., to determine their occupancy after some period of relaxation). An integral number ( , ) of evenly spaced lattice displacements given by ⁄ (where x is again the minimum periodic repeat length) and ⁄ form a grid of initial starting configurations. We typically adopted values of nx and nz which resulted in grid spacing of ~50 pm. This resolution was observed to convincingly localize regions of low-energy and high-energy configurations. Without taking any additional precautions, atoms located near the GB could potentially be positioned very near to atoms in the opposing crystal, thereby causing highly unfavorable starting configurations which can result in an exceptionally high potential energy that terminates the MD simulation. Thus, in addition to displacements in the GB plane, an artificial vacuum separation gap (dgap = 25 pm) was created at each GB. During simulations, the overall box volume was allowed to relax to an average zero pressure state, during which process the gap was closed and interfacial atoms adopted relaxed positions. 6 ACS Paragon Plus Environment

Page 7 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile

2.4 Simulation Approach Surface and grain boundary structures and energies are calculated from computational cells relaxed at low temperatures (1 K) using molecular dynamics, as implemented in LAMMPS.26-28 Low-temperature relaxation was implemented in this study for three reasons: 1) this approach avoided the tendency to get trapped in relatively high-energy local minima due to the energetic topology of the surface (observed in simple minimization-based approaches), 2) this process allowed for surface relaxation without significant bulk conformational changes, permitting calculation of γGB from the bulk bicrystal, and 3) low-temperature relaxation retained information about interfacial energy, γGB, as a function of initial in-plane lattice displacement, which is valuable for understanding energetic costs of displacing one crystal with respect to the other. In contrast, relaxing GBs at higher temperatures tended to result in systems finding more stable overall configurations, thereby losing information about those initial higher-energy displacement configurations. Low-temperature relaxation was implemented by running dynamic simulations at 1 K under isothermal-isobaric (NPT) conditions with short timesteps (0.1 to 0.05 fs) and implementing charge equilibration using the QEq method.20 Initial random atomic velocities are assigned for a 1 K ensemble.

2.5 Atomic Potentials Due to extensive interest in TiO2 (rutile) over the past decades, a number of empirical potentials have been defined to describe this oxide system. These include multiple fixed-charge models,29-31 as well as variable charge models based on the QEq charge equilibration scheme.12, 32-34

Here, we investigate the use of two reactive MD potentials, COMB3 35 and ReaxFF.36-37

TiO2 (rutile) was optimized using both potentials producing structural parameters consistent with higher fidelity calculations, and with the previously reported parameters using these potentials. 35-37

COMB3 has been used previously to describe a range of systems including polymers,23 metals,28, 38 and metal oxides;39-42 they consist of electrostatic, short-range, as well as van der Waal interaction terms. COMB3 relies on the electronegativity (QEq) method for charge equilibration.20 Bond-order calculations are coupled with dynamic charges, which are suggested to lead to more accurate calculations of charge state on ions.18 A potential set for Ti/TiO2/Cu 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile was developed to investigate redox reactions at the interface of Cu and TiO2.35 This parameter set was trained against first-principles data for 10 different AB2 phases. A weighted leastsquares method was employed to fit potential parameters based on predicted lattice parameters, cohesive energy, elastic constants, vacancy and interstitial defects; weighting factors favored rutile (TiO2) and hcp (Ti) physical properties, as these were the two phases of interest for this study. After parameterization, predictions of stacking fault energies (Ti) and surface energy formations and charge states (TiO2) are validated against values predicted from first principles. The Ti/TiO2 system has recently been re-parameterized for a new release of COMB3 following the completion of this study (private communications); this modification is expected to change the absolute values of surface and interface energies, but not the general behavior of simulations. ReaxFF was originally derived for hydrocarbons,24 but later expanded to include a number of systems relevant to combustion processes,43 transition metal oxidation,44-45 and aqueous systems,46-47 many of which included metal/metal-oxides. ReaxFF also relies on the electronegativity (QEq) method for charge equilibration.20 ReaxFF retains significant bondorder information in intermediate states, which is suggested to lead to greater fidelity of dynamic processes, and intermediate (high energy) states. A Ti/TiO2/H2O parameter set was developed to investigate water dissociation on both rutile and anatase surfaces.37 This parameter set was trained against first-principles data for 8 different AB2 phases, as well as 150 bare and hydrated TiO2 cluster structures. A weighted least-squares method was employed to fit potential parameters based on the equation of state, surface formation, water-binding energy, and structure and energy of TiO2 clusters. This original parameter set was subsequently modified;36 both versions are considered here.

3 Results 3.1 Free Surface Energies 3.1.1 Comparison with Previous Results Low Miller-index planes in rutile have been investigated previously using ab initio techniques, based on either linear combination of atomic orbitals (LCAO) techniques12, 48-53 or plane-wave (PW) approximations,54-62 classical atomistic simulations based on fixed charges,2930, 32, 63

and semi-empirical density functional-based tight binding (DFTB) approximations (Table

1).12, 33, 48 Ab initio approaches uniformly predict γ(110) < γ(100) < γ(001); γ(110) is predicted to fall 8 ACS Paragon Plus Environment

Page 9 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile within the range (0.40 to 1.14) J·m-2, γ(100) within the range (0.68 to 1.20) J·m-2, and γ(001) within the range (0.83 to 2.20) J·m-2. γ(101) is far less frequently considered; in those cases where it has been investigated, γ(110) < γ(101). However, the exact sequence of all four surface energies is unclear from previous studies. Both COMB3 and ReaxFF MD potentials result in low-index γsurf that are generally consistent with ab initio results. Both ReaxFF potential parameter sets predict the same energy sequence as calculated using ab initio techniques, γ(110) < γ(100) < γ(001). COMB3 potential predicts a slightly different ordered sequence, γ(110) < γ(001) < γ(100), however γ(001) and γ(100) are within 10% of each other (Table 1). Recent re-parameterization has modified surface energies to predict γ(110) < γ(100) < γ(001), as observed using ab initio techniques (private communications). We calculate γ(110) < γ(101) using both the COMB3 and the ReaxFF potentials, consistent with previous results. While predicted surface energies generally fall within the range observed using ab initio techniques, γ(100) calculated using COMB3, and γ(001) calculated using one of the two available force fields in the ReaxFF scheme exceed those ranges. Furthermore, the energy difference between γ(110) and γ(100) from ab initio studies is generally between (0.2 to 0.3) J·m-2; for COMB3, this difference is 0.36 J·m-2 while for ReaxFF, this difference is 0.18 or 0.11 J·m-2. The energy difference between γ(110) and γ(001) from ab initio studies is generally between (0.8 to 1.2) J·m-2; for COMB3, this difference is 0.25 J·m-2 while for ReaxFF, this difference is 1.86 or 1.36 J·m-2. Surface energies calculated here using COMB3 (γ(110) < γ(001) < γ(100)) have a different energy sequence than originally reported (γ(110) < γ(100) < γ(001)) .35 We were able to reproduce the original result by comparing surface energies before the system was fully relaxed. However, complete relaxation of free surfaces results in the energy sequence reported here. While the energy sequence of low-index planes matches that in those planes calculated for rutile in the original ReaxFF publication, our surface energies are nearly double those reported in that study. While the origin of this discrepancy is unclear, it may also be due to differences in surface relaxation between the two groups. 3.1.2 Surface Terminations: non-polar and polar surfaces For (100), (101), and (110) surfaces, two distinct planar terminations are available: a nonpolar termination (type II interface, following the Tasker classification),64-65 and a polar termination (type III interface; Fig. 2). For the (001) surface, the only possibility is a type I non9 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile polar termination. In all cases, γsurf of non-polar interfaces is less than γsurf of their polar counterparts (Table 2); this energy penalty ranges from (0.9 to 1.9) J·m-2 for the COMB3 results, and from (1.1 to 2.1) J·m-2 for the ReaxFF results, and represents the energetic cost of internal charge compensation necessary to keep the potential from diverging. Non-polar terminations lead to a symmetric distribution of charges, while polar terminations are compensated by partial charge-transfer from one side of the slab to the other (Fig. 3). Non-polar terminations are always O-terminated; these outer-most planes of O anions have less charge than interior O anions, due to the lack of charge transfer from the absent neighboring Ti cations (which would be present if the crystal was continuous). For the case of polar terminations, an additional plane of oxygen anions is essentially transferred from one side of the slab (d = 0) to the other side (d = dmax). This outermost oxygen layer is shielded from any nearby Ti cations by another plane of O-anions; thus, the outermost oxygen layer has very small charges on these atoms. This positive charge (or lack of negative charge) is compensated by a less positive charge on Ti-cations and a larger negative charge on O-anions on the other side of the slab. In general, the absolute magnitude of charge on atoms predicted using COMB3 is larger than those predicted by ReaxFF. +_

The electrostatic potential (' = − (-

* , where l is the distance into the slab) for

non-polar terminations are symmetric due to the symmetry of charge distribution in the slab, whereas the potential tends to diverge to a limited degree for slabs with polar terminations (Fig. 4). Slabs with non-polar terminations tend to have a small non-zero deflection due to surface relaxation effects. In contrast, polar terminations result in non-zero potential on one side of the slab, due to the existence of the polarized surface. We note that complete charge transfer to compensate for the polar surface would result in no excess potential. Thus, both COMB3 and ReaxFF predict partial compensation, in which the energetics of charge transfer counteract the electrostatic energy that results from a polar termination on a slab of a finite thickness. 3.1.3 {h0l} Surfaces We calculated surface structures and energies for a range of {h01} and {10l} surfaces, comparing simple planar terminations (those obtained by removing all atoms on one side of a flat plane) against terraced terminations (those consisting of combinations of (101), (100) and (001) surfaces; Fig. 5). All unique planar terminations are considered. In the cases of the (104), (102), 10 ACS Paragon Plus Environment

Page 11 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile and (401) surfaces, it is not possible to obtain terraced surfaces from simple planar terminations. That is, due to geometric constraints of atomic positions within the unit cell and the orientation of the surface plane, some ionic reconstruction is required to obtain terraced surfaces (illustrated in Fig. 5). Planar terminations for (104), (102), and (401) are necessarily type-III polar surfaces, as the position of the outermost oxygen anions change slightly between the two termination cases, although the resulting dipole is not large compared to other polar terminations considered in this work. In all other {h0l} surfaces considered here, it is possible to adjust the position of a planar termination to result in terraced surfaces. All terraced surfaces are type-II non-polar surfaces, according to the Tasker classification (Table 3). Surface energies of terraced terminations are uniformly lower than those of planar terminations by (0.5 to 1.0) J·m-2 (Table 3). This is likely due to the existence of unshielded Ti cations and under-coordinated O anions occurring along the planar termination. The COMB3 and ReaxFF potentials predict significantly different values and trends for γsurf for {h0l}. In general, γsurf (COMB3) / γsurf (ReaxFF) is between 0.5 and 2.0; this mirrors the difference in the low-index planes, where γsurf (COMB3) > γsurf (ReaxFF) for the (100) surface, but this relationship is reversed for the (101) surface. Surface energies for the two different parameterized ReaxFF sets are within 25 % of each other. 3.1.4 {hk0} Surfaces We additionally calculated surface structures and energies for {h10} surfaces, comparing simple planar terminations against terraced terminations, as above (Fig. 6). In the case of {h10}, terraced terminations consist of different ratios of (100) and (110) type surfaces (Fig. 6). All unique simple planar terminations (both polar and non-polar) are considered. It is not possible to generate terraced surfaces for any of the {h10} surfaces investigated using a simple planar termination; thus, these surfaces require ionic reconstruction. All terraced surfaces are type-II non-polar surfaces, while planar terminations may be either type-II non-polar, or type-III polar. Surface energies of terraced terminations are uniformly lower than that of planar terminations by (0.6 to 1.2) J·m-2 (Table 4). In comparison, there is much less difference between surface energies in either polar or non-polar planar terminations. Thus, the energy stabilization observed in terraced terminations is due to their composition of stable low-index planes, which minimizes the extent of ‘dangling bonds’ at the free surface. COMB3 predicts significantly larger γsurf for {h0l} than ReaxFF. In general, γsurf(COMB3) / γsurf ReaxFF) is 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile between 1.5 and 1.9; this is consistent with the observation γsurf(COMB3) > γsurf(ReaxFF) for both the (100) and (110) surfaces. Surface energies for the two different parameterized ReaxFF sets are within 25 % of each other. 3.1.5 Atomic Structure of Stable Surfaces Following fundamental surface energy minimization arguments, surfaces of unstable (higher energy) high Miller-index planes are predicted to decompose into atomic-scale terraces consisting of stable (lower energy) planes to minimize energy. The resulting surface energy of these terraced surfaces, γ′, is then given by the weighted average of the surface energies of the constituent terrace surfaces:66-67 ′ =

/ 0

∑ 2 .

(3)

In general, the most stable atomic configuration is composed of terraces of low energy surfaces (Fig. 7, Table 3, 4). For the COMB3 potential, γsurf of both {h0l} and {hk0} decrease monotonically with rotation of the surface plane orientation from γ(100) to γ(101), from γ(001) to γ(101), and from γ(100) to γ(010). γsurf of terraced surfaces is generally consistent with calculated γ′. However, along the transect γ(100) to γ(101), γsurf predicted by molecular dynamics simulations is significantly less than γ′ (Fig. 7). This suggests that there may be additional surface relaxation along this interface, presumably at the junctions of two different surface types, which is not considered in Eq. 3. Similar trends are observed for the ReaxFF potential (Fig. 7). γsurf from Ti/O/H parameterization from Kim et al. 36 are illustrated here; results using parameterization from Raju et al. 37 tended to be more susceptible to minor surface reconstruction. As before, along the transect γ(100) to γ(101), γsurf predicted by molecular dynamics simulations is significantly less than γ′. Differences between the two potentials are largely limited to differences in the relative energies of the low-index surfaces.

3.2 Grain Boundary Energies 3.2.1 Comparison with Previous Results Due to their natural abundance, (101) and (301) coherent twins in rutile have been the subject of previous study by high-resolution TEM;68-69 furthermore, γ(101)GB and γ(301)GB have been calculated using non-reactive MD potentials.63 The lowest energy configuration for (101) twins consists of a pseudo-reflection twin, a reflection twin coupled with a ½ [111] translation in 12 ACS Paragon Plus Environment

Page 13 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile the plane of the GB, resulting in reflection symmetry across the twin boundary only for the cation sub-lattice;63 in rutile, the (101) pseudo-reflection twin is equivalent to a rotation twin around the [010] rotation axis. For (301), the lowest energy configuration is believed to be a reflection twin;63 in rutile, the (301) reflection twin is equivalent to a pseudo-rotation twin around the [010] rotation axis (2-fold rotation, with a ½ [113] centering operation). In this study, both γ(301)GB and γ(101)GB are notable as having significantly lower interfacial energies than other {h01} GBs (Table 5, 6). For these two members of {h01}, the GB interface is significantly lower in energy, with respect to the free interface, as shown by low values of γGB/γsurf. For COMB3, γ(101)GB/γ(101)surf = 0.07, γ(301)GB/γ(301)surf = 0.24, while for ReaxFF, γ(101)GB/γ(101)surf = 0.20, γ(301)GB/γ(301)surf = 0.13. Furthermore, we confirm the previously identified lowest energy configurations and recover that anticipated equivalence between reflection and rotation type twins with additional centering translation. However, when comparing the relative energies of these two GBs, the results of the two potentials are inconsistent; for COMB3, γ(301)GB > γ(101)GB , while for ReaxFF, γ(101)GB > γ(301)GB . In addition, γGB has been determined using ab initio techniques for coincident site lattice (CSL) members of the {hk0} family of planes, including the Σ1 (100) reflection twin, and Σ5 (210) and Σ5 (310) symmetric tilt GBs.10-11, 63 Using both potentials, we observed γ(310)GB > γ(210)GB > γ(100)GB, which is consistent with the higher fidelity, ab initio results. γ(100), which represents a reflection twin, is on the order of the energy of the (101) and (103) twin GBs. However, ReaxFF significantly underestimates γGB in comparison to ab initio results; γGB(ReaxFF)/ γGB(ab initio) ~ 0.2 to 0.3, whereas γGB(COMB3)/ γGB(ab initio) ~ 0.6 to 1.1. 3.2.2 γGB Topography In general, incoherent twins do not have low γGB, due to over- or under-coordinated ions and changes in bond length across the GB. By introducing three additional translational degrees of freedom, consisting of in-plane displacement, as well as relaxation orthogonal to the interface, lower γGB configurations are generally observed. We have investigated γGB of select {hkl} GBs, as a function of initial displacement (Table 5, 6, 7). Systems were relaxed at 1 K for up to 5 ps, although the majority of the relaxation occurs in the initial 50 fs (see supplementary information). In general, two types of GBs are observed: 1) low-energy GBs, which are marked by a small minimum γGB, with relatively high accessibility, nearly coherent GB structures (little to no perturbation in bond lengths and bond angles), as well as 2) high-energy GBs, which have 13 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile a relatively large γGB, a relatively poor accessibility of this minimum γGB configuration, and significantly disrupted GB structures (frequent over/under coordination, perturbation in bond lenths). These two classes of GBs are observed using both COMB3 as well as ReaxFF potentials. As an example of low-energy GBs, the relaxed γGB are illustrated as a function of initial displacement for the (301) GB refl+trans (Fig. 8 a,b) and rot+trans (Fig. 8 c,d); these energy topography maps are accompanied by histograms showing the relative abundance of different relaxed γGB. In both cases, a relatively large range of starting displacements relax to a nearly identical final structure, indicating high accessibility of this configuration.14 The similarity of the two topography maps (which are related by a ½ [113] translation) is caused by the underlying symmetry of the rutile structure, and illustrates the relationship between the (301) reflection and pseudo-rotation twins. In fact, it may be demonstrated that for (h01) GB interfaces, when h is odd, the pseudo-rotation and pseudo-reflection twins are related by a face-centering transformation along the interfacial plane (½ [11h]). In contrast, a high-energy GB still has regions of relatively high and relatively low energy, however the minimum γGB is significantly larger, and the accessibility of these low-energy configurations is relatively low (Fig. 9). COMB3 and ReaxFF generally predict similar, although not identical γGB topographies, suggesting that the underlying distributions and symmetries are inherent to the GB, rather than to a specific potential. As an example, γGB of a (105) pseudo-reflection twin have very similar distributions for a [501] displacement less than ~1.5 nm; whereas relaxed γGB diverge at [501] displacements greater than ~1.5 nm (Fig. 10). It is important to note that this observation does not necessarily imply that the global minimum-energy relaxed structure of this GB is different in the two cases, rather that different initial displacements are not necessarily all able to relax to the global minimum-energy configuration. Thus, while the general energy features in γGB topography tend to be independent of empirical potential, detailed maps relying on a single potential should be used with caution. 3.2.3 {h0l} Grain Boundaries Use of reactive MD potentials allows for a rapid, systematic investigation of γGB as a function of orientation parameters that can provide additional fidelity over simplified analytical relationships. As an example, the γGB of {h0l} rot+trans tilt GBs, which are related by rotating the interfacial plane around the [010] axis, are illustrated in Fig. 11. It is clear from this 14 ACS Paragon Plus Environment

Page 15 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile distribution that special low-energy interfaces are clearly identified by both MD potentials; these show up, as expected, as cusps in the γGB as a function of the angle between [001] in neighboring grains. Due to the symmetry of the rutile structure, a reflection across a (001) plane leads to a continuous, undisturbed lattice; thus, γ(001)GB = 0. In contrast, (100) leads to a continuous cation lattice (and is therefore designated a Σ1 CSL, on a cation-basis), but contains a reflection symmetry in the anion lattice; thus, γ100)GB is non-zero. Other low-energy configurations (301) and (101) are all captured using reactive MD, suggesting that this approach is able to identify relative energy cusps that are not necessarily known a priori. Comparing the two empirical potentials, COMB3 tends to predict deeper energy cusps than ReaxFF.

4 Discussion 4.1 Validity Because direct experimental measurement of the surface energy of different free surfaces and different GBs is not readily available for the rutile structure, validity of reactive MD results is best determined by direct comparison with ab initio results, which are generally assumed to have higher fidelity. Considering previously published studies using somewhat different techniques (Table 1), it is possible to bracket potential energy ranges for γsurf, and to identify a consistent ordered sequence of γsurf. In this context, while both reactive potentials generally reproduce the sequence of free boundary energies predicted by ab initio techniques, COMB3 predicts γ(001) < γ(100), which is opposite the ab initio results. Furthermore, γ(100) calculated using COMB3 potentials falls outside the range observed by ab initio techniques. In comparison, both ReaxFF potential parameter sets studied meet the predicted sequence and fall within the expected range of γsurf, with the exception that one of the parameter sets predicted a value of γ(001) greater than that obtained from ab initio methods.37 For this reason, we focused on the other available ReaxFF parameter set in the remainder of this study.36 Ab initio predictions of select γGB are also available to validate MD predictions. In this case, while both COMB3 and ReaxFF predict the observed sequence: γ(100) < γ(210) < γ(310), ReaxFF γGB are uniformly low (only 20 to 30 % of the predicted ab initio values) with respect to COMB3 values (50 to 110 % of the predicted ab initio values). Thus, while MD results predict γsurf and γGB which are generally consistent with ab initio results, caution should be exercised in over-interpreting quantitative results; as an example, quantitative results from ReaxFF would 15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile overestimate the abundance of low-energy planes, which could alter expected collective properties of a bulk polycrystalline ceramic. Several other general characteristics of GB simulations using reactive MD potentials support the general validity of this technique. First, it is clearly observed that atomistic results converge with general thermodynamic theory in predicting lowest-energy decomposition of free surfaces into terraced low-index stable planes (Fig. 7). Second, without foreknowledge, atomistic simulations are able to recover “special” GB orientations that are either predicted through near-CSL theory, or are generally observed in natural twinning of rutile crystals (Fig. 11). Finally, considering different planar terminations, it is generally observed that the non-polar terminations are, as anticipated, the lowest energy terminations. In contrast, polar terminations lead to divergent internal potential (Fig. 4) due to the incomplete charge transfer from one interface to another, which increases the overall energy of these systems with respect to nonpolar terminations. The fidelity of predicted charge distributions with respect to ab initio predictions would be a worthwhile focus of a more detailed investigation

4.2 Comparison between Reactive Potentials There are a number of distinct differences between surface and GB energies predicted by the empirical potentials investigated in this study. In the case of free surfaces, COMB3 (γ(110) < γ(101) < γ(001) < γ(100)) predicts a different energy sequence than ReaxFF (γ(110) < γ(100) < γ(101) < γ(001)). These low-index free surfaces represent stable crystals with one or two unique bond types severed at the interface. Thus, the predicted energy sequence is a direct consequence of the energy stored in those bonds (less the stabilization energy that would result from their absence). Therefore, the distinction between these results likely results primarily from differences in bond parameterization used in both methods. It should be noted that neither of these parameter sets were trained exclusively for surface energies; GB’s were not used at all. Thus, it could be possible in both cases to improve the agreement with ab initio γsurf by explicitly including γsurf in the training set with a significant weighting. In contrast, the energies of GBs also depend strongly on differences in bond order calculation, which are approached differently in the two parameter sets. This results in a significant difference in the depth of the energy cusps observed at low-energy GBs; COMB3 predicts much deeper energy cusps than ReaxFF (Fig. 11). It is notable that γGB tends to be a 16 ACS Paragon Plus Environment

Page 17 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile factor of 2 greater in high angle GBs calculated using COMB3. While γGB for a general high angle GB is typically difficult to determine in DFT due to the large computational cell needed, it is worth noting that ab initio results for Σ5 near-CSL (210) GBs are 1.7 to 1.9 J·m-2, while Σ5 near-CSL (310) GBs are 2.4 to 2.8 J·m-2. These special (presumably lower-energy) cases are still predicted to have significantly greater γGB than all of the GBs calculated using ReaxFF (but are comparable with GB calculated using COMB3). Thus, it is conceivable that ReaxFF overestimates bond-strength, which has the net result of generally stabilizing GBs.

4.3 General Utility for GB Investigations The general challenge for simulating GBs laid out at the beginning of this work, was the accurate prediction of γGB as a function of 8 microscopic parameters. We have demonstrated the utility of this technique in terms of rapidly scanning of this multi-variant parameter space. It is worth emphasizing that this scale of investigation would not be feasible relying on ab initio techniques; in fact, even calculation of a single high-miller index GB, which involves a large number of atoms in the unit cell, is challenging for ab initio techniques. In terms of the five macroscopic parameters relating neighboring grain rotation and the orientation of the interface, MD techniques have been previously demonstrated in metallic systems.13-15 Here, we demonstrate this technique in a well-understood model oxide system, capturing the salient features of the γGB relationships, as we understand them. In the absence of more detailed information, γGB may be estimated as:

!

= ,/ + , − ! ,

(4)

where EB is a term that describes the bonding (stabilization) energy across the GB.7, 9 As a simplifying assumption, EB may be considered to be approximately constant for high-angle GBs (Fig. 11). Importantly, γGB as predicted by reactive MD adds additional fidelity over the simple analytical expression given by (4). Namely, the approach adopted in this study 1) identifies GB orientations that actually represent a continuous lattice (e.g., γ(001)GB = 0), 2) identifies lowenergy cusps rather than always predicting smooth changes in γGB between stable low-index surfaces (e.g., minima at γ(301)GB), and 3) if based on well-validated reactive MD potentials, can account for differences in EB across different low-index GBs (e.g., γ(100)GB, γ(101)GB). Additionally, in oxide systems, γGB is highly dependent on three microscopic parameters describing the translational displacement of one crystal with respect to the other; the general 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile coherence between the two empirical potentials used in this study supports the overall validity of this technique to rapidly scan over in-plane displacements. This work represents a foundation for further exploitation of MD techniques to investigate different aspects of GB structure and behavior that are very appropriate for this simulation length and time-scale. Specifically, we envision utilizing the described MD technique in a number of different fashions: 1) Pre-screening over a large range of initial GB configurations. Given the rapid nature of this approach relative to other techniques, MD could be used to identify different lowenergy configurations which could then be utilized as initial configurations for more quantitative, computationally intensive approaches. For example, we envision a hierarchical approach where the outputs of relaxed MD structures could be adopted as starting structures for ab initio techniques, thereby minimizing the overall computational cost of a particular GB relaxation. 2) Identification of unique low-energy GBs. In the general case, special low-energy GB configurations are not necessarily known ahead of time. In the case of low-symmetry structures, it is not trivial to identify all of these configurations, and to define a subset which are significantly stabilized with respect to a general high-index GB. Again, the ability to test a large number of starting configurations enables the rapid identification of these relatively important GB configurations. 3) Dynamic investigations of transient processes. The nature of molecular dynamics is well-suited to investigate evolution of GB structures over time. For example, this technique could be used to investigate slip along GBs due to shear, or coarsening of grains through motion of triple-junctions. 4) Transport processes along GBs. Mass transport by atomic diffusion along GBs has traditionally been investigated in one of two ways: either by using nudged-elastic-band techniques to investigate the energetic barriers that limit hopping of an individual atom along some particular path, or by random-walk simulations, which track the motion of individual particles at some temperature, over some period of time. While the first technique could be implemented using a number of different simulation approach, the latter is particularly well-suited to the strengths of the MD approach.


Page 19 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile Given the overall capabilities of reactive MD demonstrated in this work, it is now feasible to investigate these additional approaches in more detail.

5 Conclusions Here, we have demonstrated the use of reactive empirical potentials (ReaxFF, and COMB3) in MD simulations to investigate the structure and energy of surfaces and GBs in oxide systems. With few exceptions, reactive MD potentials successfully recover the relative energy of γsurf and γGB predicted by ab initio techniques. Furthermore, γsurf for {h0l} and {hk0} surfaces are consistent with surface energies calculated for terraced surfaces composed of low-energy terraces; these are generally the lowest-energy terminations observed. γGB is observed to be highly dependent on rigid-body translations between neighboring grains in the plane of the GB. This energy topography recovers expected symmetry of the GB plane, and is similar, but not identical between different MD potentials. Finally, reactive MD successfully identifies energy cusps associated with low-energy interfaces (e.g., those observed at coherent twin boundaries), allowing for higher-fidelity calculation of grain γGB as a function of microscopic parameters. Thus, reactive MD potentials appear capable of successfully reproducing the key features anticipated for GB structures and energetics. This study serves as a foundation to further utilization of reactive MD techniques to investigate the structure, properties, and dynamics of GBs in oxide systems. While such techniques have been demonstrated previously for cubic metallic systems, they have not yet been applied oxides, which generally have lower symmetry unit cells (tetragonal, orthorhombic, monoclinic, etc.) and the added complication of charged species which may exhibit partial charge transfer to reduce the energy of the system. This study has demonstrated the successful use of reactive MD potentials to describe such systems. We anticipate future studies employing these potentials 1) to screen GBs as a function of all eight orientation parameters coupled with secondary quantitative ab initio approaches, 2) to identify unique low-energy GBs which have special significance in polycrystalline systems, 3) to investigate dynamic transient processes which occur over timescales that are too long to simulate using higher fidelity approaches, specifically including 4) ionic and molecular transport along GBs.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile Supporting Information The Supporting Information is available free of charge on the ACS Publications website at DOI: . Figures illustrating computational details including computational cell, validation of cell size and GB spacing, and example relaxation of GB (PDF).

Acknowledgements We would like to thank Prof. S.B. Sinnott, Prof. A. van Duin, and Dr. T. Liang for informative discussions during the preparation of this manuscript. We would like to thank AFOSR for support of this work (14RX13COR), and the Texas A&M Supercomputing Facility (http://sc.tamu.edu) for computational facilities used in this study.


Page 21 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile References 1. Barrales-Mora, L. A., 2D and 3D Grain Growth Modeling and Simulation; Cuvillier Verlag, 2008. 2. Lejček, P., Grain Boundaries: Description, Structure and Thermodynamics. In Grain Boundary Segregation in Metals, Springer: 2010; pp 5-24. 3. Ohtomo, A.; Hwang, H., A High-Mobility Electron Gas at the LaAlO3/SrTiO3 Heterointerface. Nature 2004, 427, 423-426. 4. Hwang, H.; Iwasa, Y.; Kawasaki, M.; Keimer, B.; Nagaosa, N.; Tokura, Y., Emergent Phenomena at Oxide Interfaces. Nat. Mater. 2012, 11, 103-113. 5. Dikin, D.; Mehta, M.; Bark, C.; Folkman, C.; Eom, C.; Chandrasekhar, V., Coexistence of Superconductivity and Ferromagnetism in Two Dimensions. Phys. Rev. Lett. 2011, 107, 056802. 6. Saylor, D. M.; Dasher, B.; Pang, Y.; Miller, H. M.; Wynblatt, P.; Rollett, A. D.; Rohrer, G. S., Habits of Grains in Dense Polycrystalline Solids. J. Am. Ceram. Soc. 2004, 87, 724-726. 7. Saylor, D. M.; Morawiec, A.; Rohrer, G. S., Distribution of Grain Boundaries in Magnesia as a Function of Five Macroscopic Parameters. Acta Mater. 2003, 51, 3663-3674. 8. Pang, Y.; Wynblatt, P., Correlation between Grain‐Boundary Segregation and Grain‐ Boundary Plane Orientation in Nb‐Doped TiO2. J. Am. Ceram. Soc. 2005, 88, 2286-2291. 9. Pang, Y.; Wynblatt, P., Effects of Nb Doping and Segregation on the Grain Boundary Plane Distribution in TiO2. J. Am. Ceram. Soc. 2006, 89, 666-671. 10. Dawson, I.; Bristowe, P.; Lee, M.-H.; Payne, M.; Segall, M.; White, J., First-Principles Study of a Tilt Grain Boundary in Rutile. Phys. Rev. B: Condens. Matter 1996, 54, 13727. 11. Körner, W.; Elsässer, C., Density Functional Theory Study of Dopants in Polycrystalline TiO2. Phys. Rev. B: Condens. Matter 2011, 83, 205315. 12. Hallil, A.; Amzallag, E.; Landron, S.; Tétot, R., Properties of Rutile TiO2 Surfaces from a Tight-Binding Variable-Charge Model. Comparison with Ab Initio Calculations. Surf. Sci. 2011, 605, 738-745. 13. Rittner, J.; Seidman, D., Symmetric Tilt Grain-Boundary Structures in FCC Metals with Low Stacking-Fault Energies. Phys. Rev. B: Condens. Matter 1996, 54, 6999. 14. Tschopp, M.; McDowell, D., Structures and Energies of Σ3 Asymmetric Tilt Grain Boundaries in Copper and Aluminium. Philos. Mag. 2007, 87, 3147-3173. 15. Tschopp, M.; McDowell, D., Asymmetric Tilt Grain Boundary Structure and Energy in Copper and Aluminium. Philos. Mag. 2007, 87, 3871-3892. 16. Gemming, S.; Janisch, R.; Schreiber, M.; Spaldin, N., Density-Functional Investigation of the (113)[− 110] Twin Grain Boundary in Co-Doped Anatase TiO2 and its Influence on Magnetism in Dilute Magnetic Semiconductors. Phys. Rev. B: Condens. Matter 2007, 76, 045204. 17. Glassford, K. M.; Chelikowsky, J. R., Structural and Electronic Properties of Titanium Dioxide. Phys. Rev. B: Condens. Matter 1992, 46, 1284. 18. Liang, T.; Shin, Y. K.; Cheng, Y.-T.; Yilmaz, D. E.; Vishnu, K. G.; Verners, O.; Zou, C.; Phillpot, S. R.; Sinnott, S. B.; van Duin, A. C., Reactive Potentials for Advanced Atomistic Simulations. Annu. Rev. Mater. Res. 2013, 43, 109-129. 19. Shin, Y. K.; Shan, T.-R.; Liang, T.; Noordhoek, M. J.; Sinnott, S. B.; van Duin, A. C.; Phillpot, S. R., Variable Charge Many-Body Interatomic Potentials. MRS Bull. 2012, 37, 504512. 21 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile 20. Rappe, A. K.; Goddard III, W. A., Charge Equilibration for Molecular Dynamics Simulations. J. Phys. Chem. 1991, 95, 3358-3363. 21. Janssens, G. O.; Baekelandt, B. G.; Toufar, H.; Mortier, W. J.; Schoonheydt, R. A., Comparison of Cluster and Infinite Crystal Calculations on Zeolites with the Electronegativity Equalization Method (EEM). J. Phys. Chem. 1995, 99, 3251-3258. 22. Mortier, W. J.; Ghosh, S. K.; Shankar, S., Electronegativity-Equalization Method for the Calculation of Atomic Charges in Molecules. J. Am. Chem. Soc. 1986, 108, 4315-4320. 23. Liang, T.; Devine, B.; Phillpot, S. R.; Sinnott, S. B., Variable Charge Reactive Potential for Hydrocarbons to Simulate Organic-Copper Interactions. J. Phys. Chem. A 2012, 116, 79767991. 24. van Duin, A. C.; Dasgupta, S.; Lorant, F.; Goddard, W. A., Reaxff: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396-9409. 25. van Duin, A. C.; Strachan, A.; Stewman, S.; Zhang, Q.; Xu, X.; Goddard, W. A., Reaxffsio Reactive Force Field for Silicon and Silicon Oxide Systems. J. Phys. Chem. A 2003, 107, 3803-3811. 26. Plimpton, S., Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1-19. 27. Aktulga, H. M.; Fogarty, J. C.; Pandit, S. A.; Grama, A. Y., Parallel Reactive Molecular Dynamics: Numerical Methods and Algorithmic Techniques. Parallel Comput. 2012, 38, 245259. 28. Liang, T.; Shan, T.-R.; Cheng, Y.-T.; Devine, B. D.; Noordhoek, M.; Li, Y.; Lu, Z.; Phillpot, S. R.; Sinnott, S. B., Classical Atomistic Simulations of Surfaces and Heterogeneous Interfaces with the Charge-Optimized Many Body (COMB) Potentials. Mater. Sci. Eng., R 2013, 74, 255-279. 29. Lewis, G.; Catlow, C., Potential Models for Ionic Oxides. J. Phys. C: Solid State Phys. 1985, 18, 1149. 30. Matsui, M.; Akaogi, M., Molecular Dynamics Simulation of the Structural and Physical Properties of the Four Polymorphs of TiO2. Mol. Simul. 1991, 6, 239-244. 31. Collins, D. R.; Smith, W.; Harrison, N. M.; Forester, T. R., Molecular Dynamics Study of TiO2 Microclusters. J. Mater. Chem. 1996, 6, 1385-1390. 32. Swamy, V.; Gale, J. D.; Dubrovinsky, L., Atomistic Simulation of the Crystal Structures and Bulk Moduli of TiO2 Polymorphs. J. Phys. Chem. Solids 2001, 62, 887-895. 33. Hallil, A.; Tétot, R.; Berthier, F.; Braems, I.; Creuze, J., Use of a Variable-Charge Interatomic Potential for Atomistic Simulations of Bulk, Oxygen Vacancies, and Surfaces of Rutile TiO2. Phys. Rev. B: Condens. Matter 2006, 73, 165406. 34. Tétot, R.; Hallil, A.; Creuze, J.; Braems, I., Tight-Binding Variable-Charge Model for Insulating Oxides: Application to TiO2 and ZrO2 Polymorphs. EPL 2008, 83, 40001. 35. Cheng, Y.-T.; Shan, T.-R.; Liang, T.; Behera, R. K.; Phillpot, S. R.; Sinnott, S. B., A Charge Optimized Many-Body (COMB) Potential for Titanium and Titania. J. Phys.: Condens. Matter 2014, 26, 315007. 36. Kim, S.-Y.; Kumar, N.; Persson, P.; Sofo, J.; van Duin, A. C.; Kubicki, J. D., Development of a Reaxff Reactive Force Field for Titanium Dioxide/Water Systems. Langmuir 2013, 29, 7838-7846. 37. Raju, M.; Kim, S.-Y.; van Duin, A. C.; Fichthorn, K. A., Reaxff Reactive Force Field Study of the Dissociation of Water on Titania Surfaces. J. Phys. Chem. C 2013, 117, 1055810572. 22 ACS Paragon Plus Environment

Page 23 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile 38. Yu, J.; Sinnott, S.; Phillpot, S., Optimized Many Body Potential for FCC Metals. Philos. Mag. Lett. 2009, 89, 136-144. 39. Shan, T.-R.; Devine, B. D.; Hawkins, J. M.; Asthagiri, A.; Phillpot, S. R.; Sinnott, S. B., Second-Generation Charge-Optimized Many-Body Potential for Si/SiO2 and Amorphous Silica. Phys. Rev. B: Condens. Matter 2010, 82, 235302. 40. Cheng, Y.-T.; Shan, T.-R.; Devine, B.; Lee, D.; Liang, T.; Hinojosa, B. B.; Phillpot, S. R.; Asthagiri, A.; Sinnott, S. B., Atomistic Simulations of the Adsorption and Migration Barriers of Cu Adatoms on ZnO Surfaces Using Comb Potentials. Surf. Sci. 2012, 606, 1280-1288. 41. Devine, B.; Shan, T.-R.; Cheng, Y.-T.; McGaughey, A. J.; Lee, M.; Phillpot, S. R.; Sinnott, S. B., Atomistic Simulations of Copper Oxidation and Cu/Cu2O Interfaces Using Charge-Optimized Many-Body Potentials. Phys. Rev. B: Condens. Matter 2011, 84, 125308. 42. Shan, T.-R.; Devine, B. D.; Kemper, T. W.; Sinnott, S. B.; Phillpot, S. R., ChargeOptimized Many-Body Potential for the Hafnium/Hafnium Oxide System. Phys. Rev. B: Condens. Matter 2010, 81, 125328. 43. Chenoweth, K.; van Duin, A. C.; Goddard, W. A., Reaxff Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 10401053. 44. Valentini, P.; Schwartzentruber, T. E.; Cozmuta, I., Molecular Dynamics Simulation of O2 Sticking on Pt (111) Using the Ab Initio Based Reaxff Reactive Force Field. J. Chem. Phys. 2010, 133, 084703. 45. Chenoweth, K.; van Duin, A. C.; Persson, P.; Cheng, M.-J.; Oxgaard, J.; Goddard, I., William A, Development and Application of a Reaxff Reactive Force Field for Oxidative Dehydrogenation on Vanadium Oxide Catalysts. J. Phys. Chem. C 2008, 112, 14645-14654. 46. Zhang, Q.; Çaǧın, T.; van Duin, A.; Goddard III, W. A.; Qi, Y.; Hector Jr, L. G., Adhesion and Nonwetting-Wetting Transition in the Al/α-Al2O3 Interface. Phys. Rev. B: Condens. Matter 2004, 69, 045423. 47. Russo, M. F.; Li, R.; Mench, M.; van Duin, A. C., Molecular Dynamic Simulation of Aluminum–Water Reactions Using the Reaxff Reactive Force Field. Int. J. Hydrogen Energy 2011, 36, 5828-5835. 48. Labat, F.; Baranek, P.; Adamo, C., Structural and Electronic Properties of Selected Rutile and Anatase TiO2 Surfaces: An Ab Initio Investigation. J. Chem. Theory Comput. 2008, 4, 341352. 49. Swamy, V.; Muscat, J.; Gale, J. D.; Harrison, N. M., Simulation of Low Index Rutile Surfaces with a Transferable Variable-Charge Ti–O Interatomic Potential and Comparison with Ab Initio Results. Surf. Sci. 2002, 504, 115-124. 50. Scaranto, J.; Mallia, G.; Giorgianni, S.; Zicovich-Wilson, C.; Civalleri, B.; Harrison, N., A Quantum-Mechanical Study of the Vinyl Fluoride Adsorbed on the Rutile TiO2 (110) Surface. Surf. Sci. 2006, 600, 305-317. 51. Bredow, T.; Giordano, L.; Cinquini, F.; Pacchioni, G., Electronic Properties of Rutile TiO2 Ultrathin Films: Odd-Even Oscillations with the Number of Layers. Phys. Rev. B: Condens. Matter 2004, 70, 035419. 52. Muscat, J.; Harrison, N., The Physical and Electronic Structure of the Rutile (001) Surface. Surf. Sci. 2000, 446, 119-127. 53. Muscat, J.; Harrison, N.; Thornton, G., Effects of Exchange, Correlation, and Numerical Approximations on the Computed Properties of the Rutile TiO2 (100) Surface. Phys. Rev. B: Condens. Matter 1999, 59, 2320. 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile 54. Hameeuw, K.; Cantele, G.; Ninno, D.; Trani, F.; Iadonisi, G., The Rutile TiO2 (110) Surface: Obtaining Converged Structural Properties from First-Principles Calculations. J. Chem. Phys. 2006, 124, 024708. 55. Goniakowski, J.; Holender, J.; Kantorovich, L.; Gillan, M.; White, J., Influence of Gradient Corrections on the Bulk and Surface Properties of TiO2 and SnO2. Phys. Rev. B: Condens. Matter 1996, 53, 957. 56. Ramamoorthy, M.; Vanderbilt, D.; King-Smith, R., First-Principles Calculations of the Energetics of Stoichiometric TiO2 Surfaces. Phys. Rev. B: Condens. Matter 1994, 49, 1672116727. 57. Morgan, B. J.; Watson, G. W., A DFT+U Description of Oxygen Vacancies at the TiO2 Rutile (110) Surface. Surf. Sci. 2007, 601, 5034-5041. 58. Perron, H.; Domain, C.; Roques, J.; Drot, R.; Simoni, E.; Catalette, H., Optimisation of Accurate Rutile TiO2 (110),(100),(101) and (001) Surface Models from Periodic DFT Calculations. Theor. Chem. Acc. 2007, 117, 565-574. 59. Thompson, S. J.; Lewis, S. P., Revisiting the (110) Surface Structure of TiO2: A Theoretical Analysis. Phys. Rev. B: Condens. Matter 2006, 73, 073403. 60. Vijay, A.; Mills, G.; Metiu, H., Adsorption of Gold on Stoichiometric and Reduced Rutile TiO2 (110) Surfaces. J. Chem. Phys. 2003, 118, 6536-6551. 61. Lindan, P.; Harrison, N.; Gillan, M.; White, J., First-Principles Spin-Polarized Calculations on the Reduced and Reconstructed TiO2 (110) Surface. Phys. Rev. B: Condens. Matter 1997, 55, 15919. 62. Bates, S.; Kresse, G.; Gillan, M., A Systematic Study of the Surface Energetics and Structure of TiO2 (110) by First-Principles Calculations. Surf. Sci. 1997, 385, 386-394. 63. Lee, W.-Y.; Bristowe, P.; Gao, Y.; Merkle, K., The Atomic Structure of Twin Boundaries in Rutile. Philos. Mag. Lett. 1993, 68, 309-314. 64. Noguera, C., Polar Oxide Surfaces. J. Phys.: Condens. Matter 2000, 12, R367. 65. Tasker, P., The Stability of Ionic Crystal Surfaces. J. Phys. C: Solid State Phys. 1979, 12, 4977. 66. Van Hove, M.; Somorjai, G., A New Microfacet Notation for High-Miller-Index Surfaces of Cubic Materials with Terrace, Step and Kink Structures. Surf. Sci. 1980, 92, 489-518. 67. Cahn, J. W.; Handwerker, C. A., Equilibrium Geometries of Anisotropic Surfaces and Interfaces. Mater. Sci. Eng., A 1993, 162, 83-95. 68. Daneu, N.; Schmid, H.; Rečnik, A.; Mader, W., Atomic Structure and Formation Mechanism of (301) Rutile Twins from Diamantina (Brazil). Am. Mineral. 2007, 92, 1789-1799. 69. Lu, W.; Bruner, B.; Casillas, G.; Mejía-Rosales, S.; Farmer, P. J.; José-Yacamán, M., Direct Oxygen Imaging in Titania Nanocrystals. Nanotechnol. 2012, 23, 335706. 70. Schelling, P.; Yu, N.; Halley, J., Self-Consistent Tight-Binding Atomic-Relaxation Model of Titanium Dioxide. Phys. Rev. B: Condens. Matter 1998, 58, 1279.

Tables Table 1. Summary of reported surface formation energies for TiO2 (Rutile), in J·m-2, after ref. [12]. Method (110) (101) (100) (001) Ref. ab initio


Page 25 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile LCAO-based HF (Crystal03) HF (Crystal98) HF (Crystal95) Hybrid-B3LYP (Crystal06) Hybrid-B3LYP (Crystal03) Hybrid-B3LYP (Crystal03) Hybrid-PBEO (Crystal03) LDA (Crystal03) LDA (Crystal98) LDA (Crystal95) GGA-PBE (Crystal03) GGA (Crystal03) GGA (Crystal95)

0.92 1.00

1.13

0.70 0.83 1.20

0.69

1.45 1.59 1.88 1.87 1.30 1.39 1.48 0.83

48

48 48 49 53 48 51-52 53

55

1.12

1.65

56 57 57

0.70

1.25

58 59 60 61 62 55

0.68

1.36

12

0.78

1.02

1.46

70

0.95

1.24

2.20

29, 32 63

1.33

Qeq-based (charge equilibration) MS-Q SMB-Q SMTB-Q COMB3 ReaxFF

53

54

0.91 1.14 0.89 0.58 0.83 0.50 0.54 0.40-0.57 0.81 0.73 0.82 0.48

Empirical Potentials Ionic (fixed-charge) CFR CFR MA

49, 53

50

0.42 0.64

Semi-empirical, DFTB Self-consistent Tight Binding

48

12

0.49 0.40 1.10 0.55 0.90 0.90

PW-based LDA (QU-EXPR) LDA (CASTEP) LDA (?) GGA (VASP) GGA+U (VASP) GGA (VASP) GGA (VASP) GGA (VASP) GGA (CETEP) GGA (VASP) GGA (CASTEP) GGA (Crystal06)

2.08 2.20 1.21

1.77

2.07

2.40

30, 32

1.40 0.60 0.42 0.99 0.29

1.32 0.85 0.49 1.25

1.85 1.74 1.26 1.49

48

0.71

33 12 35 36



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile COMB3

1.01

1.07

1.37

1.26 this study

ReaxFF

0.55

1.52

0.73

2.41

a

2.05

b

0.69 a

Ti/O/H parameters from Raju et al. (2013)

b

Ti/O/H parameters from Kim et al. (2013) 36

1.45

0.80

, this study , this study

37


Page 27 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48


Shamberger et al. Investigating GB Structures and Energetics of Rutile Table 2. Surface formation energies for low-index TiO2 (Rutile), in J·m-2. COMB3 results

ReaxFF resultsa

ReaxFF resultsb

Termination /

surf

surf

surf

Tasker Class.

Surface

γ

-2

(100) (001) (101) (110)

(unc.) -2

γ

-2

(unc.) -2

γ

-2

(unc.) -2

J·m

J·m

J·m

J·m

J·m

J·m

1.37 2.75 1.26 1.06 2.98 1.01 1.95

0.04 0.10 0.03 0.05 0.08 0.05 0.11

0.73 2.86 2.41 1.66 3.44 0.55 2.09

0.05 0.05 0.03 0.04 0.04 0.05 0.05

0.80 2.11 2.05 1.45 3.03 0.69 1.75

0.05 0.05 0.03 0.04 0.04 0.05 0.05

a

Ti/O/H parameters from Raju et al. (2013) 37

b


c

initial, unrelaxed dimensions

II III I II III II III

non-polar polar non-polar non-polar polar non-polar polar

Asurfc

dslabc

Å2

Å

244.7

41.3

486

379.9 451.9

26.6 34.8

486 756

346.0

45.5

756

Natoms



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

Page 28 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile Table 3. Family of {h0l} free surfaces. Surface

surf

γ

-2

J·m (100) (501) (401)

d

(301) (201) (101) (102) d (103) (104) (105) (001)

d

ReaxFF resultsa

ReaxFF resultsb

surf

surf

COMB3 results (unc.) -2

J·m

γ

-2

J·m

(unc.) -2

J·m

γ

-2

J·m

J·m-2

Å2

Å

Natoms

244.7 854.0

41.3 35.1

972 1440

* 1 * *

II III II II II

non-polar polar non-polar non-polar non-polar

699.9

34.3

1152

551.0 413.0 451.9

32.6 29.0 34.8

864 576 1512

0.04 0.04 0.04

* 1 *

II III II

non-polar polar non-polar

798.2

33.8

1296

388.6

34.7

648

0.04 0.04 0.03 0.03

* 1 *

II III II I

non-polar polar non-polar non-polar

513.1

35.1

864

638.4 379.9

29.3 26.6

900 972

0.05 0.04

0.80 1.21

0.05 0.04

1.36 1.97 1.33 1.26 1.07

0.04 0.04 0.04 0.04 0.04

1.52 1.63 1.69 1.49 1.66

0.04 0.04 0.04 0.03 0.04

1.27 1.27 1.35 1.43 1.45

0.04 0.04 0.04 0.03 0.04

1.20 2.16 1.23

0.04 0.04 0.04

2.08 2.80 2.04

0.04 0.04 0.04

1.75 2.32 1.81

1.24 1.73 1.25 1.27

0.04 0.04 0.04 0.02

2.01 2.46 2.04 2.41

0.04 0.04 0.03 0.03

1.86 2.17 1.90 2.05


dslabc

non-polar non-polar

0.73 1.24

b

Asurfc

II II

0.03 0.04


Tasker Class.

(unc.)

1.38 1.37

a

termin.

*

c

initial, unrelaxed dimensions require minor restructuring to obtain the lowest energy surface. * terraced surface termination

d


Page 29 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48


Shamberger et al. Investigating GB Structures and Energetics of Rutile Table 4. Family of {hk0} free surfaces. Surface

COMB3 results

ReaxFF resultsa

ReaxFF resultsb

surf

surf

surf

γ

-2

J·m (100) (510) d

(410) d

(310) d

(210) d

(110)

1.37 1.35 2.12 2.36 2.38 1.34 2.27 2.29 2.53 1.32 2.11 2.11 2.54 1.25 1.81 2.24 2.28 2.32 1.01

(unc.) -2

J·m

0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.05 0.05 0.05 0.05 0.05 0.05 0.04 0.05 0.04 0.05

γ

-2

J·m

(unc.) -2

J·m

γ

-2

J·m

J·m-2

0.05 0.04

0.80 0.83

0.05 0.04

2.64

0.04

1.75

0.04

0.72

0.04

0.82

0.04

2.49

0.04

1.82

0.04

0.71 2.25

0.05 0.05

0.81 1.50

0.05 0.05

0.68

0.05

0.79

0.05

2.23

0.05

1.30

0.05

0.56

0.04

0.69

0.05


b


c

initial, unrelaxed dimensions.

Tasker Class.

Asurfc

dslabc

Å2

Å

Natoms

(unc.)

0.73 0.73

a

termin.

* 2 1 3 * 3 1 2 * 1 3 2 * 3 1 2 4

II II III II III

non-polar non-polar polar non-polar polar

III III III II II III III II III II III II II

polar polar polar non-polar non-polar polar polar non-polar polar non-polar polar non-polar non-polar

244.7 831.8

41.3 36.0

972 1440

672.6

35.7

1152

515.8

43.6

1080

547.1

41.1

1080

346.0

45.5

1512

d

require minor restructuring to obtain the lowest energy surface. * terraced surface termination



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile Table 5. Grain Boundary Formation Energies for TiO2 (Rutile), in J·m-2. GB GB Σ Axis Plane GB Type desig.a Method b

Empirical (CFR - fixed charge) Empirical (COMB3 - QEq) Empirical (ReaxFF - QEq) Empirical (CFR - fixed charge) Empirical (COMB3 - QEq) Empirical (ReaxFF - QEq) ab initio (PW, LDA) Empirical (COMB3 - QEq) Empirical (ReaxFF - QEq)

0.42 0.32 0.18 0.47 0.51 0.12

[010]

(101)

rotation twin / pseudo-reflection twin

Σ2

[010]

(301)

reflection twin / pseudo-rotation twin

Σ2 b

reflection twin

Σ1

[010]

(100)

γGB

Ref.

0.124 63 0.074 this study 0.294 this study 63

this study this study 11

this study this study

[001]

(210)

tilt boundary

Σ5 c

ab initio (PW, LDA) Empirical (fixed charge) ab initio (PW, SIC-LDA) Empirical (COMB3 - QEq) Empirical (ReaxFF - QEq)

1.72 10 1.70 63 1.92 11 1.16 this study 0.40 this study

[001]

(310)

tilt boundary

Σ5 c

ab initio (PW, LDA) ab initio (PW, SIC-LDA) Empirical (COMB3 - QEq) Empirical (ReaxFF - QEq)

2.37 11 2.84 11 1.34 this study 0.70 this study

a

index on Ti-basis.

b

near-CSL designation.

c

CSL index prior to in-plane translation.


Page 31 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Shamberger et al. Investigating GB Structures and Energetics of Rutile Table 6. Family of {h0l} refl+trans and rot+trans GBs. GB axis GB plane COMB3 ReaxFF a refl+trans rot+trans refl+trans

γGB

γGB

γGB

J·m-2

J·m-2

J·m-2

[010] [010]

(100) (501)

0.51 1.82

0.12 0.59

[010]

(401) b

0.53

[010] [010] [010]

(301) (201) (101)

1.44 1.83 0.32 2.03 0.08

0.32

termin.

* 1

0.18 0.75 0.29

0.08

b

1.50 0.99 1.35

2.29

1.13 0.89 0.83

* 1

1.97 1.67 1.43 0.00

1.84

1.22

* 1

[010]

(102)

[010]

(103)

[010]

(104) b

[010] [010]

(105) (001)

0.87 0.00

a


b

requires minor restructuring to obtain the lowest energy surface.

Table 7. Family of {hk0} refl+trans GBs. GB axis

GB plane

COMB3

γ [001] [001]

(100) (310)

[001]

(210)

a

ReaxFFc

GB

γ

termin.

GB

-2

J·m

J·m-2

0.51 1.34 1.44 1.80 1.16 1.28 1.28 1.28

0.12 0.79 0.70 0.83 0.40 0.62 0.62 0.62

1 2 3 2 1 3 4




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 53

Shamberger et al. Investigating GB Structures and Energetics of Rutile Figures

Figure 1. Computational cell used in this study. Cell illustrates refl+trans across a plane {hk0}, where the interfacial plane is determined by some rotation around [001], accompanied by rigidbody displacement (& , & ).

Figure 2. Surface termination in low Miller-index planes. Illustrated terminations are non-polar terminations. Polar terminations are illustrated with solid blue lines. 32 ACS Paragon Plus Environment

Page 33 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Figure 3. Average charges on oxygen atoms (-qO), and on titanium atoms (+qTi) for both nonpolar terminations (blue solid lines) and polar terminations (red dashed lines), as a function of slab depth, d, using a) COMB3 potentials, and b) ReaxFF potentials.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 53


Figure 4. Electric fields (solid black lines) and potential energy distributions (dashed red line) within non-polar and polar terminations, as a function of slab depth, d, using a) COMB3 potentials, and b) ReaxFF potentials.


Page 35 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Figure 5. Terminations investigated for {h0l} surfaces. Periodic in-plane repeat distances are illustrated with horizontal brackets. Termination with (101) and (001) or (100)-type surfaces are illustrated with thin grey lines. In those planes marked with an asterix, geometry necessitates some surface restructuring in order to have a terraced surface. The original position of these atoms (indicated by shaded circle) and the final position of these atoms (illustrated by arrows) are both illustrated.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 53


Figure 6. Terminations investigated for {hk0} planes. Periodic in-plane repeat distances are illustrated. Termination with (110) and (100)-type surfaces are illustrated with thin grey lines. In those planes marked with an asterix, geometry necessitates some surface restructuring in order to have a terraced surface; this is accomplished by removing those atoms which fall above the grey line. In contrast, different planar terminations considered are showed with blue lines; these are numbered in descending order from outermost termination (1) to innermost (either 3 or 4, based on the number of unique terminations).


Page 37 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Figure 7. γsurf of a) {h0l} and b) {hk0} surfaces. x’s indicate calculated surfaces using COMB3 MD potentials, o’s indicate calculated surfaces using ReaxFF MD potentials. Multiple terminations are considered for most surface orientations (Table 3, 4); those terminations composed of terraced stable low Miller-index surfaces are connected by a solid line. Dashed (red) line illustrates those surface energies calculated from weighted averages of low Millerindex surfaces as described in the text (eq. 3).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 53


Figure 8. Histogram of γGB (a, c) and distributions of γGB with initial displacement in the plane of the grain boundary (b, d) along the (301) plane using COMB3 MD potentials; (a, b) refl+trans, (c, d) rot+trans.


Page 39 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Figure 9. a) Histogram of γGB and b) distributions of γGB with initial displacement in the plane of the grain boundary along the (401) plane for refl+trans using COMB3 MD potentials.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 53


Figure 10. Histogram of γGB (a, c) and distributions of γGB with initial displacement in the plane of the grain boundary (b, d) along the (105) plane using the COMB3 (a, b), and ReaxFF (c, d) MD potentials. Relaxation pathways of those initial displacements in the dashed box in the lower left corner are illustrated in the supplementary information section.


Page 41 of 53

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Figure 11. Grain boundary energies of {h0l} pseudo-reflection twin grain boundaries as a function of the angle between [001] in each lattice. Blue x’s are calculated using COMB3 potentials, red circles are calculated using ReaxFF potentials. Lines connect those GB terminations with the lowest γGB. Low-energy boundaries (energy cusps) are identified. Solid and dashed lines are calculated using Eq. 4 in text.



(301)

(001)

(101)

(301)

1.5

(100)

2.0

γGB (J/m2)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1.0 0.5 0.0 0

(105) 20 40 60 80 100 120 140 160 180

θ (deg)

ACS Paragon Plus Environment

Page 42 of 53

Page 43 of 53


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

GB1 y

z

a2′

a2

x

Cell Reference Frame

a1′ a1

a3′

GB2 δz

Rotation Axis

δx a1 a1′′

a3′′ a2′′

a2

GB1



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19[100] 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(100)

(101)

(001)

[001] [001] [010]

Page 44 of 53

[100]


(110)

[100]

[100]

[010]

Page 45 of 53

-qO, qTi

2.0

-qO, qTi

b)

1.0

2.0

(100)

(100) Nonpolar Termination -q O -q Ti Polar Termination -q O -q Ti

1.5 1.0 0.5 2.0

-qO, qTi

a)

1.5

0.5

(001)

(001)

1.5 1.0 0.5 2.0

-qO, qTi

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(101)

(101)

(110)

(110)

1.5 1.0 0.5 0.0

0

1

2

d / nm

3

4

5 0


1

2

d / nm

3

4

5


a) non-polar

polar

E, V

E, V

polar (100)

0 0

0

E, V

E, V

0

(001)

E V

0

E, V

E, V

0

(101)

0

0

0

E, V

0 E, V

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 53

b) non-polar

(110)

0

0

0

0 0

1

2

3

d / nm

4

5 0

1

2

3

d / nm

4

5

0

1

2

3

d / nm


4

5 0

1

2

3

d / nm

4

5

Page 47 of 53


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(105) (104) * (103) (102) * (101) (201) (301) (401) * (501)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1 3

(510) * 1 3

(410) *

1 3

(310) *

1 4

(210) *

(110)


Page 48 of 53

Page 49 of 53

3.0

a)

(001)

(101)

2.0

(100)

2

γsurf (J/m )

2.5

1.5 1.0 0.5

0

3.0

10

20

30

10

15

40

50

60

70

80

20

25

30

35

40

θ (deg)

90

b)

2.5 2.0

(100)

1.5

(010)

2

γsurf (J/m )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.0 0.5

0

5

θ (deg)

45



100

[010] displacement (nm)

frequency

γ (J/m2)

a)

γmin = 0.319 J/m2

80 60 40 20 0 100

γmin = 0.320 J/m


80 60 40 20 0

0.5

1

1.5

2

γGB (J/m2)

2

0.2

1.5

0.1

1 0.5 0

2.5

3

3.5

0.2

0.4

0.6

0.8


1.0

γ (J/m2) 3

d)

0.4

2.5 0.3

2

0.2

1.5

0.1

1

0 0

2.5

0.3

c)

2

3

b)

0.4

0

frequency

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 50 of 53

0.5 0

0.2

0.4

0.6

0.8



1.0

Page 51 of 53

12

γmin = 1.443 J/m2

0.5


10 8

frequency

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6 4 2

0.4

0

0.5

1

1.5

γGB

2

(J/m2)

2.5

3

2.5

0.3 0.2

2

0.1 0

0

3

1.5 0

3.5


0.5

1.0


1.5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

a)

Page 52 of 53

γmin = 1.43 J/m

2


b) 0.4 0.3 0.2 0.1 0

c)

0

0.5

1.0

1.5

2

0

0.5

1.0

1.5

2

γmin = 0.87 J/m2


d)

0.5

0.4 0.3 0.2 0.1 0

1

1.5

γGB

2

(J/m2)

[501] displacement (

2.5ACS Paragon 3 Plus 3.5Environment

[501] displacement (

Page 53 of 53

3.0 2.5 2.0

(001)

1.5

(101)

(301)

(100)

2

γGB (J/m )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.0 0.5 0.0

0

20

40

60

80 100 120 140 160 180

θ (deg)


Investigating Grain Boundary Structures and Energetics of Rutile with

Recommend Documents