Comparative first-principles structural and vibrational properties of rutile and anatase TiO2

The structural and vibrational properties of two polymorphs of TiO2, rutile and anatase, have been investigated by first-principles methods at different levels of exchange-correlational (XC) energy functionals in density functional theory (DFT). Reports in the literature to date are contradictory regarding the stability of the rutile phase using DFT XC-functionals more sophisticated than simple local-density approximation. Here the PBEsol generalized gradient approximation (GGA), TPSS meta-GGA, and HSE06 hybrid functionals have been employed to demonstrate the XC-functional effects on the calculated structural, phonon and thermodynamic properties of rutile and anatase TiO2. Lattice and elastic parameters correctly calculated with these XC-functionals show good agreement with the experimental values. Calculated phonon frequencies generated stable phonon dispersion relations for both rutile and anatase TiO2 when correctly converged, in agreement with the experimental observations. The phonon frequencies along high symmetry Brillouin zone paths and their corresponding phonon density of states showed sensitivity to different levels of XC-functional employed in phonon dispersion prediction. Nevertheless, the thermodynamic properties of rutile and anatase TiO2 estimated by harmonic approximations are in excellent experimental agreement and are effectively invariant to the level of theory employed in the DFT XC-functional.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Titanium dioxide (TiO 2 ) is one of the most investigated materials in materials science and surface science research communities due to its novel structural, electronic, optical and thermal properties, and their potential applications.It has been widely used as a heterogeneous catalyst in a variety of energy and environmental applications including photocatalysis [1][2][3] and photoelectrochemical water splitting [4,5], solar cells for hydrogen and electric energy production [6,7], gas sensors [8], air [9] and water purification [5,10], high-performance dielectric materials [11], optical coating for self-cleaning glasses [12], white pigment [13] just to name a few.The enhanced performance of TiO 2 in such various areas derives from its strong oxidizing abilities, chemical stability, nontoxicity, transparency to visible light and low cost [14,15].TiO 2 is mainly found in nature in three major crystal forms: rutile, anatase and brookite, with varying degree of distortion in their TiO 6 octahedral arrangements [16].Among these three phases, rutile is thermodynamically the most stable to 1800 K [17], while the other two phases, anatase and brookite are meta stable and both covert to rutile upon heating [18,19].Furthermore, Zhang and Banfield [18] reported size-dependent phase transformation of these three phases in equally sized particles, where rutile is found to be the most stable phase at above 35 nm and anatase being thermodynamically stable below 11 nm.Brookite forms the most stable crystal phase between 11 and 35 nm.However, the preparation of brookite form is challenging as it is always accompanied by secondary phases such as anatase and/or rutile with a low structural symmetry [20][21][22].
The phonon modes of a solid material determine its thermodynamic properties including temperature-dependent volumetric and constant-pressure heat capacities, thermal conductivity, entropy, relative molar enthalpy, Helmholtz and Gibbs free energies [54].They also provide important optical properties such as dielectric function, infrared and Raman spectra [55].Additionally, phonon modes are often used to describe the physical properties of a material, for instance structural stability, elastic constants, superconductivity, ferroelectricity, and phase transitions.Therefore, precise measurement of phonon modes is the key in order to understand any of these phonon-intimated properties of rutile and anatase TiO 2 .
Traylor et al [56] experimentally measured stable phonon dispersion relations of rutile TiO 2 along the principal directions in the Brillouin zone (BZ) using coherent inelastic neutron scattering method and further compared those relations by force-field approaches such as rigid-ion model and shell model.They found qualitative agreement for some phonon modes only with the shell model.In recent years, a number of studies based on first-principles density functional theory (DFT) methods have been performed to reproduce the experimentally observed phonon dispersions of rutile TiO 2 .For instance, Lee et al [57] used variational density-functional perturbation theory within local-density approximation (LDA) to perform a full analysis of vibrational modes of rutile TiO 2 at the Γ point of the BZ, and found that the calculated phonon frequencies are in agreement with infrared and Raman spectroscopic measurements.Similarly, Sikora [58] used the direct method of ab initio lattice dynamics by employing LDA with ultrasoft pseudopotentials (PPs), and also found a good agreement with the experimentally measured phonon dispersion relations of rutile TiO 2 by neutron scattering.They further reported thermodynamic properties of rutile TiO 2 such as entropy, free energy, and heat capacity as a function of temperature, derived from the experimentally comparable phonon density of states (DOS).
On the other hand, Montanari and Harrison [59] performed a comparative DFT study of rutile TiO 2 using LDA, PW91 and PBE generalized gradient approximation (GGA) exchangecorrelational (XC)-functionals and the sensitivity of the lattice dynamics to the choice of XC-functional.They observed a discrepancy between the phonon frequencies predicted by LDA and GGA functionals, where LDA agrees with the experimental frequencies but GGA underestimates the measured frequencies.In particular, PBE generates imaginary transverseoptical A 2u phonon modes at the Γ point, leading to structural instabilities and to a ferroelectric transition which disagree with all previously reported literature.They attributed this discrepancy to the overestimation of the equilibrium volume by PBE GGA.Another similar study using LDA, PBE GGA, and with additional HSE06 hybrid functionals also show inconsistencies in the predicted phonon dispersion relations, where surprisingly most accurate hybrid functional HSE06 along with PBE GGA functional predict rutile TiO 2 as unstable structure [60].Although Fu et al previously obtained stable phonon dispersion relations in high-symmetry BZ paths for rutile TiO 2 using PBE GGA functional [61].
Meanwhile, both rutile and anatase TiO 2 were studied by LDA and PW91 GGA functionals, where PW91 underestimates the LDA predicted phonon bands for both phases of TiO 2 [62].This downshift of phonon frequencies subsequently affected the predicted heat capacity, where larger values are obtained by PW91 compared to LDA due to the larger occupation of states in lower frequencies regime with PW91.Another approach of studying the vibrational properties of both phases were reported in literature using PBE+U approach [63], where the choice of Hubbard U parameter affected the predicted phonon properties with larger U values leading to deviation from the experimental phonon frequencies and thermodynamic properties for both structures.The earliest first-principles work on lattice dynamics of anatase TiO 2 were reported by Mikami et al [64] using LDA XC-functional.They predicted stable phonon dispersion relations for anatase TiO 2 and evaluated dielectric tensor and Born effective charges consistent with experimental values.
In summary, LDA, PW91 and PBE GGA, PBE+U, HSE06 hybrid functionals have been used to predict the first-principles vibrational properties of rutile and anatase TiO 2 .Despite a great deal of theoretical investigation, only LDA showed reasonable experimental agreement with the predicted phonon dispersion relations, whereas PBE GGA showed contradictory results.Surprisingly higher level of theories used in XCfunctionals seemed to poorly estimate the phonon modes compared to experiment, but the underlying causes are not fully understood.Therefore, these studies should be revisited to properly understand the impact of XC-functionals on phonon properties of TiO 2 .As theoretical prediction of phonon modes in crystals within harmonic approximation is strongly influenced by a wide range of parameters, including the choice of a suitable XC-functional [65,66], adequate description of long-range bonding interactions and interpolation artifacts, and calculation convergence [67].Hence modeling TiO 2 with taking into consideration of all these parameters may improve the description of electronic structures and chemical bonding of TiO 2 , and thereby exploring the lattice dynamics of rutile and anatase TiO 2 to derive their accurate phonon properties.
In this study, we revisited the phonon properties of rutile and anatase TiO 2 within harmonic approximations using DFT methods within PBEsol GGA, TPSS meta-GGA, and HSE06 hybrid XC-functionals, which contradicts with the previously reported literature [60].Our findings suggest that stable phonon dispersion relations can be predicted for rutile and anatase TiO 2 within the harmonic approximation using all these XC-functionals.Calculated structural and thermodynamic properties of rutile and anatase TiO 2 are in good agreement with available experimental data.Corresponding vibrational DOS analysis showed the individual atomic contributions to the thermodynamic properties of both phases.A comparison of the thermodynamic properties such as entropy, specific heat, and relative molar enthalpy predicted by different levels of XC-functional theory showed that the thermodynamic properties of rutile and anatase TiO 2 are invariant on the level of DFT XC-functionals used.

