A DFT + U study on structural, electronic, vibrational and thermodynamic properties of TiO2 polymorphs and hydrogen titanate: tuning the Hubbard ‘U-term’

Structural, electronic, vibrational and thermodynamic properties have been tested when Hubbard parameter U is implemented in density functional theory calculations for TiO2 polymorphs: anatase, rutile, TiO2–B and for hydrogen titanate (H2Ti3O7) bulk. Optimum U parameter values were found for each system, balancing geometric changes and electronic properties, namely, U = 4 eV for anatase and TiO2–B, U = 5 eV for rutile and hydrogen titanate. Although the addition of this parameter improves the prediction of electronic properties, with no significant structural changes, we found that it would not be adequate for predicting vibrational properties.


Introduction
Titanium oxide (TiO 2 ), has been widely studied both experimentally and theoretically [1][2][3][4][5] due to its potential applications like high-efficiency solar cells, photocatalysts, and storage capacitors in dynamic random access memories [6][7][8]. TiO 2 presents various polymorphs among which the most studied ones are anatase and rutile. They exhibit exceptional physical and chemical properties, specially for photochemical and photoelectrochemical applications. Sunlight can be used to produce electricity or to drive chemical reactions by photocatalysis and photovoltaics. In order to find efficient materials for solar energy conversion we need to study their properties i.e. structure, surfaces and interfaces, light absorption, charge transport, electron and hole trapping, among others.
In systems with d-and f-localized electrons the overestimation of electron delocalization is a known weakness in DFT methodology. In the 90 s the DFT+U method was developed, which consists in an explicit treatment of electronic correlation with a Hubbard-like model for a subset of states in the system [15]. The essential idea is to treat the strong on-site Coulomb interaction of localized electrons using an additional Hubbard-like term. The strength of the on-site interaction is described by two parameters U and J, the on-site Coulomb and the on-site exchange interaction, respectively. The DFT+U method has been widely applied in several systems [16][17][18][19][20]. This implementation improves the results in the calculation of energetic, electronic and magnetic properties of insulating and semi-conducting materials which contain transition metals. Nevertheless, in early transition metal compounds (like Ti) the more extended orbitals decrease the electron correlation, so it is needed to study the performances of the DFT+U method. In the last decade, several studies have considered the U-term in predicting the properties of anatase and rutile [4,[21][22][23].
Here, we report how the structural, electronic, vibrational and thermodynamic properties behave when the U-term is taken into account during calculation. As it is known, theoretical modeling plays a very important role Original content from this work may be used under the terms of the Creative Commons Attribution 3.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.
in materials research, since it provides information at an atomic level which is inaccessible via experiments, helping to identify which factors control the behavior of a specific material.

Methodology
Electronic structure and energy calculations have been performed within the frame of the density functional theory (DFT) [24] implemented in the Vienna ab initio simulation package (VASP) [25]. The projector augmented wave pseudopotential [26] was used to account for the electron−ion core interaction, using the Perdew-Burke-Ernzerhof functional as generalized gradient approximation (GGA) [27] for the exchange −correlation term. A cutoff energy of 400 eV was used to expand the Kohn-Sham orbitals into plane wave basis sets. The employed pseudopotentials correspond to the following configuration: 3s 2 3p 6 3d 2 4s 2 for Ti, 2s 2 2p 4 for O and 1s 1 for H. This choice amounts to a valence of 12 for Ti, 6 for O and 1 for H. In every case the K-point mesh taken was equivalent to 4×4×4 Monkhorst-Pack grid [28] for the full (reducible) Brillouin zone, allowing total convergence. We have implemented the Hubbard U model to improve the calculated band gap width [15,[29][30][31]. In this implementation, the on-site coulombic (U) and exchange (J) terms are combined into a single effective U parameter (U eff ) to account for errors in exchange correlation on Ti 3d orbitals.
We have modeled the following TiO 2 polymorphs: anatase (unit cell containing 12 atoms), rutile (6 atoms) and TiO 2 -B (24 atoms), as well as hydrogen titanate (21 atoms), as it can be seen in figure 1 (for better visualization the unit cell was multiplied in some directions).
As a first step the bulk structures (i.e. the cell parameters and atomic positions) were fully relaxed in the absence of the U parameter (U=0 eV). Then, from these relaxed structures, the bulks were allowed to relax again, but within the GGA+U approximation, setting the U term which takes the following values: U=3, 4 and 5 eV.
In order to proceed with the vibrational and thermodynamics analysis, we performed first a phonon analysis by means of density functional perturbation theory as implemented in VASP. For these calculations, we increased the accuracy of the grids calculations by selecting a better grid for the FFT-grid and the fine FFT-grid, corresponding to 2×G cut and 4×G cut , correspondingly (PREC=accurate and ADDGRID). Additionally, we doubled the k-point sampling in each direction, assuring an 8×8×8 grid sampling. It guarantees a better accuracy for phonons and second derivatives calculations. With this purpose, we have utilized the pristine unit cell (see table 1) for all the considered polymorphs. Due to possible long range dipole-dipole interactions, the non-analytical term correction was considered in the calculation. In every case the absence of soft modes were checked by the analysis of the phonon dispersion diagrams. The phonon analysis and the further thermodynamical properties were determined by using the post-processing code Phonopy [32].

Results and discussion
Geometric optimization As it was said above, anatase, rutile, TiO 2 -B and hydrogen titanate were relaxed in order to find their optimum geometries. In table 1 the lattice parameters a, c and c/a ratio at different U values are listed.
The incorporation of U term affects the lattice parameters, expanding the cell while the U value increases. That is the reason why it is needed to find the right U value, which balances the accuracy of the predicted geometrical structure and electronic properties for each system. In table 2 it can be seen that Ti-O bond lengths in every system increase when U takes higher values. However, the changes are not significant or drastic. Blue, red and white spheres represent titanium, oxygen and hydrogen atoms, respectively.

Electronic properties
In table 3 a comparison of the band gap widths at different U term values and experimental results obtained from literature are summarized [33][34][35]. As the U term value increases, the band gap increases as well.
The density of states is plotted in figure 2, each system at promising optimum U term value, namely, U=4 eV for anatase and TiO 2 -B, U=5 eV for rutile and hydrogen titanate. Band gap energies are in good agreement with those found in previously cited literature for TiO 2 -B and hydrogen titanate. However for the cases of anatase and rutile a selection of higher U values leads to a major structural deviation, that is why we have selected values of U=4 eV for anatase and U=5 eV for rutile (see table 4). Despite the band gap energies are far from experimental reported data an improvement can be seen when U term in taken into account in the calculation.
In these curves a significant contribution of O 2p states to the valence band, and an important contribution of Ti 3d states in the conduction band, can be seen. This is a well-known feature for bulk TiO 2 polymorphs. In order to determine the appropriate values of U for each system, we have compared the calculated band gap widths for a set of U values (3, 4 and 5 eV) with experimental band gap values. However, it is not only necessary to fit the band gap width but also it is needed to maintain the experimental lattice parameters. This information can be seen in table 4, for TiO 2 -B polymorph and hydrogen titanate there is no structural experimental information available and for this reason the values calculated for different U parameters were compared with the results calculated using U=0 eV.

Vibrational properties
In order to obtain the zone-center phonons of the TiO 2 polymorphs and hydrogen titanate, the dynamical matrix should be diagonalized for q=0. Through factor analysis [36], in the I4 1 amd, P4 2 /mnm, C12/m1 and C2/m space groups for anatase, rutile, TiO 2 -B and hydrogen titanate bulks respectively, the irreducible representations of the optical vibrations are as follows: The modes with subscript g (gerade) are Raman-active and those with subscript u (ungerade) are infraredactive, while the modes represented with the letter E are degenerate.
In all the cases, the phonons density of states have been obtained, and Ti, O and H atoms are projected in black, red and blue curves, respectively, as it can be seen in figure 3. In TiO 2 polymorphs the low values of frequency correspond mainly to titanium atoms, and high values to oxygen atoms, anatase and rutile plots agree with those from literature [37]. In the case of hydrogen titanate we can see a wide gap among curves, there are two picks at high energy values which correspond to the light hydrogen atoms.
In figure 4 the histograms of each bulk system are presented, the bars represent the group of vibration modes in a specific range of frequency. The vertical black lines are the vibration modes calculated at DFT U=0 eV level for comparison. From these charts, it can be noticed that the implementation of the U term has an effect on the vibrational properties, developing different vibration modes at different frequencies. How this influences the thermodynamic properties will be analyzed in the following section.
For bulk rutile, there are lots of data reported: in table 5 a comparison of this information and our results is shown [38][39][40][41][42]. Here it can be seen clearly how the Hubbard implementation changes the vibrational properties, frequency values of vibrational modes differ and shift significantly when different values of U term are taken into account in the calculation. Reported theoretical results without U term approach better to experimental information.

Thermodynamic properties
Once the phonon frequencies set is obtained, thermodynamic properties of the systems can be calculated by quasi harmonic approximation. According to the following expressions, entropy, vibration energy and Helmholtz free energy can be determined: Also, specific heat at constant volume can be calculated as follows: Properties of solids can be extrapolated at temperatures higher than 0 K by vibration harmonic approximations. Based on these quasi-harmonic approximations, we have obtained the phonon frequencies at five different volumes, namely, 3% and 1.5% reduced volume, 0%, and 3% and 1.5% expanded volume. In this manner, we have calculated a set of thermodynamic parameters for TiO 2 polymorphs and hydrogen titanate bulks. In figures 5-8 we have plotted enthalpy and entropy versus temperature for the four systems as well as C p (specific heat at constant pressure) and C v (specific heat at constant volume) versus temperature, respectively. In black crosses experimental data for anatase and rutile from literature was added [43].
It can be seen in these four figures that the three TiO 2 polymorphs behave in a similar way, while hydrogen titanate takes higher values in every thermodynamic property studied here. In the enthalpy versus temperature plot (figure 5), we can see that up to 200 K the three TiO 2 polymorphs take very similar values at standard DFT, however at DFT+U level the enthalpy values are overestimated. In the case of entropy ( figure 6) the curves of each system can now be distinguished, nevertheless the curves fit better with experimental information when  rutile takes a U=5 eV value and anatase U=4 eV than DFT standard calculation. In C v (and C p ) versus temperature plot (figures 7 and 8), from 400 K up once again, the three TiO 2 polymorphs behave similarly, approaching an asymptotic behavior at higher temperatures. Differences in these three TiO 2 polymorphs are not significant.
Thermal volumetric expansion can be determined and are plotted as shown in figure 9. TiO 2 polymorphs have the same behavior; hydrogen titanate does not experience any significant volumetric change when temperature is increased. Instead, a volume contraction is observed near 200 K which is correlated to the presence of several imaginary vibrational modes, that only appear for the cell volumes associated to that temperature, but that does not appear for the equilibrium values associated at 0 K. This fact could be connected to a possible phase transformation from hydrogen titanate to TiO 2 (B) in the same range of temperatures [44]. Anatase and rutile lead to bigger volumetric variations.
Thermodynamic properties fluctuate significantly when the Hubbard term is considered in our analysis. The previous differences noticed for phonon frequencies have an effect on the results of thermodynamic properties calculation.
In table 6 bulk modulus and volume variation using DFT and DFT+U at 0 K and 298 K is listed to compare the effect of addition of Hubbard parameter in our calculations [13,45,46].
The addition of Hubbard parameter leads to an increase in cell volume being less significant for TiO 2 -B and hydrogen titanate. The bulk modulus varies when U term is considered, being larger for anatase and TiO 2 -B, and smaller for rutile and hydrogen titanate. In the cases of TiO 2 -B and hydrogen titanate few data in literature is found, in order to compare our results, we have added two rows called 'DFT U=0' which correspond to our results at a DFT standard level.