Electronic Properties of Realistic Anatase TiO2 Nanoparticles from G0W0 Calculations on a Gaussian and Plane Waves Based Scheme

The electronic properties of realistic (TiO2)n nanoparticles (NPs) with cuboctahedral and bipyramidal morphologies are investigated within the many body perturbation theory (MBPT) G0W0 approximation using PBE and hybrid PBEx (12.5% Fock contribution) functionals as starting points. The use of a Gaussian and plane waves (GPW) scheme reduces the usual O4 computational time required in this type of calculations close to O3 and thus allows considering explicitly NPs with n up to 165. The analysis of the Kohn-Sham energy orbitals and quasiparticle (QP) energies shows that the optical energy gap (Ogap), the electronic energy gap (Egap) and the exciton binding energy (ΔEex) values decrease with increasing TiO2 NP size, in agreement with previous work. However, while bipyramidal NPs appear to reach the scalable regime already for n = 84, cuboctahedral NPs reach this regime only above n =151. Relevant correlations are found and reported that will allow one predicting these electronic properties at the G0W0 level in even much larger NPs where these calculations are unaffordable. The present work provides a feasible and practical way to approach the electronic properties of rather large TiO2 NPs and thus constitutes a further step in the study of realistic nanoparticles of semiconducting oxides.


INTRODUCTION
Accurate and affordable computational methods are required to properly describe the electronic properties of semiconducting materials. The prediction of these properties usually relies on band structure calculations carried out at different levels of theory. Since the 1970s, density functional theory (DFT) has become the standard tool to describe the physical and chemical properties of extended solids and finite systems. 1,2 Unfortunately, the broadly used local density approximation (LDA) and Generalized Gradient Approximation (GGA) implementations of DFT fail to predict the electronic properties of semiconductor materials resulting in a systematic underestimation of the calculated energy gap with respect to experiment. [3][4][5] To solve this critical issue several, more or less computational expensive, alternatives have been proposed including the use of the on-site Hubbard correction (the so called +U approach) [6][7][8] or the use of hybrid functionals including a part of non-local Fock exchange. [9][10][11][12] Both of these approaches exhibit a semiempirical flavor since U or the amount of Fock exchange need to be chosen, usually by comparing to experiment. Different approaches have been proposed, and efficiently implemented, in the past few years that are both physically well-grounded and accurate in terms of comparing to experiment. 13,14 In particular, we mention the recent implementations of the many-body perturbation theory based GW formalism proposed by Hedin more than 50 years ago. [15][16][17] GW based techniques constitutes actually the state-of-the-art due to a good performance in predicting the energy gap in molecules and nanoparticles and also the band gap of bulk solids. [18][19][20][21][22] However, the high accuracy of the GW methods is accompanied by a very unfavorable scaling with respect to the size of the system that makes it difficult to approach the size of realistic systems composed by hundreds/thousands of atoms. For the first-order perturbative approximation of the GW method, usually referred to as G0W0, the scaling with size is as high as O(N 4 ). [23][24][25] The O(N 4 ) dependence constitutes a real bottleneck when attempting to use these methods to approach the electronic properties of realistic systems such as semiconducting nanoparticles (NPs). These limitations have driven to reformulate the GW method for large scale calculations eliminating numerical approximations in the calculation of polarization, using efficient contour deformation techniques and providing a parallel implementation of the algorithm based on the splitting of both the single particle Green function and the screened Coulomb interaction. [26][27][28] Recently, an efficient algorithm to compute quasiparticle (QP) energies in the G0W0 formalism has displayed a reduction of the scaling of one order of magnitude to O(N 3 ), thus allowing to tackle more realistic systems. 29 This novel formulation based on a Gaussian and plane waves (GPW) scheme with Resolution of the Identity (RI) 30 is selected herein to investigate the electronic properties of realistic titanium dioxide (TiO2)n (n = 29-165) NPs composed of up to 495 atoms including cuboctahedral (or truncated bipyramidal) and octahedral (or bipyramidal) shapes as displayed in Figure 1 and described in more detail in a subsequent section. The results in the present work evidences that the scalable regime of these NPs is morphology dependent and unveils interesting correlations between Kohn-Sham orbital and G0W0 quasiparticle energies.