Electronic structure methods
The experimental lattice structures of rutile (P4 2 /mnm) [35] (figure 1(a)), and anatase (I4 1 /amd) [36] (figure 1(b)) TiO 2 were obtained from the Inorganic Crystal Structure Database [68] for DFT calculations reported in this work.Non-spinpolarized DFT calculations were performed using Vienna Ab-initio Simulation Program (VASP) for the non-magnetic rutile and anatase TiO 2 [69].In order to find any potential dependency of TiO 2 lattice structures and vibrational properties on the choice of XC-functionals, PBEsol GGA, TPSS meta-GGA, and HSE06 hybrid XC-functionals were employed with the projector augmented wave (PAW) method [70].The valence configurations employed by these PAW potentials are specifically: 3d 3 4s 1 (Ti) and 2s 2 2p 4 (O).All calculations employed a 500 eV plane wave cut-off with accurate 'precision' mode in VASP, and added Gaussian smearing (with a width of 0.05 eV) to alleviate potential density convergence issues [71].The BZ of the rutile and anatase TiO 2 unit cells were sampled using 6 × 6 × 8 and 6 × 6 × 2 Monkhorst-Pack k-point meshes, respectively.The total DFT energies were fully optimized with respect to ionic positions, cell shape and cell volume, using ionic and electronic convergence criterion of ⩽|0.01| eV Å −1 and 1.0 × 10 −8 eV, respectively.

Elastic properties calculations
The elastic constants (C ij ) of tetragonal rutile (P4 2 /mnm) and anatase (I4 1 /amd) TiO 2 were calculated by strain-stress relationship with finite differences using VASP [72].Six finite distortions of the crystal lattice determined the elastic tensor and corresponding stress [73].The six independent elastic constants (C 11 , C 12 , C 13 , C 33 , C 44 , C 66 ) of tetragonal (I) determine the mechanical stability of rutile and anatase TiO 2 and are required to satisfy the following four necessary and sufficient Born elastic stability criteria [74], (2)

