KID Procedure Applied on the [(PY5Me2)MoO]+ Complex.

The KID (Koopmans in DFT) protocol usually applies in organic molecules of the closed-shell type. We used the KID procedure on an open-shell Mo-based system for the first time to choose the most suitable density functional to compute global and local reactivity descriptors obtained from the conceptual density-functional theory (DFT). From a set of 18 density functionals, spread from the second until the fourth rung of Jacob's ladder: BLYP, BP86, B97-D, MN12-L, MN15-L, M06-L, M11-L, CAM-B3LYP, PBE0, B3LYP, N12-SX, M06-2X, MN15, MN12-SX, ωB97X-D, M11, LC-ωHPBE, and APFD, we concluded that CAM-B3LYP provides the best outcome, and in the second place, M06-2X. Because the vertical first ionization potential and vertical first electron affinity in the ground state (gs) are defined as follows I = E gs(N - 1) - E gs(N) and A = E gs(N) - E gs(N + 1), where E gs(N - 1), E gs(N), and E gs(N + 1) correspond to energies of the system bearing N, N + 1, and N - 1 electrons, along with Koopmans' theorem (KT) given by I ≈ -εHOMO (εHOMO, highest occupied molecular orbital energy) and A ≈ -εLUMO (εLUMO, lowest unoccupied molecular orbital energy), the deviation from the KT was performed by the use of the index, such that J I = |E gs(N - 1) - E gs(N) + εHOMO| and J A = |E gs(N) - E gs(N + 1) + εLUMO|, which are absolute deviations from the perspective of I and A, respectively. Furthermore, the εSOMO (SOMO: singly-occupied molecular orbital energy) leads us to another index given by |ΔSL| = |εSOMO - εLUMO|. Therefore, J HL and |ΔSL| are indexes defined to evaluate the quality of the KT when employed within the context of quantum chemical calculations based on DFT and not the Hartree-Fock theory. We propose the index that could be more suitable to choose the most proper density functional because the J HL and |ΔSL| are independent indexes.

Introduction

Molybdenum, a second-row transition-metal, has become a promising candidate to synthesize new catalysts. Specifically, the scientific literature revealed findings involving this metal in H2 production from water scission.[1,2] Karunadasa et al.(3) found a Mo–N-based metal-organic compound whose molecular formula is given by [(PY5Me2)Mo(CF3SO3)]+, withPY5Me2 being the pentadentate ligand 2,6-bis[1,1-bis(2-pyridil)ethyl]-pyridine.[4]

In the beginning of this investigation performed by Karunadasa et al.,[3] the original global process was suggested by the following chemical reaction

However, these researchers experimentally confirmed[5] that [(PY5Me2)MoO]2+ undergoes two proton-coupled reductions, leading it to the [(PY5Me2)Mo(H2O)]2+ complex. This complex is reduced into [(PY5Me2)Mo(H2O)]+ followed by an intramolecular rearrangement turning it into the [(PY5Me2)Mo(H) (OH)]+ complex. The latter releases H2 and [(PY5Me2)MoO]+. This resulting complex has the feature of capturing electrons through multiple steps, leading to the aquo-Mo [(PY5Me2)Mo(H2O)]+ complex again, thus repeating the process of H2 release in neutral aqueous media and seawater. Therefore, the [(PY5Me2)MoO]2+ of eq is considered as a precursor that turns into the [(PY5Me2)Mo(H) (OH)]+ complex releasing H2 through an elementary step as depicted by eq

For that reason, we focused our attention on the [(PY5Me2)MoO]+ (Figure ) complex that takes part of this last and crucial elementary step involved in the catalytic cycle to release molecular hydrogen.

Figure 1

[(PY5Me2)MoO]+ complex.

