Energetic Study of Clusters and Reaction Barrier Heights from Efficient Semilocal Density Functionals

The accurate first-principles prediction of the energetic properties of molecules and clusters from efficient semilocal density functionals is of broad interest. Here we study the performance of a non-empirical Tao-Mo (TM) density functional on binding energies and excitation energies of titanium dioxide and water clusters, as well as reaction barrier heights. To make a comparison, a combination of the TM exchange part with the TPSS (Tao–Perdew–Staroverov–Scuseria) correlation functional—called TMTPSS—is also included in this study. Our calculations show that the best binding energies of titanium dioxide are predicted by PBE0 (Perdew–Burke–Ernzerhof hybrid functional), TM, and TMTPSS with nearly the same accuracy, while B3LYP (Beck’s three-parameter exchange part with Lee-Yang-Parr correlation), TPSS, and PBE (Perdew–Burke–Ernzerhof) yield larger mean absolute errors. For excitation energies of titanium and water clusters, PBE0 and B3LYP are the most accurate functionals, outperforming the performance of semilocal functionals due to the nonlocality problem suffered by the latter. Nevertheless, TMTPSS and TM functionals are still good accurate semilocal methods, improving upon the commonly-used TPSS and PBE functionals. We also find that the best reaction barrier heights are predicted by PBE0 and B3LYP, thanks to the nonlocality incorporated into these two hybrid functionals, but TMTPSS and TM are obviously more accurate than SCAN (Strongly Constrained and Appropriately Normed), TPSS, and PBE, suggesting the good performance of TM and TMTPSS for physically different systems and properties.


Introduction
Kohn-Sham density functional theory (DFT) [1,2] is the most popular electronic structure theory, and has achieved practical success in the first-principles prediction of the electronic structure of matter.In this theory, only the exchange-correlation energy component-which accounts for all many-body effects-remains unknown and has to be approximated as a functional of the electron density.As such, the predictive power of DFT largely depends on the density functional approximation to the exchange-correlation energy.Therefore, the development of exchange-correlation density functionals with better accuracy and wider applicability has been the central topic of the field [3][4][5][6][7][8].More importantly, achieving high accuracy for wide-ranging properties and diverse systems will greatly enhance the computer design of new materials and devices.
The exchange-correlation energy is defined by where n(r) is the electron density and ρ xc (r, r ) is the exchange-correlation hole at r around the electron at r (spin index has been suppressed and atomic units (h = e = m = 1) are used).Density functionals can be developed by imposing the exact constraints on the exchange-correlation energy E xc , or on the exchange-correlation hole ρ xc (r, r ), or on both.The constraints on the exchange-correlation energy are global, while the constraints on the exchange-correlation hole are local and thus are more difficult to impose.Most density functionals have been proposed from the energy constraints.While the functional form is rather restricted for the local spin-density approximation (LSDA) [3,4], the GGA (Generalized Gradient Approximation) [5] is more flexible and can incorporate more constraints.A meta-GGA [6][7][8] is even more flexible and can incorporate more local ingredients or inputs than LSDA and GGA.The simplest construction is the LSDA [3,4], in which the exchange part (i.e., Slater exchange) is calculated from the noninteracting plane wave, while the correlation part is calculated numerically with the quantum Monte Carlo (QMC) method.The LSDA is uniquely defined, although different parameterizations of the QMC data have been proposed.Since the typical valence electron density in solids is nearly uniform, the LSDA has been widely used in solid-state calculations.However, the LSDA tends to overestimate the short-range interaction arising from the orbital overlap of valence electrons in a molecule or solid and thus displays strong overbinding tendency.This tendency becomes more serious when the electron density becomes more inhomogeneous.For example, the LSDA yields moderately too-high cohesive energies of solids, but it significantly overestimates the atomization energies of molecules, in which the electron density is not slowly varying, leading to too-short bond lengths and lattice constants.GGA [5] is designed to correct the overbinding tendency of the LSDA.It was shown that GGA can significantly reduce the error of the LSDA in molecular atomization energy by a factor of about 10.However, this correction leads GGA to produce too-long bond lengths and lattice constants, and worsens the good performance of the LSDA for molecular solids, physisorption, and layered materials, in which the long-range van der Waals (vdW) interaction is important.A fundamental reason for this is that the excess short-range interaction of the LSDA effectively compensates for the absent long-range part of the vdW interaction.GGA removes the excess short-range interaction of the LSDA, but it is unable to provide the nonlocal long-range vdW interaction, because it is semilocal after all.This leads GGA to yield too-long bond lengths and lattice constants.In other words, GGA approaches experiments from the opposite side of the LSDA.To cure the shortcoming of GGA, meta-GGA functionals have been developed by including the orbital kinetic energy densities τ σ = ∑ occ i |φ iσ (r)| 2 as an additional input [6][7][8] beyond those of GGA.Using the kinetic energy density as another local ingredient makes the functional form of meta-GGA very flexible.This flexibility enables a meta-GGA to recover the correct slowly varying gradient expansion of the exchange-correlation energy up to fourth order in density gradient, and at the same time can reduce the self-interaction error of GGA.As a result, meta-GGA can extend the short-range part of the vdW interaction to a wider range.The former can improve the GGA description of usual bulk solids and molecules, while the latter can greatly improve GGA functionals for systems in which the long-range vdW interaction is important (e.g., molecular complexes and solids).There are two ways to achieve this goal.One is to include many parameters fitted to many sets of properties such as M06L (Minnesota 06 Local functional) meta-GGA [6], and the other is to incorporate many known constraints on E xc , such as TPSS (Tao-Perdew-Staroverov-Scuseria) [7] and SCAN meta-GGAs [8].
A more appealing approach for developing meta-GGA functionals is to incorporate exact constraints not only via E xc , but also via the hole ρ xc (r, r ).This will allow one to build more exact constraints into a meta-GGA than a meta-GGA constructed by only the constraints on the energy E xc , and thus has a better chance of achieving high accuracy.The starting point of this method is to approximate the hole directly, rather than approximating the hole with the reverse engineering approach [9], as is done in most density functionals.From the hole, one can calculate the underlying functional via Equation (1).Then, one can establish a desirable one-to-one relationship between a density functional and the associated hole.This is a natural method.The difficulty is that the conventionally defined exchange hole is nonlocal in nature.It was shown that this difficulty can be largely reduced with a general coordinate transformation [10], because the conventional exchange hole becomes localized under the general coordinate transformation.This offers a way to approximate the exchange hole with a semilocal functional.We have recently illustrated this approach by developing a nonempirical Tao-Mo (TM) meta-GGA functional [11].Numerical tests show that this functional can achieve remarkable accuracy for wide-ranging properties of molecules and solids [12,13].
Here we apply the TM meta-GGA to study the excitation and binding energies of titanium dioxide and water clusters, as well as reaction barrier heights.To have a better understanding of the performance of the TM functional, a combination of the TM exchange part with the original TPSS correlation (i.e., TMTPSS) is also included in this study.The selection of these systems is based on two considerations.First, these systems are formed with different bonding types (e.g., ionic bonds, covalent bonds, and hydrogen-bonds), and are thus physically different systems.Second, these clusters by themselves are of great interest, due to their remarkable properties and applications.All of these properties (i.e., excitation energies of clusters and reaction barrier heights) are closely related to electronic nonlocality.They represent a stringent test for semilocal DFT methods.Our calculations show that in all the cases, the TM and TMTPSS functionals consistently show better accuracy, improving upon commonly-used semilocal functionals.It is competitive with computationally more expensive hybrid functionals, pushing cost-effective first-principles prediction to higher accuracy.