Vibrational properties calculations
First-principles phonon calculations of rutile (P4 2 /mnm) and anatase (I4 1 /amd) TiO 2 were performed by the finite displacement method [75].Vibrational properties of the optimized unit cells were calculated using the supercells approach with 2 × 2 × 4 and 3 × 3 × 1 expansion of the rutile and anatase unit cells, respectively, with the number of k points reduced accordingly.The force constants were calculated via self-consistent crystallographic supercells approach in terms of the Hellmann-Feynman forces induced by the atomic displacements of 0.015 Å in the supercells under PBEsol GGA, TPSS meta-GGA, and HSE06 hybrid functionals, as implemented in VASP [76,77].The non-analytical term corrections (NACs) were added to the dynamical matrices from the derived Born effective charges and dielectric constants that explicitly include the long-ranged anisotropic dipole-dipole interactions [77].Phonon modes (ω) were then calculated using these force constants via the Phonopy software package [54] using the Monkhorst-Pack grid of 48 × 48 × 48 q points for the phonon wave vectors.These phonon modes were then used to predict the phonon dispersion relations and phonon DOS of rutile and anatase TiO 2 .Once the phonon dispersion relations are defined over the BZ, the thermal properties such as the temperature-dependent entropy (S), constant volume heat capacity (C v ) and Helmholtz free energy (A) are computed within the harmonic approximation using the thermodynamic relations at a temperature T and constant volume V [78], Here the summations run over wave vectors q and band indices j.On the other hand, LDA [79] reported a and c parameters differed by ∼0.9 and 1.2%, respectively, while PW91 [79] and From the comparison of predicted C ij values by all employed functionals, it showed that PBEsol generated close to experimental elastic moduli compared to the other functionals, while HSE06 showed highest deviation from experiment.Contrary, predicted C ij values by TPSS have less deviation than HSE06 but more than PBEsol.Therefore, the deviation of theoretical elastic moduli from experiment increased for the higherrung XC-functionals tested.Notably, the reported C ij values by LDA, PW91 and PBE GGA are found to be largely underestimated the experimental values, as well with the predicted values by PBEsol, TPSS and HSE06 functionals in the present work.However, C ij values calculated by all XC-functionals met Born stability criteria in equations ( 1) and (2) [74], indicating that all employed functionals predicted rutile TiO 2 as a mechanically stable structure.

Results and discussion
Calculated bulk modulus, Young's modulus and shear modulus of rutile TiO 2 (P4 2 /mnm) are presented in table 3. PBEsol  PW91 showed slightly poor performances in predicting these properties.

Phonon dispersion relations.
The phonon dispersion relations of rutile TiO 2 (P4 2 /mnm) predicted by PBEsol, TPSS and HSE06 functionals are compared with the experimental phonon dispersion relations along the high symmetry BZ paths of Γ-X, Γ-M, and Γ-Z in figure 2. Overall, TPSS meta-GGA showed excellent agreement with the experimentally measured acoustic phonon branches at these selected BZ paths within low-frequency regimes below 300 cm −1 , and with the optical phonon branches at intermediate frequency regimes between 400 and 450 cm −1 .Although HSE06 has better agreement with optical phonon branches at high-frequency regimes above 750 cm −1 .Notably, HSE06 showed reasonable agreement with the experimental optical phonon branches at intermediate frequency regimes between 400 and 450 cm −1 , and PBEsol prediction agreed quite well at 450 cm −1 between Γ-X and Γ-M paths.Shojaee and Mohammadizadeh [62] predicted stable phonon dispersion relations of rutile TiO 2 using LDA and PW91 GGA functionals, which agreed with the prediction of PBE GGA functional by Fu et al [61].However, Wu et al [60] lately reported imaginary phonon modes along the acoustic branch around R-point of rutile TiO 2 by PBE GGA and HSE06 hybrid functionals, and identified this structure as unstable.Notably, they also obtained stable phonon dispersion curves at R-point of rutile TiO 2 using LDA functional.This contradictory prediction makes it questionable whether the calculated phonon modes by GGA and HSE06 at R-point has any physical significance or it reflects any technical problems in the calculations (e.g.convergence issues, size of supercell).Therefore, PBEsol GGA, TPSS meta-GGA, and HSE06 hybrid functionals have been employed in this work to predict the phonon dispersion relations of rutile TiO 2 along the high symmetry BZ paths of Γ-X-M-Γ-Z-R-A-Z as shown in figure 3. Interestingly, all employed functionals predicted stable phonon dispersion curves for rutile TiO 2 , specifically no soft modes were observed at R-point which disagreed with the previously reported literature [60] but agreed with other theoretical works [61,62].This discrepancy might arise from the unconverged supercell size used in the previous literature [60] (for instance, in this present work, we observed unstable phonon modes at R-point with 2 × 2 × 2 and 2 × 2 × 3 expansion of rutile unit cell for specific XC-functionals, however, 2 × 2 × 4 and 2 × 2 × 5 expansion generated stable phonon modes in all considered BZ paths, as shown in figures S1-S5).However, the level of theory employed in the XC-functionals in the present work impacted the distribution of phonon dispersion along the BZ path.For instance, the predicted phonon dispersion curves by TPSS are at higher energy levels compared to the prediction by both PBEsol and HSE06 in the lowfrequency regime of 0-200 cm −1 .On the other hand, TPSS showed reasonable agreement with PBEsol and HSE06 in the intermediate-frequency regime (450-600 cm −1 ), but showed inconsistency with HSE06 in high-frequency regimes (600-850 cm −1 ).TPSS has good agreement with PBEsol's prediction in high-frequency regimes, while at low-frequencies PBEsol agreed with the HSE06 prediction of phonon dispersion.In summary, PBEsol GGA, TPSS meta-GGA, and HSE06 hybrid functionals confirmed rutile TiO 2 as a stable structure with no imaginary phonon modes at high-symmetry points, however, the distribution of phonon bands along the BZ paths are impacted by the DFT XC-functionals used.