NANOPARTICLE MODELS AND COMPUTATIONAL STRATEGY
TiO2 NPs have been widely studied from theoretical and experimental points of view due to their applications related to various technologies, including catalysis [31][32][33] and photocatalysis. [34][35][36] In particular, photocatalytic activity is related to the electronic structure based properties of the TiO2 NPs 37,38 since these vary with size, shape and morphology as shown for models of realistic nanoparticles. 39 Clearly, high-level computational methods are required to provide accurate results and, hence, come out with reliable predictions.
A rather recent work investigated the performance of the G0W0 method for the prediction of the electronic energy gap of small to medium size (TiO2)n (n = 1-20) NPs, 40 as defined below, with structures obtained by a global optimization algorithm. 41 In spite of interesting trends, the TiO2 NPs considered in previous work, containing up to 60 atoms are too small to provide information about realistic systems and hence, analysis of larger TiO2 NPs is required to better link theoretical predictions with the experiment. Here, we focus on the electronic properties considering optical (Ogap), electronic (Egap) energy gaps, and exciton binding energy (ΔEex) of realistic titanium dioxide (TiO2)n NPs with n = 29-165. As in previous works, 39,40 Ogap is estimated from the difference between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) defined from the Kohn-Sham orbital energies. This is a reasonably good approximation as shown in a recent study comparing the thus defined Ogap to results arising from frequencydependent dielectric function; 42 and is in agreement with the interpretation proposed by Baerends et al. 43 On the other hand, Egap is calculated from differences in the QP energies (G0W0) for HOMO and LUMO (QPHOMO and QPLUMO) and thus corresponds to the difference between vertical ionization potential and electron affinity.
Finally, ΔEex is estimated as the difference between Egap and Ogap and provides information related to the electron-hole pair interaction energy. Note that, in absence of optical excitons, ΔEex vanishes when the system size becomes sufficiently large.
All calculations presented here have been carried out by using a GPW scheme for solving Kohn-Sham equations as implemented in the CP2K code. 44,45 The hybrid scheme employs a Gaussian basis to expand molecular orbitals and an auxiliary plane-wave basis for the expansion of the electronic density. In order to have an efficient expansion of the density in plane waves, core electrons are replaced by pseudopotentials.
Dual-space norm-conserving and relativistic pseudopotentials of the Goedecker-Teter-Hutter (GTH) type 46 parametrized for PBE 47 have been selected for oxygen (O) and titanium (Ti) atoms. 48 Importantly, we have adopted a recently implemented algorithm which for G0W0 calculations has a computational cost scaling close to O(N 3 ). 49 This makes it possible to tackle significantly larger systems than with the standard algorithm exhibiting a O(N 4 ) scaling.
In principle, GW calculations can be carried out in a self-consistent way and in that case final results do not depend on the initial density used to build the Green's function and screening potential. 50 However, those calculations are excessively demanding, even for simple systems containing a few atoms only. In practice, one relies on the G0W0 approach which leads to reasonable results 51 although these do depend on the initial starting density. 52 To investigate the effect of the initial density, the present G0W0 calculations are 4 carried out starting from the PBE density (G0W0@PBE) and also from that obtained from the PBEx hybrid functional including a 12.5% of non-local Fock exchange (G0W0@PBEx) which has proven to appropriately describe stoichiometric and reduced bulk anatase and rutile. 39 To carry out the G0W0@PBEx calculations, auxiliary density matrix method (ADMM) basis sets are used in order to deal with the hybrid part of the functional. 53 The ADMM basis sets have only been employed for the double-zeta primary basis sets, because they make the convergence of Hartree-Fock exchange (HFX) difficult with the standard HFX approach. For the primary Gaussian basis sets, we employ valence-only basis sets of double-zeta (DZ) and triplet-zeta (TZ) quality. The double-zeta basis sets are labeled as DZVP-MOLOPT-GTH for O atom and DZVP-MOLOPT-SR-GTH for Ti atom, 54 and triple-zeta basis sets are denoted cc-TZ for both oxygen and titanium. 55 In addition, the solution of the GW equations adopted here is made through the RI approximation. 56 For both primary basis sets, an auxiliary basis set denoted RI_TZ has been used for oxygen, and both the cc-TZ basis set and another basis set denoted RI have been used as auxiliary bases for titanium. The coefficients of the primary, RI, and ADMM basis sets, along with a sample CP2K input file, are presented in Supporting Information (SI).