One of the most recent investigations concerning compounds based on [(PY5Me2)MoO]+ complexes was aimed to analyze the electron-withdrawing and electron-donating nature of substituent groups along with their locations at the para-position of equatorial/axial pyridine rings in order to study its influence on the 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 whether it is located on equatorial or axial pyridine rings.[6] Furthermore, instead of using spin-density, the difference between electron-withdrawing and electron-donating substituent groups at axial para-positions can be distinguished visually utilizing the so-called dual descriptor.[7]

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.[5,8] Previous theoretical studies focused on the energy barrier, overall energy, substituent effect, and 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,6] However, all these theoretical approaches have assumed no dependence of the density functional, that is, up to now we assumed that the BP86[9,10] 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 and in agreement with electrochemical measurements. Nevertheless, reactivity descriptors usually do not link in the same way, so it becomes necessary to assess a representative set of density functionals to unveil whether BP86 is still suitable to evaluate reactivity or not in the framework of the conceptual density functional theory (CDFT).[11,12] Besides, this study is required because condensed values of local reactivity descriptors implemented in some quantum chemical or post-SCF software are based on the frontier molecular orbital approximation (FMOA), meaning that the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) drive the complete reactivity of a molecule. On the other hand, the electronic chemical potential and molecular hardness, the main two global reactivity descriptors coming from the CDFT, are usually computed according to the Koopmans’ theorem (KT).

In the present work, we used the KID (Koopmans in DFT) procedure proposed by Glossman-Mitnik and coworkers,[13−16] where a set of indexes called J quantify deviations when using the KT, that is, JI quantifies how much different I ≈ −εHOMO is when compared against I = Egs(N – 1) – Egs(N) and JA quantifies how much different A ≈ −εLUMO is when compared against A = Egs(N) – Egs(N + 1).

In order to determine a suitable density functional to be used in computing global reactivity descriptors because the number of available density functionals is too large so as to perform an analysis functional by functional,[17−19] we analyzed 18 density functionals spread between the second and fourth rung of Jacob’s ladder. It is worthwhile to mention that this is an alternative procedure to validate the application of the KT different than the one proposed by Bellafont et al.(20) While this methodology relies on the identification of the vertical first ionization potential with the negative of the HOMO energy and its validation through the comparison with known experimental values for a group of molecules, the KID procedure assess this validation by resorting to a comparison between the calculated values obtained through a ΔSCF and those given by the KT prescription, looking for internal coherence of the density functional and avoiding the need of the experimental results, which are not always available for large molecular systems. Moreover, this is based on the proven validity of this identification within the generalized Kohn–Sham theory,[21−23] which is valid not only for the vertical first ionization potential but in all cases. This means that the vertical first electron affinity must be equal to the negative of the energy of the last occupied orbital of the radical anion, which is usually called the singly-occupied molecular orbital energy (SOMO). Thus, if the orbital energies of the LUMO and the SOMO are similar for a given density functional, or in turn, the |ΔSL| = |εSOMO – εLUMO| descriptor is close to 0, then the practical extension of the KT to the calculation of the vertical first electron affinity based on the LUMO orbital energy is justified. This paves the route for the accurate determination of HOMO–LUMO gaps (H–L gaps), which are useful in numerous practical applications and whose validity is assessed through the JHL descriptor which has been formulated in such a way that good results by potential cancellation of errors can be avoided.

Theoretical Background

The review of the results obtained in the study was aimed at confirming the performance of the KID 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 the results obtained 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. The simplest conformity to the theorem of Koopmans is connecting εHOMO to −I and εLUMO to −A through three key indexes as follows: JI = |εHOMO + Egs(N – 1) – Egs(N)|, JA = |εLUMO + Egs(N) – Egs(N + 1)|, and , where at a given optimized geometry in the ground state, the term Egs(N) stands for the system energy with its original number of electrons (N), while Egs(N + 1) and Egs(N – 1) are the system energies with one more electron (N + 1) and one less electron (N – 1), respectively, of the same system with the same optimized geometry as it was obtained in the ground state.

It should be noticed that the JA 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 system with its original number of electrons, then the use of JA makes sense while LUMO tending to SOMO. For this reason, another descriptor |ΔSL| = |εSOMO – εLUMO| has been designed[13−16] to help in the verification of the accuracy of the approximation. Note that a |ΔSL| → 0 does not guarantee the fulfillment of the KT by itself, |ΔSL| is a control parameter for the JA, so discarding a possible JA → 0 which could be just a coincidence when |ΔSL| ≠ 0, in other words, JA → 0 does not make sense when SOMO is too different from LUMO. In molecules like mirabamides[24] |ΔSL| < 0.03, nevertheless, we cannot expect similar magnitudes when dealing with metal-organic or organometallic molecules as the analyzed Mo-oxo complex in the present work. On the contrary, we expect larger deviations, and hence, it is even more critical to minimize the error when applying the KT by selecting the most suitable density functional. Consequently, this article marks a precedent in the use of the KID protocol on open-shell transition-metal complexes because it will provide an idea about the order of magnitude of this deviation expressed in the J values.