Phonon DOS.
The influence of XC-functionals on the individual atomic contributions to the phonon dispersion of rutile TiO 2 (P4 2 /mnm) has been examined.The projected  phonon DOS of Ti 4+ and O 2− in rutile TiO 2 were calculated by PBEsol, TPSS and HSE06 functionals and are shown in figure 4 (with NAC) and figure S6 (without NAC).HSE06 hybrid functionals expanded the total phonon DOS over the frequency ranges (0-850 cm −1 ), while PBEsol and TPSS limited the DOS range within 0-820 cm −1 .The projected phonon DOS can be categorized into three different regimes: low-frequency (0-160 cm −1 ), intermediate-frequency (160-600 cm −1 ), and high-frequency (600-850 cm −1 ).The range of phonon frequencies associated with the Ti 4+ in rutile TiO 2 occurred between low-and intermediate-frequency regimes, where the principle vibrational contributions by Ti 4+ were observed within low-frequency between 0 and 160 cm −1 .The principle vibrational mode of Ti 4+ predicted consistently by HSE06 and PBEsol was observed at ∼90 cm −1 , whereas TPSS shifted it to the higher energy level near ∼120 cm −1 .Consistent with the phonon dispersion curves, the phonon frequency of Ti 4+ shifted around 40 cm −1 higher by TPSS at low-frequencies compared to PBEsol and HSE06.The overall amplitude of Ti 4+ vibrational DOS consistently decreased moving into the intermediate-and high-frequency regimes for all functionals.Conversely, the phonon DOS associated with the O 2− vibrations in the rutile TiO 2 were more pronounced within the intermediate frequency range, with a peak amplitude at ∼460 cm −1 for both TPSS and HSE06, and at ∼450 cm −1 with PBEsol functional.A secondary peak associated with the O 2− vibrations was observed near ∼480 cm −1 for all XC-functionals.
Notably, TPSS predicted the peak amplitude of Ti 4+ vibrational frequencies at higher energy level in the low-frequency regimes compared to other XC-functionals, but in case of O 2− the peak was observed at a similar energy level to that of the HSE06 functional.In high frequency regimes, HSE06 lifted the O 2− frequencies around 40 cm −1 higher than PBEsol and TPSS values, which were also observed in their respective optical phonon bands.To sum up, TPSS showed deviation from the PBEsol and HSE06 predicted vibrational contributions of Ti 4+ at the low-frequency regimes, as a result acoustic phonon branches predicted by TPSS were observed at higher energy level compared to PBEsol and HSE06 (figure 3).However, the vibrational contributions of O 2− predicted by TPSS agreed with HSE06 at the intermediate frequency regimes and as a consequence, TPSS and HSE06 predicted optical phonon branches were observed at the similar energy level in this regime.On the other hand, since the vibrational contributions of O 2− predicted by TPSS and PBEsol at the high-frequency regimes agreed with each other, their optical dispersion curves were also found at the same energy level in this regime.

