Metastability Triggered Reactivity in Clusters at Realistic Conditions: A Case Study of N-doped (TiO$_2$)$_n$ for Photocatalysis

Here we report a strategy, by taking a prototypical model system for photocatalysis (viz. N-doped (TiO$_2$)$_n$ clusters), to accurately determine low energy metastable structures that can play a major role with enhanced catalytic reactivity. Computational design of specific metastable photocatalyst with enhanced activity is never been easy due to plenty of isomers on potential energy surface. This requires fixing various parameters viz. (i) favorable formation energy, (ii) low fundamental gap, (iii) low excitation energy and (iv) high vertical electron affinity (VEA) and low vertical ionization potential (VIP). We validate here by integrating several first principles based methodologies that consideration of the global minimum structure alone can severely underestimate the activity. As a first step, we have used a suite of genetic algorithms [viz. searching clusters with conventional minimum total energy ((GA)$_\textrm{E}$); searching clusters with specific property i.e. high VEA ((GA)$_\textrm{P}^{\textrm{EA}}$), and low VIP ((GA)$_\textrm{P}^{\textrm{IP}}$)] to model the N-doped (TiO$_2$)$_n$ clusters. Following this, we have identified its free energy using ab initio thermodynamics to confirm that the metastable structures are not too far from the global minima. By analyzing a large dataset, we find that N-substitution ((N)$_\textrm{O}$) prefers to reside at highly coordinated oxygen site to maximize its coordination, whereas N-interstitial ((NO)$_\textrm{O}$) and split-interstitial ((N$_2)_\textrm{O}$) favor the dangling oxygen site. Interestingly, we notice that each types of defect (viz. substitution, interstitials) reduce the fundamental gap and excitation energy substantially. However, (NO)$_\textrm{O}$ and (N$_2)_\textrm{O}$ doped clusters are the potential candidates for overall water splitting, whereas N$_\textrm{O}$ is congenial only for oxygen evolution reaction.