To sum up, the lower the J values, the more suitable the KT application is. Because J = 0 is too hard a condition to meet, we propose a J threshold value equal to 0.5 to select the most proper density functional, while |ΔSL| is expected to be as small as possible. According to the latter, we can sum up all J indexes in one term, . Whence, the smaller JHLS is, the better density functional is.

The KT[25] should be applied only within the context of Hartree–Fock’s (HF) theory and not DFT. Nevertheless, Toro-Labbé and Zevallos demonstrated that the 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 HF theory.[26] The exact physical meaning could assign to the Kohn–Sham HOMO using the “Kohn–Sham analogue of KT in the HF theory”, meaning that the Kohn–Sham HOMO approximately equals the opposite of the vertical first ionization potential, −I.[27,28] Janak’s theorem[29] backs up this statement.

To end, it is essential to remark that it is not our interest to claim that this method can replace the comparison against experimental parameters or contrast with a high level of theoretical calculations as CCSD(T). The goal of the KID protocol is to unveil the density functional from a set of functionals that best abides the requirements of the KT because the local reactivity commonly is computed by first resorting to the FMOA, thus meaning that the HOMO and LUMO densities are employed to build up the local reactivity descriptors often used in the framework of CDFT.

Results

Table contains all results estimated from our quantum chemical calculations.

Table 1

HOMO (εHOMO), LUMO (εLUMO), and SOMO (εSOMO) Orbital Energies, H–L Gap, and the KID Descriptors (JI, JA, JHL, |ΔSL|, and JHLS all in eV) Tested in the Verification of the Koopmans-like Behavior of the Selected GGA, MGGA, HGGA, HMGGA, and LC-DFT Density Functionals for the Mo Complexa,b

DF	εHOMO	εLUMO	εSOMO	H-L gap	JI	JA	JHL	|ΔSL|	JHLS
GGA1
BLYP[9,30,31]	–5.6883	–5.4338	–2.3402	0.2544	1.681	1.567	2.298	3.094	3.854
BP86	–5.9764	–5.7160	–2.6020	0.2604	1.683	1.581	2.309	3.114	3.877
B97-Dc[32]	–5.7773	–5.5187	–2.3968	0.2544	1.681	1.600	2.321	3.122	3.890
MGGA2
MN12-Ld[33]	–5.7114	–5.2379	–2.2844	0.4735	1.593	1.499	2.187	2.954	3.676
MN15-Ld[34]	–5.7767	–5.3876	–2.3883	0.3891	1.636	1.570	2.268	2.999	3.760
M06-L[35]	–5.7223	–5.3604	–2.3081	0.3619	1.639	1.555	2.259	3.052	3.797
M11-Ld	–6.0526	–5.7661	–2.66664	0.2865	1.681	1.608	2.326	3.100	3.876
HGGA3
CAM-B3LYP[36]	–6.8780	–3.4014	–3.1005	3.4765	0.148	0.198	0.247	0.301	0.390
PBE0[37]	–6.2271	–4.5623	–2.6420	1.6648	1.043	0.976	1.428	1.927	2.399
B3LYP[30,38,39]	–6.1035	–4.7302	–2.5932	1.3734	1.181	1.070	1.594	2.137	2.666
N12-SX[40]	–5.8034	–4.8733	–2.2735	0.9301	1.397	1.305	1.912	2.599	3.226
HMGGA4
M06-2X[41]	–6.7261	–3.5927	–3.4975	3.1334	0.266	0.291	0.394	0.095	0.406
MN15[42]	–6.6140	–4.0357	–2.9636	2.5783	0.621	0.554	0.832	1.072	1.357
MN12-SXe[40]	–5.9438	–5.0224	–2.4727	0.9214	1.391	1.274	1.887	2.550	3.172
LC-DFT5
ωB97X-Dc[43]	–7.3297	–2.7815	–3.5620	4.5481	0.397	0.349	0.528	0.780	0.942
M11e[44]	–7.4763	–2.5905	–4.0619	4.8858	0.576	0.619	0.846	1.471	1.697
LC-ωHPBE[45−48]	–7.8494	–2.4667	–4.4371	5.3827	0.787	0.676	1.038	1.970	2.226
APFDc[49]	–6.2418	–4.6700	–2.6768	1.5717	1.099	0.993	1.481	1.993	2.483