Calculated vibrational entropy (S), specific heat (C v ) and relative molar enthalpy [H(T)-H(298 K)
] of rutile TiO 2 (P4 2 /mnm) by PBEsol, TPSS, and HSE06 functionals are compared between 0 and 2000 K with experimental data [92] in figures 5(a)-(c) with NAC, and figures S7(a)-(c) without NAC, respectively.The calculated S values by different XCfunctionals are in good agreement with experimental values between 298 and 2000 K (figure 5(a)).For instance, the experimental S value at room temperature is reported around 49.83 J K −1 mol −1 .At equivalent temperature, TPSS predicted the S value close to experiment at 50.75 J K −1 mol −1 .Comparatively, HSE06 overestimated the experimental value by 4 J K −1 mol −1 , and PBEsol overvalued by 5 J K −1 mol −1 .The S values increased almost exponentially with all functionals with increasing temperature up to 2000 K.However, the predicted S values by HSE06 at high temperature showed significant improvement compared to their prediction at low temperature.For example, HSE06 predicted an identical experimental S value of ∼185.5 J K −1 mol −1 at 2000 K.In contrast, TPSS underestimated this value by 3 J K −1 mol −1 and PBEsol overestimated by 2 J K −1 mol −1 .
Predicted C v values of rutile TiO 2 (P4 2 /mnm) by PBEsol, TPSS, and HSE06 are consistent with each other between the temperature range 0-2000 K in figure 5(b).The variation of C v value with temperature obeyed the Debye's T 3 law at low temperature range (0-300 K).At 298.15 K, the experimental C p value is 55.18 J K −1 mol −1 , while theoretical C v values calculated by HSE06 (55.68 J K −1 mol −1 ) and TPSS (55.99 J K −1 mol −1 ) are almost identical with the experimental value at equivalent temperature.Conversely, PBEsol overestimated the experimental value by 1 J K −1 mol −1 .Above 500 K, the theoretical C v values deviated from experimental C p due to the classical limit of Dulong-Petit's law, indicating theoretical C v values predicted by different XCfunctionals provided a reliable and consistent measurement.At high temperatures the deviation between the predicted C v values by different XC-functionals became narrower.For example, PBEsol, TPSS, and HSE06 predicted more or less identical C v values ∼74 J K −1 mol −1 at 2000 K, a 3 J K −1 mol −1 less than experimental C p value.The predicted [H(T)-H(298 K)] values of rutile TiO 2 (P4 2 /mnm) by PBEsol, TPSS, and HSE06 are shown in figure 5(c).The calculated [H(T)-H(298 K)] values varied almost linearly with temperature and showed excellent experimental agreement between 298 and 2000 K.The theoretical [H(T)-H(298 K)] values by different XC-functionals are effectively invariant at all temperature range considered here.For instance, at 300 K, the predicted [H(T)-H(298 K)] value is ∼0.1 kJ mol −1 by PBEsol, TPSS and HSE06, reached close to the experimental value of 0.09 kJ mol −1 , while at 2000 K all XC-functionals predicted this value ∼121.5 kJ mol −1 , a 4 kJ mol −1 less than the experimental value.In summary, the predicted S, C v and H(T)-H(298 K)] values of rutile TiO 2 (P4 2 /mnm) by PBEsol, TPSS, and HSE06 suggested that the thermodynamic properties of rutile TiO 2 has little, if any, dependence on the choice of level of theory used in DFT XC-functionals.
3.2.Anatase (I4 1 /amd) TiO 2 3.2.1.Lattice structure and elastic stability.PBEsol, TPSS and HSE06 predicted optimized lattice parameters of anatase (I4 1 /amd) TiO 2 are presented in table 4 and compared with available theoretical and experimental data.PBEsol and HSE06 underestimated the experimental lattice vector (a/Å) by ∼0.8% and ∼0.9%, respectively, while TPSS only underestimated by less than 0.1%.Interestingly, HSE06 overestimated the other experimental lattice vector (c) by ∼0.5%, while PBEsol and TPSS underestimated the experimental c parameter by ∼0.9 and ∼1%, respectively.The resulting lattice parameters by PBEsol, TPSS, and HSE06 compressed the four shorter and two longer Ti-O bond lengths of TiO 6 octahedra in anatase TiO 2 .For example, compared to the experimentally measured shorter Ti-O bond lengths of 1.9351 Å (4×) and longer bond lengths of 1.9866 Å (2×), PBEsol reduced these values to 1.9194 Å and 1.9714 Å, respectively.Relatively, TPSS generated closer to the experimental bonds lengths with a value of 1.9314 Å and 1.9768 Å, respectively, for shorter and longer Ti-O.HSE06 moderately decreased the shorter bond lengths to 1.9239 Å, and longer bond lengths to 1.9741 Å.In general, all three XC functional displayed an overbinding tendency of anatase TiO 2 , similar to the rutile one.
As a consequence, HSE06 underestimated the experimental lattice volume only by ∼1.2%, whereas the underestimation of a parameter and overestimation of c parameter minimized the calculated volume discrepancy by HSE06 from the experiment.Comparatively, TPSS produced the closest experimental lattice volume, with a least deviation of ∼1% than experimental value.However, PBEsol largely underestimated the lattice volume of anatase TiO 2 by ∼2.5% due to the underestimation of both a and c lattice vectors.Comparably, LDA [79] reported lattice volume underestimated by ∼2.7% from experiment, while PW91 [79] overestimated by ∼2.5%, however, PBE [93] is found to produce close experimental value with a deviation of just ∼0.9%.All of these optimized lattice parameters by different XC-functionals suggested that PBE and TPSS are reasonable and economical to predict accurate lattice parameters of anatase TiO 2 .
The elastic parameters of anatase TiO 2 determine its structural stability.Therefore, the elastic moduli of anatase TiO 2 is computed by PBEsol, TPSS, and HSE06 functionals (table 5) and assessed the impact of different XC-functionals (including LDA and GGA from literature) on these properties.Similar to rutile, anatase TiO 2 also has the tetragonal form, so it has six unique elastic moduli (C ij ).The deviation between present work and other theoretical values of six C ij for anatase TiO 2 varied from functional to functional.For instance, TPSS generated C 12 , C 33 and C 66 values close to LDA [79] compared to the predicted values by PBEsol and HSE06, whereas both PW91 [79] and PBE [93] GGA showed reasonable agreement with each other for all predicted C ij values.In contrast, HSE06 largely overestimated all C ij values (except C 11 and C 44 ) with respect to all other XC-functionals, consistent with the prediction of elastic moduli for rutile TiO 2 .Notably, irrespective of their deviation, all XC-functionals generated standard elastic moduli that satisfied the Born stability criteria in equations ( 1) and ( 2) [74], indicating anatase as a stable TiO 2 structure.PBEsol, TPSS, and HSE06 predicted bulk modulus, Young's modulus and shear modulus of anatase (I4 1 /amd) TiO 2 , are compared with available theoretical and experimental values in table 6. TPSS showed dominant performance in predicting experimental bulk modulus compared to PBEsol and HSE06.For instance, TPSS predicted bulk modulus of anatase TiO 2 to be 188 GPa, only 9 GPa less than experiment.Contrary, PBEsol overestimated this value by 17 GPa, while HSE06 overestimated by 30 GPa.No experimental data is found to compare the Young's and shear modulus of anatase TiO 2 .It can be seen that the predicted Young's and shear modulus by PBEsol, TPSS and HSE06 functionals are consistent with each other, while calculated shear modulus have good agreement with the previous predictions by PW91 [79] and PBE [93] GGA functionals, but largely differed from the reported LDA [79] value.Similarly, LDA calculations provided larger value for Young's modulus of anatase TiO 2 than present (PBEsol, TPSS, HSE06) and previous calculations (PW91, PBE).

