Assessment of the KID Procedure on a Mo–oxo Complex, an Open–Shell System

The KID (Koopmans in DFT) procedure usually applies in organic molecules of the closed–shell type. We used the KID procedure in an open–shell system for the first time to choose the most suitable density functional to compute global and local reactivity descriptors coming from the Conceptual Density–Functional Theory. From a set of 18 density functionals spread from the second until the fourth rung of the Jacob’s ladder: BP86, B97-D, BLYP, CAM-B3LYP, M06-L, M11-L, MN12-L, B3LYP, PBE0, N12-SX, M06-2X, M11, MN12-SX, CAM-B3LYP, LC-ωHPBE, ωB97X-D, APFD, MN15 and MN15-L, we concluded that CAM-B3LYP provides the best outcome.

The global process is summarized by the following chemical reaction: These researchers experimentally confirmed that [(PY 5 Me 2 )Mo(H 2 O)] 2+ is able to turn water into H 2 through a suggested catalytic cycle [3] where the [(PY 5 Me 2 )Mo(H)(OH)] + complex can release H 2 , as a result the [(PY 5 Me 2 )MoO] + complex has the feature of capturing electrons through multiple steps, leading to the aquo-Mo complex which makes it an electrocatalyst to release H 2 in neutral aqueous media and seawater. Therefore the [(PY 5 Me 2 )MoO] + complex is the focus of our attention as Figure 1 depicts. The most recent investigation concerning compounds based on [(PY 5 Me 2 )MoO] + complexes was aimed to analyze the electron-withdrawing and electron-donating nature of substituent groups along with their locations at para-position of equatorial/axial pyridine rings in order to study its influence on kinetics (∆E ‡ ) and thermodynamics (∆E • ) of the molecular hydrogen release process. The results revealed an opposite effect given by each type of substituent group depending on if located on equatorial or axial pyridine rings [5]. Furthermore, a difference between electron-withdrawing or electron-donating substituent groups at axial para-positions can be distinguished visually utilizing the so-called dual descriptor better than spin-density [6].
An analysis of the global and local reactivity of Mo-oxo complexes as electrocatalysts is mandatory to attain a deep understanding of the molecular hydrogen release process [7,8]. Previous theoretical studies focused on: the energy barrier, overall energy, substituent effect, explicit and implicit solvent effect, along with the use of few local reactivity descriptors. All of them allowed to gain more insights concerning this type of electrocatalyst [3,5,7]. However, all these theoretical approaches have assumed none dependence of the density functional, i.e. up to now we assumed the BP86 density functional is accurate enough to carry out these analyses. Karunadasa et al. [3] suggested the use of BP86 owing its best performance to reproduce energetic values, which are the closest to those obtained experimentally in agreement with electrochemical measurements. Nevertheless, reactivity descriptors usually do not link in the same way, so being necessary to assess a representative set of density functionals to unveil whether BP86 is still proper to evaluate reactivity or not in the framework of the Conceptual-DFT (CDFT) [9,10]. Besides, this study is required because condensed values of local reactivity descriptors implemented in some quantum chemical or post-SCF softwares are based on the frontier molecular orbital approximation, meaning that HOMO and LUMO drive the complete reactivity of a molecule. On the other hand, electronic chemical potential and molecular hardness, the main two global reactivity descriptors coming from the CDFT, are usually computed according to the Koopmans' theorem.
In the present work, we used the KID procedure proposed by Glossman-Mitnik et al. [11][12][13][14] where a set of indexes called J quantify deviations when using the Koopmans' theorem. In order to quantify a suitable density functional to be used in computing global reactivity descriptors and since the number of available density functionals is too large as to perform an analysis functional by functional [15][16][17], we analyzed 18 density functionals spread in the first four rungs of Jacob's ladder. It is worthwhile to mention this is an alternative procedure to validate the application of the Koopmans' theorem different that one proposed by Bellafont et al. [18] where the biding energies assess the validity of the Koopmans' theorem.