Introduction
Accurate prediction of the structure of clusters as a catalyst, under reaction conditions, is the most fundamental challenge to get a detailed understanding of the active sites and their importance. Determining the catalyst structures at various reaction conditions is still a great challenge even for modern experimental methods. First principles based state-ofthe-art global optimization methods viz. genetic algorithm (GA), 1-6 basin-hopping (BH), 7 parallel tempering, 8 particle-swarm optimization (PSO), 9 stochastic tunneling, 10 simulated annealing (SA) 11 etc. can predict the catalysts' structure. Moreover, if we know some primary information such as the elemental constituents in the catalyst and its reaction conditions (e.g. temperatures and pressures, doping concentration, etc.), accurate prediction of the equilibrium state is in principle possible via ab initio thermodynamics. 12 However, this situation becomes complicated for real catalysis if we look deep into practical correlation of the predicted structure and its relevance with the concerned catalytic reactivity. The structures closer to the global minimum (based on ground state total energies) have higher occurrence probability at a finite temperature, but that does not ensure these structures are responsible for the observed activity. 13 On the other hand, metastable isomers of the catalyst are definitely not as stable as its global minimum but may lead to having higher activity due to presence of active sites. Moreover, under reaction conditions, the catalyst comprises of a wide range of structures all of which could be active to some extent in the catalytic reaction. 14 Over the past, in most of the theoretical studies of clusters, it is assumed that experimentally probed clusters are in their ground state conditions due to thermodynamic equilibrium. [15][16][17][18][19] Contrary to this, few experimental and theoretical studies have revealed that the experimentally detected clusters are the metastable isomers, rather than the ground-state ones. 20- 23 Here we present a robust theoretical approach to study the active sites of a cluster by taking a prototypical model system for Transition-metal (TM) oxides, in particular titanium dioxide (TiO 2 ), owing to its ubiquity, low cost, stability, nontoxicity, catalytic activity and environment friendly nature. TiO 2 has a great significance in photocatalysis from the perspective of industrial applications. [24][25][26][27][28][29][30][31][32][33] However, the wide band gap of TiO 2 that only absorbs the UV light of the solar spectrum, limits its efficiency in technological applications. Previous works suggest that the non-metal doping enhances the photoactivity of nanoclusters of TiO 2 . [34][35][36][37][38][39] The higher photocatalytic activity and stronger optical response are noticed with an increase of the N-doping concentration in TiO 2 nanophotocatalyst. 35,38,[40][41][42][43][44][45][46][47][48] Xiaobo et al. have revealed additional electronic states for non-metal dopants (N, C and S), above the valence band edge of pure TiO 2 nanomaterials using X-ray photoelectron spectroscopy (XPS), 38 which lead to the substantial modification in the optical properties. On the contrary, in theoretical study of Shevlin et al., no response in visible region is reported for N-doped TiO 2 clusters. 49 The viable reason of this discrepency of theoretical and experimental finding is, in the theoretical work they have addressed only the most stable substitutional N-defects, whereas in the XPS spectra both substitutional and interstitial (meta)stable defects are detected. 38 Therefore, it's well known these days that electronic properties of nanoclusters highly depend on structural configuration and particularly for the catalysis purpose, metastable structures are promising choice rather than the global minimum structure. 13,18,20,22 Note that previous studies have suggested that clusters possessing a high vertical electron affinity (VEA) or a low vertical ionization potential (VIP) are the promising choice as a photocatalyst. This is due to their ability to accept or donate an electron more readily. 20,22, 23 We have, therefore, implemented a suite of massively parallel cascade genetic algorithms (GA).
The first is the conventional energy-based GA [viz. (GA) E ] as described in detail in Ref. 12 This will give us all the local isomers close to energy based global minimum. The second GA is tailored explicitly to find metastable structures having specific bias for a property [viz. (GA) P ]. This specific property is used to evaluate the fitness function for (GA) P . If this property is high vertical electron affinity, we call it (GA) EA P , whereas if it's low vertical ionization potential we represent it (GA) IP P . As a test case we have shown the performance of these three GAs [viz. (GA) E , (GA) EA P , (GA) IP P ] for pristine (TiO 2 ) n clusters at various sizes in Fig S1 of supporting information (SI). It's clearly shown that while (GA) E searches low energy clusters, (GA) P focusses some metastable part of the PES to optimize some specific properties. More details and validity of this implementation can be found in Ref. 23 Here, we have applied these three GAs to build a database of pristine as well as doped (TiO 2 ) n clusters with n = 4 -10, 15, 20. Note that we have investigated three different configurations In this article, as a first step from an exhaustive scanning, we have considered three types of (un)doped (TiO 2 ) n clusters: (i) clusters having the minimum ground state total energy (ii) clusters with high vertical electron affinity (VEA), and (iii) clusters possessing the low vertical ionization potential (VIP). Note that despite (meta)stable structures are promising candidates for catalysis, their free energy of formation should not be too far away from the free energy based global minimum. Therefore, we determine the thermodynamic stability of these structures by minimizing its Gibbs' free energy of formation as a function of charge state at realistic conditions (e.g. temperature (T ), oxygen partial pressure (p O 2 ), doping). 12, [50][51][52][53] This facilitates us to estimate the probability of occurrence of these (meta)stable structures. Following this, a few clusters, that are thermodynamically stable as well as possess active sites, are selected and their electronic structures are accurately analyzed using GW calculations. This is how we have systematically studied doped (TiO 2 ) n clusters for application in photocatalysis.

Methodology
All density functional theory (DFT) calculations have been performed using FHI-aims code, which is an all electron code with numeric, atom-centered basis set. 54    We have found that O-site with high coordination is preferable for substitutional defect.