Phonon dispersion relations.
PBEsol, TPSS, and HSE06 predicted phonon dispersion curves of anatase (I4 1 /amd) TiO 2 have been examined along the high symmetry BZ paths Γ-X-P-N-Γ-M-S, with and without NAC, and are presented, respectively, in figures 6 and S8.A negligible soft mode (<5 cm −1 ) was observed between Γ-X with NAC (figure 6), which is far below the threshold to be considered as imaginary phonon frequencies [96].However, stable phonon modes were obtained without NAC in all BZ paths by PBEsol, TPSS, and HSE06 functionals (figure S8), consistent with the previous observation by LDA and PW91 GGA functionals [40,62].At a very low-frequency range (0-150 cm −1 ), all functionals produced almost similar acoustic phonon frequencies over the selected BZ paths.Above 150 cm −1 , PBEsol started to underestimate the phonon bands from other two functionals up to 650 cm −1 .Over 650 cm −1 , PBEsol optical phonon frequencies were observed at identical energy level with TPSS phonon frequencies along the high-symmetry BZ paths.
In contrast, HSE06 generated optical phonon bands are at marginally higher energy levels compared to PBEsol and TPSS in the frequency ranges 650-875 cm −1 , with the gap being more pronounced from intermediate-to highfrequency regions (500-875 cm −1 ).As observed previously for rutile TiO 2 , TPSS showed variable performance in predicting phonon bands over the low-, intermediate-and highfrequency regimes.TPSS agreed with HSE06 and PBEsol at low-frequency regimes below 100 cm −1 , but largely overestimated them at intermediate regions between 400 and 500 cm −1 , and agreed only with PBEsol at high-frequency regions above 650 cm −1 .