RESULTS AND DISCUSSION
We start this section discussing the influence that primary basis set (see Tables S1 and S2 in SI) has on the Ogap and Egap of TiO2 NPs by using the mean absolute error (MAE) incurred by using the smaller DZ basis set with respect to the larger cc-TZ one. The MAE descriptor indicates 0.05 (0.07) and 0.19 (0.21) eV for Ogap and Egap, respectively, corresponding to the PBE and the hybrid PBEx exchange correlation functionals with the latter one in brackets. Note that there are no significant deviations in either energy gap (±0.02 eV) between PBE and PBEx. However, MAE is larger for Egap than for Ogap showing that the size of the primary basis set has a larger effect on QP energies regardless of the functional employed. Following previous convergence tests 57 and given the large dimensions of TiO2 NPs studied here, we assume that primary TZ basis set shows the best performance to describe properly the electronic properties of TiO2 NPs. Primary DZ basis overestimates Egap leading consequently, to overestimation of ΔEex compared with those values predicted with TZ basis (see Tables S1 and S2 for further details).
Once the effect of primary basis set is assessed, we tackle the analysis of the trend in the HOMO and LUMO orbital energies, and the subsequent Ogap with increasing NP size by using the Kohn-Sham (KS) orbital energies obtained by PBE and hybrid PBEx functionals. Later, an analogous discussion will be presented considering QP energies associated to HOMO and LUMO orbitals, and the subsequent Egap derived from G0W0@PBE and G0W0@PBEx based calculations. The former analysis depicted in Figure S1 shows two clear decreasing trends for both HOMO and LUMO KS orbital energies with TiO2 NP size regardless of the employed functional. The large values of the smallest sized TiO2 NPs correspond to the pronounced quantum confinement effect, which vanishes smoothly by increasing the TiO2 NP size. 58,59 An interesting result is the marked effect of the NP morphology. In general, TiO2 NPs with cuboctahedral morphology feature systematically larger HOMO and LUMO orbital energies (in absolute values) than NPs with bipyramidal morphology (see Figure S1). We attribute this systematic trend to the exposed surfaces; bypiramidal TiO2 NPs exhibit six (101) facets, whereas, cuboctahedral NPs expose simultaneously (101) and (001) facets. Similar tendencies are observed for values arising from calculations with the hybrid PBEx functional (see Figure S1, bottom panels). However, note that the PBEx functional promotes an opening of the energy gap by inducing downshifts and upshifts in HOMO and LUMO orbital energies, respectively. 60 In this way, the resulting Ogap shows a clear reduction of its value with increasing TiO2 NP size, in agreement with previous results. 61 Note, however, that the differences observed in the HOMO and LUMO orbital energies induced by the different morphology cancel out and, hence, a unique trend encompasses both families of TiO2 NPs. Only in the case of the (TiO2)97 NP the orbital energies are partially downshifted relative to this general trend, an effect due to its particular structure directly related with exposed (101) and (001) surfaces ratio. This is the only nanoparticle where the ratio between exposed surfaces is favorable for (001) and hence its Ogap has a lower value, as expected. 61 By using the first-order perturbation G0W0 approach, the ionization energies and electron affinities are calculated from QPHOMO and QPLUMO energies, respectively (see Figure 2). As mentioned earlier, these calculations demand the use of the larger auxiliary RI basis set raising significantly the computational cost of GW calculations. The effect of going from the auxiliary cc-TZ to the RI basis set (for Ti atoms) is analyzed in detail on the smallest (TiO2)n (n=29 and 35) NPs comparing the QP energies and Egap for cc-TZ as primary basis and with cc-TZ or RI as auxiliary basis sets by using PBE and hybrid PBEx functionals (see Tables S3 and S4). Note that each selected TiO2 NP is a representative case for each one of the TiO2 morphologies investigated here (see Figure 1). These calculations show that the change of auxiliary basis set from cc-TZ to RI promotes systematic downshifts and upshifts on QPHOMO and QPLUMO, respectively. Focusing on  (Table   S4); a similar behavior is observed for calculations performed with PBE functional (Table S3). The choice of the RI auxiliary basis set instead of the cc-TZ in the CP2K calculations is justified by the very good match between Egap values from the present G0W0 calculations for these two TiO2 NPs to those obtained by using a relativistic all-electron description with numerical atomic-centered orbital basis set as implemented in FHIaims code 62 (Tables S3 and S4 Figures S1 and 2 is not surprising because QP energies depend noticeably on the starting point, that is, on the "arbitrary" set of Kohn-Sham eigenvalues and orbitals 6 calculated by either PBE or PBEx functionals. The resulting Egap values follow the expected trend reducing its value with the number of TiO2 units (Figure 2). This result is fully consistent with a previous study where Egap is estimated by using ΔSCF approach based on all electron relativistic DFT including numerical atomcentered orbitals. 40 Next, we analyze the exciton binding energy (ΔEex) that provides the strength of electron-hole interaction. The trend on the calculated ΔEex provides a direct measure of the convergence of electronic properties toward the bulk where it vanishes due to the low exciton binding energy in bulk TiO2. This behavior is clearly shown in Figure 3 for both PBE and PBEx calculations and is in line with previous results based on ΔSCF calculations. 61 Looking at PBEx results, ΔEex exhibits a decreasing trend with the number of TiO2 units from 2.6 to 2.1 eV. The bipyramidal TiO2 NPs (red squares in Figure 3) have an exciton binding energy difference of 0.09 eV corresponding to the difference between (TiO2)84 and (TiO2)165. This result is an indication that, at this size, TiO2 NPs with bipyramidal shape are close to the scalable regime limit. This is, where properties of NPs scale linearly towards bulk limit values. On the other hand, for the cuboctahedral TiO2 NPs (blue squares in Figure 3) the ΔEex difference between (TiO2)78 and (TiO2)151 is as large as 0.16 eV.
This result indicates that TiO2 NPs with bipyramidal shape reach the scalable regime limit at lower sizes than truncated octahedral NPs. Note also that in the case of the cuboctahedral (TiO2)97 NP, ΔEex is far away from the general trend observed (see Figure 3), a clear consequence of the exposed (001) facets. In spite of that, for a given NP size, themagnitude of ΔEex distinguishes TiO2 NPs with different morphology, a new and important observation, following the tendencies towards the scalable regime limit.

CONCLUSIONS
In this work, a Gaussian and plane wave scheme has been used to approach the electronic properties of (TiO2)n (n = 29-165) NPs composed by up to 495 atoms at the G0W0 level of theory considering either PBE or hybrid PBEx functionals as starting points. These calculations are possible thanks to the reduction in the usual O 4 computational time to close to O 3 .
The cuboctahedral and bipyramidal NP morphologies exhibit different trends in the Kohn-Sham orbital and QP energies suggesting that this provides an electronic structure descriptor to discriminate them from photoemission experiments. In agreement with previous work, 61 Ogap and Egap decrease with the NP size regardless of the exchange correlation functional employed. The introduction of non-local Fock exchange in the functional induces a gap opening when going from PBE to PBEx, consequently Egap(G0W0) undergoes a similar gap opening giving larger Egap(G0W0@PBEx) value larger than Egap(G0W0@PBE). Also, a larger contribution of Fock exchange is observed in the calculation of Ogap, as expected from previous works. 63,64 The observed shift between PBE and PBEx functionals is also reflected on the caclulated ΔEex values.
From the more reliable PBEx results, it turns out that the bipyramidal TiO2 NPs reach the scalable regime at smaller sizes with ΔEex trending to bulk faster than for the NPs with cuboctahedral morphology.
Moreover, the correlation obtained between Egap and Ogap will be useful to predict the former quantity, at either G0W0@PBE or G0W0@PBEx level, for large TiO2 NPs just from the KS orbital energies through Ogap rather than from the QPs which would require large computational resources.
Finally, one must point out that the differences in the QPs predicted from G0W0@PBE and G0W0@PBEx indicate that calculated results are not fully converged within the G0W0 formalism. Nevertheless, it is expected that results arising from the use of a hybrid functional as starting point are closer to the experiment as inferred from such comparison in some systems. 65,66 8

ASSOCIATED CONTENT
The Supporting Information is available free of charge on the ACS Publications website at DOI: Tables S1 and S2 list the electronic properties of TiO2 NPs for DZ and TZ basis set with PBE and hybrid PBEx functionals, respectively. Figure S1 correspond to the evolution of HOMO and LUMO orbital KS energies with TiO2 NP size. Tables S3 and S4 compiles the comparison between CP2K and FHI-AIMS codes. Finally, the CP2K basis sets employed are described by using the CP2K inputs.

Notes
The authors declare no competing financial interest. Figure 1. Scheme of (TiO2)n (n = 29-165) nanoparticles, which are obtained by using top-down approach.
Anatase-type truncated octahedral (cuboctahedron) and anatase-type octahedral (bipyramid) shapes are distinguished by orange and green labels, respectively. Red and blue spheres stand for to O and Ti atoms, respectively. Note that bipyramid TiO2 NPs feature only the (101) surface and the cuboctahedron TiO2 ones exhibit both (101) and (001) surfaces. Further details about the atomic structure can be found in ref. 61.  . Exciton binding energy of (TiO2)n NPs with cuboctahedral (blue circles and squares) and bipyramidal (red circles and squares) shapes. Note that circles and squares correspond to PBE and PBEx calculations, respectively. Green and gray regions are included to distinguish the tendencies obtained with PBE and hybrid PBEx functionals, respectively.