Results and Discussions
This finding is in good agreement with the previous simulation study. 49 Note that from the perspective of electronic configuration, N-atom has three unpaired electrons in the outermost shell, whereas O-atom has only two electrons. Therefore, the N-atom will favor more folded O-site to attain the maximum coordination number in the clusters. For interstitial case, we have observed that interstitial N-atom prefers to form the bond with dangling O-atoms because these oxygen atoms have the localized states at the HOMO level (see Fig S3).
The interaction of (NO) O doping affects the HOMO states, which may lead to substantial modification in their electronic properties. Moreover, since N-atom has unpaired electrons, it has tendency to share the electrons with more electronegative atom (dangling O-atom is having less coordination number). Likewise (NO) O , (N 2 ) O also favors the dangling O-site (see (i) The formation energy of (N) O in the charge state q is given by: (ii) The formation energy for interstitial N i.e. (NO) O is: (iii) Similarly, the formation energy for (N 2 ) O can be written as: where m is the mass, I A is the moment of inertia of O 2 molecule, M is the spin multiplicity and σ is the symmetry number. Similarly, µ N is estimated by the formation of N 2 molecule, i.e., ∆µ N = -0.25 eV at ambient condition. 19,64 3D phase diagrams of (GA) E clusters for size n = 5, 10, 15, 20 are shown in Fig 2(a), O defect is also observed at higher range of pressure for size n = 10, 15, 20. At higher values of µ e (n-type doping), (N 2 ) 0 O defect is stable at lower range of pressure except for n = 10 case. However at feasible pressure, O defect is found to be the most stable. Hence, we can summarize that interstitial defects [(NO) O and (N 2 ) O ] are most stable in (GA) E based (TiO 2 ) n clusters. However, substitutional defect N O has higher formation energy at ambient condition (T = 300 K, p O 2 = 1 atm), that can be seen clearly in 2D phase diagrams as shown in SI (Fig S4(a-d)). Note that until date, theoretical calculations are limited to non-metal doping (that to substitution only) in TiO 2 clusters closer to energy based global minimum. 37,49 It is therefore of profound interest to address the stability and electronic structure of property based doped clusters (i.e.  We have then determined the fundamental gap (E g ) and excitation energy of all the N-doped (TiO 2 ) n clusters and compared with pristine counterpart (see Fig 3a). These are computed at the level of G 0 W 0 @PBE0. Note that the difference of vertical electron affinity (VEA) and vertical ionization potential (VIP) gives the fundamental gap. 65 We can also de-  Table S1). This means these are metastable isomers and any conventional total energy based global minumum search algorithm will miss them to detect.
In Fig 3a, we have shown the fundamental gap vs. excitation energy of (un)doped clus-   Using the information of VIP and VEA, one can calculate the reduction potentials associated with the free charge carriers. Specifically, we use many body perturbation theory to calculate the thermodynamic driving force for the water splitting half-reactions 5 and 6. Experimental potential values are given relative to the Standard Hydrogen Electrode (SHE) (pH = 0). In practice, the required overall potential difference is larger than 1.23 eV to overcome energetic losses and kinetic barriers. In Fig 4(b, c, d), we show the G 0 W 0 @PBE0 predicted vertical potentials relative to the SHE to obtain the potential candidates for photocatalytic water splitting. The colored area (as in Fig 4(b, c, d)  In (N) O dopant cases, the spectrum has shifed to the lower frequency in comparison to the undoped cluster (see Fig 6(a-d)). Contrary to this, in the IR spectra of (NO) O and (N 2 ) O doped clusters, we have noticed the additional peaks above the highest peak of undoped cluster as shown in Fig 6(e-h).  I. VEA, VIP and Relative energy of pristine (TiO 2 ) n clusters In Fig S1(a, b,  has scanned the clusters with minimum energy (see Fig S1a), O is the most stable for n = 10 and 15, whereas (N 2 ) 0 O is stable for n = 20 (see the Fig S4(b-d)). In (GA) P EA case, near p-type doping region (N 2 ) +1 O is the most sable phase for n = 5, 10, and near n-type doping region (NO) 0 O becomes favorable as shown in Fig S4( V. Excitation energy of undoped (TiO 2 ) n clusters