Binding Energies and Excitation Energy of Titanium Dioxide Cluster
Titanium dioxide (TiO 2 ) has attracted great attention in the past few decades [14][15][16][17], due to its low cost, environmental compatibility, and experimentally proven potential for photocatalytic and photovoltaic applications [14,15].In particular, the use of TiO 2 as heterogeneous catalysts for photochemical splitting of water to produce renewable hydrogen and dye-sensitized solar cells has been the subject of intense research [14][15][16][17].
To understand the microscopic structure of TiO 2 , the binding energies of (TiO 2 ) n with n = 2-4 clusters in the gas phase are investigated.Figure 1 shows the chemical structure of (TiO 2 ) n , while Table 1 displays the binding energies of (TiO 2 ) n clusters from density functional calculations.From Table 1, we observe that the binding energies between TiO 2 molecules are in the range of a normal chemical bond.However, due to the large practicabilities of the O 2− ion and transition-metal ions, there exist strong vdW interactions in these ionic solid clusters.Because the long-range vdW interaction is missing in standard density functionals, this effect leads to large mean absolute errors (MAE) of all conventional DFT methods studied for this kind of cluster.Nevertheless, the TM and TMTPSS functionals give good improvement upon other semilocal DFT methods studied, such as PBE (Perdew-Burke-Ernzerhof) and TPSS.They are comparable to the more expensive hybrid PBE0 or even more accurate than B3LYP, as shown in Table 1.This performance is consistent with other recent studies [18,19].Table 1 also shows that the binding energy per molecule (TiO 2 ) increases with cluster size, as predicted by all DFT methods.This trend qualitatively agrees with MP2 prediction.For all methods including MP2, the increasing rate becomes smaller and should be gradually saturated to the cohesive energy of the bulk titanium dioxide when n > 7 [20].
Next, we study the lowest singlet excitation energies of titanium dioxide clusters with time-dependent DFT (TDDFT) within the adiabatic approximations.These excitation energies correspond to the HOMO-LUMO energy gaps of the cluster.In the large-size (n > 7) or bulk limit, they should converge to the band gap of bulk titanium dioxide [20,21].Figure 2 shows the lowest excitation energies of the (TiO 2 ) n (n = 1-4) clusters calculated from PBE, TPSS, TM, TMTPSS, B3LYP, and PBE0 within the adiabatic TDDFT.As shown in Figure 2, for small (TiO 2 ) n clusters with n = 1-4, all these excitation energies show strong odd-even oscillation in accordance with previous DFT results [20], and the predicted values are still far above the experimental results [20,21].All of these clusters show a large energy gap consistent with the previous calculations [20].However, the calculations overestimated the gaps considerably, probably because relaxation effects were not included [21,22].From Figure 2 we can see that all these DFT methods exhibit a similar trend for the lowest excitation energies of the (TiO 2 ) n (n = 2-4) clusters.All DFT lines show a drop of the excitation energy from 2 to 3 TiO 2 molecules, but MP2 shows an increase.There is a noticeable excitation energy shift between semilocal and hybrid functionals, leading to the excitation energies of hybrid functionals being closer to the accurate MP2 values.This shift clearly arises from the Hatree-Fock (HF) exchange built in these two hybrid functionals.The magnitude of the shift is dependent upon the amount of the HF exchange included in a hybrid functional.For example, PBE0 contains 1/4 HF exchange, while B3LYP contains 1/5 HF exchange.The HF exchange not only provides electronic nonlocality, which is necessary for some properties discussed above, but also improves the density tail behaviour of the exchange potential and reduces self-interaction error.As a result, hybrid functionals are more accurate than semilocal DFT for excitation energies, as argued by Becke and demonstrated in the literature [23].However, this does not mean that the more HF exchange contained in a hybrid functional, the more accuracy we can achieve, because one needs to consider the balance between the exchange and correlation parts or error cancellation.
Figure 2 also shows that the deviation of semilocal functionals from hybrid functionals becomes slightly larger as cluster size increases.Similar to hybrid functionals, the lowest singlet excitation energies of all the semilocal functionals will eventually saturate to the band gap of the bulk titanium dioxide.We observe that among the semilocal functionals, TM and TMTPSS are more accurate than others.In contrast, PBE gives the lowest excitation energies.While PBE0 does correct the exchange description relative to PBE, there is no correction whatsoever in PBE0 when it comes to dispersive interactions-a nonlocal correlation effect.