The 6-311+G(d,p) basis set was employed.

1GGA: generalized gradient approximation. 2MGGA: meta generalized gradient approximation. 3HGGA: hybrid generalized gradient approximation. 4HMGGA: hybrid-meta generalized gradient approximation. 5LC-DFT: long-range corrected DFT.

D refers to addition of molecular mechanics.

Meta-non-separable gradient approximation.

Range-separated hybrid meta generalized gradient approximation.

Discussion

This section refers to all the results quoted in Table . In the first place, we have focused our attention on JHLS because it permits us to discard or accept the application of the KID protocol according to JHL and |ΔSL| values, thus revealing the most suitable density functional to be considered from the selected set.

The GGA (generalized gradient approximation), MGGA (meta generalized gradient approximation), HGGA (hybrid generalized gradient approximation), HMGGA (hybrid meta generalized gradient approximation), and LC-DFT (long-range corrected DFT) density functionals (DF) present the following average deviations from the KT: ⟨JHLS⟩GGA = 3.874 eV, ⟨JHLS⟩MGGA = 3.777 eV, ⟨JHLS⟩HGGA = 2.170 eV, ⟨JHLS⟩HMGGA = 1.645 eV, and ⟨JHLS⟩LC-DFT = 1.622 eV indicate that some types of HGGA, HMGGA, and LC-DFT density functionals can offer certain consistency with the KT but others still exhibit important deviations.

The best improvement is found within the HGGA set, where the CAM-B3LYP is the best DF from all considered in the present study. One explanation could be based on the construction of this functional; on the one hand, CAM-B3LYP is composed of 0.19 HF plus 0.81 Becke 1988 (B88) exchange interaction at a short-range. On other hand, it has 0.65 HF plus 0.35 B88 at the long-range, and, the standard error function describes the intermediate region. This result is related to the presence of an amount of HF exchange interaction. Evidence that favors this claim is detected through the difference between MN15 and MN15-L, clearly the inclusion of the HF exchange improves the performance of the J indexes. All of those density functionals lacking of HF exchange (0%) lead to worse JHLS values; BP86, B97-D, BLYP, M06-L, M11-L, MN12-L, and MN15-L. After including the HF exchange, JHLS values improve as: B3LYP (20%), ωB97X-D (22.2–100%), PBE0 (25%), N12-SX (25–0%), MN12-SX (25–0%), M11 (42.8–100%), MN15 (44%), and M06-2X (54%). Notice that we also included results coming from the use of the APFD density functional having a 23% HF exchange, it shows us better results than given by the DFs lacking the HF exchange but similar to that offered by density functionals of the HGGA-type.

Besides, it seems the improvement in the performance of the KT by a density functional that includes a certain percentage of HF exchange is accompanied by the split of the interelectronic Coulomb operator into a short-range part and a long-range part, the latter having a higher percentage of HF exchange. Consequently, we focused our attention on those density functionals where the separation of the interelectronic Coulomb operator is taken into account by differently weighting the HF exchange at short- and long-ranges. In the case of M11, we have 42.8% (short-range) and 100% (long-range); CAM-B3LYP, 19% (short-range) and 65% (long-range); LC-ωHPBE, 0% (short-range) and 100% (long-range); and ωB97X-D, 22.2% (short-range) and 100% (long-range).

Additional evidence supporting the use of the CAM-B3LYP density functional has to do with the H–L gap, such analysis results from the benchmark study performed by Tecmer et al.(50) Because 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, and CAM-B3LYP. The vertical excitation energies obtained were compared with the reference data obtained using accurate wave-function theory methods for the electronic spectrum of the UO22+. As a result, they found that CAM-B3LYP was able to produce the best agreement compared 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, JHLS = 0.390.