Phonon DOS.
In order to analyze the discrepancies observed in the phonon dispersion curves of anatase TiO 2 predicted by different XC-functionals, the individual atomic contributions to the phonon DOS have been illustrated in figure 7 with NAC, and in figure S9 without NAC.The phonon DOS associated with the vibration of Ti 4+ in anatase TiO 2 were found to be more evident between 0 and 200 cm −1 , with a principle vibrational mode at ∼165 and ∼175 cm −1 , respectively, for PBEsol and TPSS (figure 7).Unlike the other two functionals, HSE06 produced the principal mode at a lower energy level of around 140 cm −1 , while still having the highest peak amplitude (figure 7).The vibrational modes of Ti 4+ became weaker thereafter, as it continued up to the energy level ∼550 cm −1 .Concurrently, the DOS associated with O 2− vibrations in the anatase TiO 2 were visible between 50 and 870 cm −1 .The peak amplitude of O 2− vibrational modes were observed at ∼430 cm −1 with TPSS, whereas HSE06 shifted the peak amplitude at ∼630 cm −1 .Contrary, PBEsol generated the peak amplitude of O 2− vibrational modes at slightly lower energy level than HSE06 at ∼615 cm −1 but with the highest amplitude than other two functionals.A secondary and tertiary peak amplitudes were observed at ∼510 and ∼630 cm −1 , respectively by TPSS, whereas the same peaks were recognized by PBEsol at ∼505 and ∼420 cm −1 , respectively, and by HSE06 at ∼525 and ∼430 cm −1 , respectively.The mutual contributions of Ti 4+ and O 2− vibrations between the frequency ranges of 100-525 cm −1 mostly reflected the discrepancies observed in the phonon dispersion curves (figure 6) by different XC-functionals.For instance, the maximum vibrational DOS amplitudes of Ti 4+ and O 2− between the frequency regime 400-525 cm −1 were realized by TPSS, overpopulating phonon bands in that particular energy level.On the other hand, HSE06 was found to shift the O 2− vibrations to the higher energy level close to 870 cm −1 , which ultimately led to the distribution of optical phonon bands at upward energy level compared to TPSS and PBEsol.For example, between the temperature ranges 0-100 K, the increase in S was nearly 10 J K −1 mol −1 predicted by all XC-functionals, whereas between 100 and 200 K it yielded additional 20 J K −1 mol −1 in S. At 298 K, the S value predicted by HSE06 is 50.17J K −1 mol −1 , agreeing with the experiment value of ∼49.82 J K −1 mol −1 .At an equivalent temperature, TPSS predicted the near exact value of 50.59 J K −1 mol −1 , while PBEsol slightly overestimated this value (∼52.54J K −1 mol −1 ).HSE06 generated 130 J K −1 mol −1 of S value at a temperature 1000 K, only 2 J K −1 mol −1 less than the experimental value.On the other hand, TPSS and PBEsol predicted this value as ∼131 and ∼133 J K −1 mol −1 , respectively.Interestingly, PBEsol reached close to the experimental value of S (∼185.5 J K −1 mol −1 ) at 2000 K by generating an S value of 184 J K −1 mol −1 .Surprisingly, at the same temperatures, the predicted S by HSE06 and TPSS are equally 4 J K −1 mol −1 less than experiment.
The variation of C v with temperature is shown in figure 8(b).The trends in C v of anatase TiO 2 are the same as those observed for rutile.C v obeyed Debye's T 3 law at low temperature below 300 K, and respected the classical limit above 500 K.At room temperature TPSS and HSE06 produced identical experimental C v of 55 J K −1 mol −1 .PBEsol overestimated this value by only 1 J K −1 mol −1 .The discrepancy between the predicted values by different XC-functionals diminishes at high temperature.For instance, all used functionals here predicted almost identical C v of ∼72.5 ± 0.2 and 74 J K −1 mol −1 at 1000 K and 2000 K, respectively.Linear variation with temperature was observed in the predicted [H(T)-H(298 K)] values of anatase TiO 2 as shown in figure 8(c).At 300 K, all XC-functionals generated [H(T)-H(298 K)] around 0.1 kJ mol −1 , close to the experimental value of ∼0.08 kJ mol −1 .A sharp increase of 38 kJ mol −1 was observed between 300 and 1000 K by all XC-functionals, which differed only by 1 kJ mol −1 from the experimental value.At 2000 K, the theoretical [H(T)-H(298 K)] reached to ∼121 kJ mol −1 with a deviation of 5 kJ mol −1 with experiment.Overall, the theoretical S, C v , and H(T)-H(298 K)] values of anatase TiO 2 obtained by PBEsol, TPSS, and HSE06 showed consistency with each other and agreed with experimental data, despite the differences observed in the calculated phonon band structures for these XC-functionals.

Conclusion
The structural stability, phonon and thermodynamic properties of rutile (P4 2 /mnm) and anatase (I4 1 /amd) TiO 2 were reported using first-principles DFT methods with the PBEsol GGA, TPSS meta-GGA, HSE06 hybrid XC-functionals.Theoretical lattice parameters predicted with these functionals were in good agreement with the experimental values.The elastic modulus, bulk modulus, Young's modulus and shear modulus predicted by PBEsol and TPSS were close to the experimental values, however, the hybrid functional HSE06 almost always overestimated these values.All applied functionals predicted rutile as a stable structure with and without NAC and no imaginary phonon frequencies were observed along the high-symmetry BZ paths by any of these applied functionals.However, the phonon dispersion curves of anatase TiO 2 with NAC were found to be slightly unstable along the Γ-X path in the BZ, although without NAC it showed stability in all high-symmetry points of the BZ.Nevertheless, the predicted phonon bands over the BZ were impacted by different XC-functionals used.The predicted phonon frequencies of rutile TiO 2 between 0 and 160 cm −1 were at higher energy levels when calculated with TPSS compared to PBEsol and HSE06, originating from the vibrational contributions of Ti 4+ .However, TPSS showed reasonable agreement with hybrid functionals between 450-600 cm −1 , where O 2− ions have substantial vibrational contributions along with Ti 4+ .Notably, TPSS showed very good agreement with the optical phonon bands predicted by PBEsol between 650 and 800 cm −1 , where only O 2− ions contribute.In contrast, for anatase TiO 2 , all tested functionals showed reasonable agreement in generating phonon bands at low-frequency regimes below 100 cm −1 .Still, TPSS overestimated the phonon frequencies compared to all other functionals at intermediate regions between 400 and 500 cm −1 , due to the net vibrational modes of Ti 4+ and O 2− .Overall, HSE06 shifted the optical phonon bands of both rutile and anatase TiO 2 to higher energy levels than PBEsol and TPSS, due to the shifting of O 2− vibrations.The temperaturedependent thermodynamic properties of rutile (P4 2 /mnm) and anatase (I4 1 /amd) TiO 2 predicted by these XC-functionals are in excellent experimental agreement and thermodynamic properties are not appreciably influenced by the level of theories employed to the applied DFT XC-functionals.