Binding Energies and Excitation Energy of Water Cluster
Water is the most abundant matter that covers about four-fifths of the Earth's surface.It has a complicated structure and anomalous properties that have not been well understood [24].Despite considerable efforts, many of the properties of water remain puzzling [25].Study shows that liquid water is not homogeneous at the nanoscopic level [26].To understand the detail of the structure, a fundamental understanding of water clusters is important [27,28].A feature of water clusters is their intermolecular interaction consisting of hydrogen bonds and a smaller part of the vdW interaction.Accurate description of these weak interactions presents a great challenge to DFT [29].Figure 1 shows the optimized structures of various water clusters, while Table 2 shows the comparison of the binding energies of water clusters (H 2 O) n (n = 2-5) between DFT methods and MP2.From Table 2, we can observe that TM can achieve high accuracy, substantially improving upon other DFT methods.
From Table 2, one can find that PBE0 performs worst.This is because PBE0 corrects the exchange description relative to PBE, but there is no correction whatsoever in PBE0 when it comes to dispersive interactions-a nonlocal correlation effect.Dispersion plays a vital role in the energetics of water systems.Table 2 also shows the results of PBE-D3BJ and TPSS-D3BJ developed by Grimme [30,31] for comparison.It is shown that the addition of dispersion (DFT-D3) to a semi-local functional always increases binding energy , as expected [29,32], but is not correctly described by the semi-local functionals.The same results are obtained by others [29,32].The addition of dispersion to a semi-local functional can bring large changes in the structure and equilibrium density of the liquid, so that errors are large [29].
To have a better understanding of water clusters, we calculate the lowest vertical excitation energies of a few small water clusters, in which the accurate ab initio values are available for comparison.The results are displayed in Figure 3. From Figure 3, we observe that the lowest vertical excitation energy of a water cluster quickly saturates to the band gaps of ice, which are smaller than experiment (8.8 ± 0.4 eV) [33][34][35][36][37][38][39][40], but with PBE0 being the closest one, followed by B3LYP and semilocal DFT methods for the same reason discussed above.We also observe that TM and TMTPSS give good match results with MP2 value among semilocal DFT methods, although the discrepancy is still too large due to the nonlocality issue.This indicates the robustness of the TM exchange.From Figure 3, we can also observe that all the calculated excitation energies tend to be red-shifted with respect to the MP2 values, as expected.