The M06-2X density functional gives us JHLS = 0.406 because of a noticeable resemblance between SOMO and LUMO energies leading to |ΔSL| = 0.095, which overcomes similar SOMO–LUMO given by the CAM-B3LYP functional. However, the closeness of the compliance of the KT JHL = 0.394 is not better than the one given by the CAM-B3LYP functional (JHL = 0.247).

JHLS indicates that GGA and MGGA density functionals should not be the first choice to apply the global reactivity descriptors from the CDFT under the KT approximation. Conversely, density functionals of the type HGGA, HMGGA, and LC-DFT conditioned to the inclusion of HF exchange partitioned into short- and long-ranges, where the latter range is favored with a higher amount of HF exchange or with a total of 54% HF exchange as given by M06-2X are suitable choices to be tested. The reader should keep in mind that using the KT allows one to get fast insights concerning global, and even local, reactivity. The JHLS index provides us a quantitative criterium to validate whether the KT is suitable to be used within the DFT context or not.

Conclusions

Currently, there are more than 200 density functionals available to execute quantum chemical calculations. Depending upon the users’ goals, certain functionals are more appropriate than others to run quantum chemical calculations. In the context of CDFT, 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 molecule reactivity. In this work, we have applied the KID procedure on an open-shell metal-organic cation complex for the first time to decide the best density functional from the KT point of view to use reactivity descriptors from the CDFT. 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 CDFT. An energetic difference equal to 9% between the SOMO and LUMO is acceptable to validate the applied criterion given by the KID protocol. Because this is not the only transition-metal complex to be analyzed, the next calculations will involve more Mo-based complexes participating in the catalytic cycle that supports the mechanism explaining the molecular hydrogen release. A similar analysis will be performed on the reactant and transition state from the chemical reaction shown in Figure . Another reason for using the KID protocol lies in the capabilities of most available Quantum Chemistry software. There are some Quantum Chemistry programs or others for performing a post-SCF analysis that use the FMOA to compute values of local reactivity descriptors and KT to get global reactivity values. When the user employs DFT as a level of theory to carry out the quantum chemical calculations, these programs do not warn the user whether the density functional in use is the most suitable to perform a reactivity analysis through FMOA and KT. This implies computing global reactivity descriptors’ values as chemical potential and molecular hardness[11,12] based on the HOMO and LUMO energies. These are Kohn–Sham and not HF molecular orbitals. Consequently, the KID protocol will help the user decide the best density functional to provide these requested HOMO and LUMO energies. Because JHL and |ΔSL| are independent J parameters, we defined that the JHLS encompasses both J indexes which makes it easier to choose the most suitable density functional.

Figure 2

Reactant (R), transition state (TS), and product (P) in the key step that releases molecular hydrogen. The [(PY5Me2)MoO]+ complex is produced and corresponds to the system under analysis in the present work.

Settings and Computational Methods

The [(PY5Me2)MoO]+ structure corresponds to a d3 complex with a net charge equal to +1 and a doublet spin-multiplicity. We used the unrestricted method based on the aforementioned density functionals provided by the Gaussian09[51] suite of programs, except the MN15 and MN15-L density functionals which were employed through the use of the Gaussian16[52] quantum chemistry package. The following exchange–correlation functionals[17−19] were used

GGA-type: BLYP,[9,30,31] B97-D,[32] and BP86.[9,10]

MGGA-type: MN15-L,[34] MN12-L,[33] M06-L,[35] and M11-L.[53]

HGGA-type: CAM-B3LYP,[36] PBE0,[37] B3LYP,[30,38,39] and N12-SX.[40]

HMGGA-type: M06-2X,[41] MN15,[42] and MN12-SX.[40]

LC-DFT-type: ωB97X-D,[43] M11,[44] and LC-ωHPBE.[45−48]

Exclusively of Gaussian: APFD[49] (Austin–Petersson–Frisch–Dobek density functional).

Besides we included the effective core potential, MWB28,[54−57] for the metal center and the standard Gaussian-type orbitals 6-311+G(d,p)[58−63] basis sets for N, C, O, and H atoms. For each density functional, 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). Then, we compared couples I and A against couples HOMO and LUMO energies through the use of the J indexes. We carried out harmonic vibrational frequency analyses to confirm that the obtained structures correspond to the minima on the Born–Oppenheimer surface; we did not find imaginary frequencies on the N-electron geometrical optimized structure for each used density functional.