Theoretical Background
The review of the results obtained in the study was aimed at confirming the fulfillment of the KID (Koopmans in DFT) protocol on an open-shell system. On doing it previously, several descriptors associated with the results that the HOMO and LUMO calculations obtained are related with results coming from the vertical first ionization potential (I) and vertical first electron affinity (A) following the ∆SCF procedure, where SCF refers to the Self-Consistent Field technique. There is a connection between the three key descriptors and the simplest conformity to the theorem of Koopmans by connecting H to −I, L to −A, and their actions by defining the HOMO-LUMO gap It should be noticed that the J A descriptor consists of an approximation which is only valid if the HOMO of the radical anion (the SOMO in the system with N + 1 electrons) resembles the LUMO of the neutral system, then the use of J A makes sense while LUMO tending to SOMO. For this reason, another descriptor |∆ SL| has been designed [11][12][13][14], to help in the verification of the accuracy of the approximation. Notice that a ∆ SL −→ 0 does not guarantee the fulfillment of the Koopmans's theorem, |∆ SL| is a control parameter for the J A so discarding a possible J A = 0 which could be just a coincidence when |∆ SL| = 0. In molecules like mirabamides [19] |∆ SL| < 0.03; however, we cannot expect similar magnitudes when dealing with metal-organic or organometallic molecules as the Mo-oxo complex analyzed in the present work; in consequence, this article marks a precedent in the use of the KID protocol on open-shell metal-transition complexes.
It should be noticed that if the HOMO of the radical anion (the SOMO in the system with N + 1 electrons) resembles the LUMO of the neutral system, then J A descriptor is a valid approximation, hence the use of J A makes sense while LUMO tending to SOMO. To sum up, the lower the J values, the more suitable the Koopman's theorem application is. Since J = 0 is a too hard condition to meet, then we propose a J threshold value equals 0.5 to select the most proper density-functional; while |∆ SL| is expected to be so small as possible.
The Koopmans' theorem [20] should be applied only within the context of the Hartree-Fock's theory and not DFT. Nevertheless, Toro-Labbé and Zevallos demonstrated that electronic chemical potential, µ = −0.5(I + A), and molecular hardness, η = I − A, based on Kohn-Sham's frontier molecular orbitals yield the same trends given by the same global reactivity descriptors based on frontier molecular orbitals of the Hartree-Fock theory [21]. The exact physical meaning could assign to the Kohn-Sham HOMO using the "Kohn-Sham analog of Koopmans' theorem in the Hartree-Fock theory", meaning that the Kohn-Sham HOMO approximately equals the opposite of the first vertical ionization potential, −I [22]. The Janak's theorem[23] backs up this statement.

Results
All results are quoted in Table 1.

Discussion
First of all, we must focus our attention on the |∆ SL| value because it permits us to discard or accept the application of the KID protocol according to J(HL) values for a certain density functional. Then we can define a threshold value for |∆ SL|, e.g., 1.0 so that all those density functionals resulting in |∆ SL| > 1.0 should be discarded because, in such cases, J A values make no sense since each one indicates that LUMO differs from SOMO as established by the threshold aforementioned. Even so, we included in discussion all obtained values since this is the first time that this protocol being applied on a transition-metal complex of the open-shell type. Notice that we also included results coming from the use of the APFD density functional having a 23% of HF exchange; it shows us better results than given by DFs lacking HF exchange, but similar that is offered by density-functionals of the HGGA-type.
However, the improvement in the performance of the Koopmans' theorem by a density functional must be constrained the split of the interelectronic Coulomb operator into a short-range part and a long-range part, having the latter a more percentage of HF exchange. Now the focus of our attention should be paid on the those density functionals including this separation of the interelectronic Coulomb operator and weighting the HF exchange at short-and long-range. In the case of MN11 we have 42.8% (short-range) and 11% X (long-range); CAM-B3LYP, 19% (short-range) and 65% (long-range); LC-ωHPBE, 0% (short-range) and 40% (long-range); and ωB97X-D, 22.2% (short-range) and 100% (long-range). The LC-ωHPBE presents the worst performance among the LC-DFT set of density functionals due to its lacking of short-range HF exchange (J(HL) = 2.034 with the 6-31+G(d,p) basis set and J(HL) = 1.038 with the 6-311+G(d,p) basis set). When one increases the basis set, the performance of LC-ωHPBE worsens since that despite being J(HL) = 1.038, the |∆ SL| = 1.970.
Additional evidence supporting the use of the CAM-B3LYP density functional that has to do with the HOMO-LUMO gap comes from the benchmark study performed by Tecmer et al. [24] Since Time-Dependent DFT (TD-DFT) results in excited states that sometimes appear lower in energy when compared against the reference state, they applied TD-DFT with the following set of density functionals: LDA, PBE, BLYP, B3LYP, PBE0, M06, M06-L, M06-2X, CAM-B3LYP. Vertical excitation energies obtained were compared with reference data obtained using accurate wave-function theory (WFT) methods for the electronic spectrum of the UO 2+ 2 . As a result, they found that CAM-B3LYP was able to produce the best agreement in comparison to coupled-cluster data, thus revealing the importance of proper inclusion of HF exchange as offered by a hybrid functional like CAM-B3LYP in transition-metal compounds. In fact, for the CAM-B3LP functional, |∆ SL| = 0.301, implies that SOMO differs from LUMO by about less than 9%, so that the variation from |∆ SL| = 0.271 up to |∆ SL| = 0.301 when changing the basis set from 6-31+G(d,p) to 6-311+G(d,p) is not significant.
Finally, as |∆ SL| indicates, GGA, MGGA, and HMGGA density functionals cannot be used to apply global reactivity descriptors coming from the Conceptual DFT when the Koopmans' theorem is applied. The reader should have in mind that the use of the Koopmans' theorem allows one to obtain values of global or even local reactivity descriptors because a fast analysis of reactivity through the use of global descriptors of the CDFT is carried out firstly through of the Koopmans' theorem. Then J indexes provide us quantitative criteria to validate whether the Koopmans' theorem is suitable to be used within this context or not.