Figure 1 .
Figure 1.Crystal structures of TiO 2 polymorphs: (a) rutile (P4 2 /mnm) and (b) anatase (I4 1 /amd).The light blue and red spheres represent Ti and O atoms, respectively.The light blue surfaces point to the TiO 6 octahedra.Solid square black lines in both structures indicate the individual unit cell.

Table 1 .
Comparison of optimized PBEsol, TPSS and HSE06 lattice parameters (a and c), and volume (V) of rutile (P4 2 /mnm) with experimental and other theoretical values.
[35]ta from[79].bDatafrom[80].cThiswork.dExp.datafrom[35].predicted by PBEsol, TPSS and HSE06 functionals, and are compared with available theoretical and experimental values in table1.The optimized lattice vectors (a/Å and c/Å) computed by PBEsol, TPPS, and HSE06 slightly underestimated the experimental values.For instance, the computed a and c values by TPSS are close to the experiment with a deviation of just ∼0.2% and ∼0.3%, respectively, while PBEsol calculated values of a and c have the maximum deviation of ∼0.6% and ∼0.9%, respectively.HSE06 predicted a and c parameters differed by only ∼0.5% and ∼0.6%, respectively, to the experimental values.These marginal underestimation of lattice parameters by PBEsol, TPSS, and bond lengths to 1.9468 Å, and the longer one to 1.9693 Å. HSE06 lowered the bond lengths to 1.937 Å and 1.9719 Å, respectively, for the shorter and longer Ti-O bonds.Overall, PBEsol, TPSS, and HSE06 showed overbinding of rutile TiO 2 .

Table 2 .
Comparison of PBEsol, TPSS and HSE06 predicted elastic constants (C ij ) of rutile TiO 2 (P4 2 /mnm) with experimental and other theoretical values.The structural stability of rutile TiO 2 is assessed in terms of elastic behavior.As rutile TiO 2 has the tetragonal form with space group of P4 2 /mnm, it has six unique elastic moduli (C ij ).These C ij values of rutile TiO 2 are predicted by PBEsol, TPSS, and HSE06 functionals and are compared with their corresponding experimental and theoretical values in table 2.
[81]ta from[79].bDatafrom[80].cThiswork.dExp.datameasuredat 4 K, from[82].eExp.datameasuredat80K,from[83].PBE[80]GGA both overestimated the a and c parameters of rutile TiO 2 as shown in table 1.Overall, PBEsol, TPSS, and HSE06 estimated the lattice parameters of rutile TiO 2 with less than 1% deviation of experimentally reported values by Burdett et al[35].As a result, the lattice volume predicted by TPSS has the least deviation with experiment (∼0.8%), while PBEsol generated the highest deviation (∼2%).The HSE06 hybrid functional produced a deviation of ∼1.5% from the experiment lattice volume.Contrary, LDA underestimated the lattice volume by ∼3%, while PW91 and PBE both largely overestimated the lattice volume of rutile TiO 2 by ∼2.3 and ∼5%, respectively, as reported in literature.This result suggested that PBEsol, TPSS, and HSE06 hybrid functional are competitive in predicting experimental lattice parameters of rutile TiO 2 .It should be noted, however, that incorporating zero-point energy effects in the quasiharmonic approximation may increase the calculated volumes slightly[81].

Table 3 .
Comparison of PBEsol, TPSS, and HSE06 predicted bulk modulus (B/GPa), Young's modulus (E/GPa) and shear modulus (G/GPa) of rutile TiO 2 (P4 2 /mnm) with experimental and other theoretical values.TPSS predicted bulk modulus are close to the experimental and other theoretical values, whereas calculated values by HSE06 are within the experimental error range.Similarly, the computed Young's modulus of rutile TiO 2 by PBEsol and TPSS differed by only 8 and 6 GPa, respectively, from the experimental value, whereas this difference is more pronounced in the case of HSE06.On the other hand, TPSS predicted the exact experimental shear modulus of rutile TiO 2 , while PBEsol underestimated this value by only 1 GPa.HSE06 overestimated the shear modulus by 15 GPa compared to the experimental reported value.Overall, PBEsol and TPSS showed reasonable agreement in predicting bulk, [84][85][86][87]bDatafrom[80].cDatafrom[88].dThiswork.eExp.datafrom[89].fExp.datafrom[90].gExp.datafrom[91].hExp.datafrom[17].andYoung's, and shear modulus of rutile TiO 2 , while HSE06 largely overestimated these values.This is in line with the observation of other solid materials using HSE06 functional[84][85][86][87].Interestingly, the reported bulk, Young's, and shear modulus by other theoretical works showed that, LDA and PBE prediction have better agreement with experiment, while