BH76 Barrier Heights
Prediction of reaction barrier heights presents a great challenge to semilocal DFT, due to the presence of stretched bonds in a transition state, which is closely related to electronic nonlocality.
The BH76 (76 Barrier Heights) test set consists of 38 hydrogen transfer barrier heights and 38 non-hydrogen transfer barrier heights.It is a standard test set that has been widely used to benchmark the theoretical methods [41,42].Table 3 shows the statistical errors of the reaction barrier heights evaluated with the TM and TMTPSS functionals, compared with other density functionals.
Table 3. Summary of deviations of the calculated reaction barrier heights from best values [41] for the BH76 test set.Results are taken from Reference [42] for M06L, Reference [8] for SCAN, and Reference [41] for local spin-density approximation (LSDA), PBE, TPSS, B3LYP, and PBE0.All values are in kcal/mol.ME = theory − best values from Reference [41].From Table 3, we can see that like other semilocal functionals, the TM and TMTPSS functionals tend to underestimate reaction barrier heights.The order of MAE for the BH76 set with different methods is PBE0 < M06L < B3LYP < TMTPSS < TM < SCAN < TPSS < PBE < LSDA.Reaction barrier heights are generally predicted more accurately with nonlocal density functionals, because such functionals exhibit a small delocalization error [43] for species with stretched bonds (i.e., transition states).From Table 3, we can also see a large reduction of errors in reaction barrier heights from non-hybrid to hybrid functionals (e.g., from PBE to PBE0).Since BH76 was used to parameterize the M06L meta-GGA functional, we have not included it for comparison, but list it in Table 3 for reference.

Materials and Methods
All calculations of clusters were carried out using the modified Gaussian 09 [44] with 6-311++G(3df,3pd) basis set.The ground-state geometries of molecules and clusters were performed self-consistently with respective density functionals.Based on the optimized geometries, we evaluated the lowest vertical excitation energies of clusters with TDDFT within the adiabatic approximation.The binding energies of (TiO 2 ) n clusters were calculated as energy difference between cluster and constituent molecules.Since our calculations involve the treatment of both the ground state (optimization of molecular geometry) and excited states, for consistency, the same basis was used in all calculations.The ultrafine grid (Grid = Ultrafine) in numerical integration and the tight self-consistent field convergence criterion (SCF = Tight) were used.To make comparison consistent, we performed the calculations with all density functionals, rather than taking the data directly from the literature for (H 2 O) n and (TiO 2 ) n clusters.Since reaction barrier height calculations are highly sensitive to the basis set, in this work, TM and TMTPSS methods were used to perform the geometries optimization and frequency analysis of the reactants, transition states, and products with MG3 (Multi-coefficient Gaussian-3) basis set by using their QCISD(Quadratic configuration interaction single and double excitation methods)/MG3 structure [41].Throughout the paper, ME and MAE are mean error (deviation) and absolute mean error (deviation) respectively.

Conclusions
We have studied the binding energies and excitation energies of clusters, as well as the BH76 barrier heights, which cover all bonding types, with recently proposed semilocal functionals.In this study, commonly-used semilocal and hybrid density functionals have also been included for comparison.Our calculations show that the TM and TMTPSS functionals can achieve good performance for these properties and systems, improving upon commonly-used semilocal DFT methods.We find that in some cases, they are competitive to the most popular hybrid functionals, but with lower computational cost.
A striking feature of the TM functional is that it incorporates many exact constraints through the underlying exchange hole, such as uniform coordinate scaling [45], negativity, correct uniform-gas limit, and spin scaling relationship [46].The high accuracy of the TM functional may greatly benefit from its exchange enhancement factor, which has the same small oscillatory behavior as M06L [47].Therefore, the TM functional can capture or extend the short-range part of the van der Waals interaction due to the de-enhancement (relative to the LSDA) in some regions, leading to an improvement for most properties.The TM correlation functional was developed separately from the exchange part, which respects all the exact conditions that the TPSS correlation satisfies, and improvement over TPSS in the low-density (strong-interaction) limit.Therefore, it is not surprising that the TMTPSS (TMx+TPSSc) functional can give the best results of energy calculation in many systems studied.