Settings and Computational Methods
We used the Gaussian09 [25] suite of programs for 16 density functionals excepting MN15 and MN15-L that were employed through the use of the Gaussian16 [26] quantum chemistry package. The following exchange-correlation functionals [15][16][17]  Exclusively of Gaussian: APFD (Austin-Frisch-Petersson functional with dispersion) [27] Besides we included the effective core potential MWB28 [28][29][30][31] for metal center and the standard Gaussian-type orbitals 6-31+G(d,p) and 6-311+G(d,p) [32][33][34][35][35][36][37] basis sets for N, C, O, and H atoms. For each level of theory, we performed geometrical optimizations on the molecular system with its original number of electrons (N). Afterward, each optimized geometry was used in single-point calculations with N + 1 and N − 1 electrons in order to compute the vertical first ionization potential (I) and the vertical first electron affinity (A). Finally, we compared couples I and A against couples HOMO and LUMO energies through the use of the J indexes. We carried out harmonic vibrational frequencies analyses to confirm the obtained structures correspond to minima on the Born-Oppenheimer surface: we did not find imaginary frequencies on the N-electron geometrical optimized structure for each used density functional.

Conclusions
Nowadays, there are more than 200 density functionals available to execute quantum chemical calculations. Depending upon the goal established by users, certain functionals are more appropriate than others to execute quantum chemical calculations. In the context of Conceptual DFT, the choice of a density functional sometimes is an issue of personal preference based on specific criteria that could have nothing to do with the analysis of molecules reactivity. In this work, we have applied the KID (Koopmans' in DFT) procedure on an open-shell metal-organic cation complex for the first time to decide the best density functional from the Koopmans' theorem point of view to use reactivity descriptors from the Conceptual DFT. From the set of density functionals we selected, the best one corresponds to CAM-B3LYP, thus discarding the use of BP86 to perform any analysis concerning reactivity and selectivity through the use of Conceptual-Density-Functional Theory. An energetic difference as 9% between SOMO and LUMO is acceptable to validate the applied criterion given by the KID protocol. Since this is not the only transition metal complex to be analyzed, the next calculations will involve more based-Mo complexes participating in the catalytic cycle that supports the mechanism explaining the molecular hydrogen release. Alike, a similar analysis will be performed on the reactant and transition state from the chemical reaction given by Figure 2. Another reason for using the KID protocol is based on the fact that some programs of Quantum Chemistry use the Frontier Molecular Orbital Approximation to compute values of local reactivity descriptors. That means that HOMO and LUMO energies are not only input data to compute global reactivity through global reactivity descriptors as chemical potential and molecular hardness [9,10], but also Fukui functions [9,10] and dual descriptor [38][39][40] employ the squares of HOMO and LUMO as a first approximation to reveal the most susceptible sites to undergo nucleophilic and electrophilic sites of molecules. Funding: This research was funded by FONDECYT grant number 1181504, whose support is acknowledged by Jorge I. Martínez-Araya.