Benchmark Density Functional Theory Approach for the Calculation of Bond Dissociation Energies of the M-O2 Bond: A Key Step in Water Splitting Reactions.

A very fascinating aspect in quantum chemical research is to determine the accurate and cost-effective methods for the calculation of electronic and structural properties through a benchmark study. The current study focuses on the performance evaluation of density functional theory methods for the accurate measurement of bond dissociation energies (BDEs) of chemically important M-O2 bonds in water splitting reactions. The BDE measurement has got noteworthy attention due to its importance in all areas of chemistry. For BDE measurements of M-O2 bonds in five metal complexes with oxygen molecules, 14 density functionals (DFs) are chosen from seven classes of DFs with two series of mixed basis sets. A combination of pseudopotential and Pople basis sets [LANL2DZ & 6-31G(d) and SDD & 6-31+G(d)] are used as a series of mixed basis sets. The B3LYP-GD3BJ functional with LANL2DZ & 6-31G(d) gives outstanding results due to low deviations, error, and the best Pearson's correlation (R) between the experimental and theoretical data. Our study suggested an efficient, low-cost, precise, and accurate B3LYP-GD3BJ/LANL2DZ & 6-31G(d) level of theory for BDE of the M-O2 bond, which may be useful for chemists working in the field of energy generation and utilization.

Introduction

Petroleum, coal, and natural gas are classical non-renewable sources of energy and cause environmental pollution by emitting CO2 to the atmosphere.[1] Therefore, the scientific community has been struggling over the past few decades to develop new and alternative renewable sources of energy and fuel.[2] The search for renewable sources of energies is further demonstrated by the facts that the power consumption requirement would be doubled by 2050, whereas the fossil fuels are depleting rapidly.[3] H2 is a green fuel, which produces water as byproduct after combustion. Water constitutes two third of the earth surface. It would be ideal if we can use water to produce hydrogen and then combust hydrogen to regenerate water.[4] In this perspective, catalytic water splitting using sunlight provides an attractive solution for a renewable energy source as well as a cleaner and greener future. Water splitting includes water oxidation and reduction. Water oxidation produces protons and electrons required to make renewable fuels.[5] Water is a plentiful and attractive candidate to be used as raw material. In this perspective, establishing a simple and superior catalytic system for efficient water oxidation is a challenging task. A number of systems including metal oxides to composite materials, noble metal complexes to transition-metal organometallics, mono to multinuclear site catalysts, and various water oxidation complexes have been investigated in a homogeneous environment and on the surfaces of photo or electrochemical conditions for water splitting.[6] This true catalytic system for efficient water splitting operates with four consecutive proton-coupled electron transfer steps to generate oxygen and hydrogen.[7]

Naturally some of the examples are present for production of clean fuel generation, for example, during photosynthesis, the tetra manganese oxygen evolving complex undergoes water oxidation, (see Figure for the ruthenium complex). This process involving four step consecutive proton coupled transfer cycle results in the generation of oxygen molecule, four protons, and four electrons. Synthetic chemists are using this idea for water splitting in artificial solar energy conversion complexes for fuel production.[7] Kurz screened a set of six multinuclear manganese complexes for catalyzing the oxygen evolution reactions under coherent experimental condition.[8] Zong and Thummel synthesized a series of three well-organized mono- and di-nuclear Ru complexes and added acetonitrile solution of Ru-catalysts to an aqueous (Ce(IV)-CF3SO3H) solution at 24 °C. Oxygen evolution is observed for both mono- and di-nuclear Ru-systems.[9] The experimental study is based on hit and trial, many times it did not give the desired products. Recently, theoretical studies are being used for studying water splitting reactions using various transition-metal complexes. Theoretically, water oxidation has been investigated by using different methods and softwares.[10] Baran and Hellman analyzed metal hangman-porphyrines as catalysts for the electrochemical reduction of O2 and oxidation of H2O.[11] The literature reveals that the rate-determining step in water oxidation involves the formation of oxygen molecule. So far, a well-established theoretical method, which accurately predicts a new catalytic system for water oxidation, is missing.

Figure 1

Catalytic water oxidation and oxygen molecule evolution mechanism by Ru complexes.

In this study, we aim to search an accurate method for the calculation of bond dissociation energies (BDEs) of the M–O2 bond, which is a key step in water splitting reactions. The literature reveals that the best way to explore an accurate method is a benchmark study (cost effective and quality evolving study).

Results and Discussion

Evaluation of DFs with LANL2DZ & 6-31G(d) Basis Sets

LANL2DZ & 6-31G(d) series of mixed basis sets are implemented with 14 density functionals (DFs), and their statistically analyzed results are given in Table and graphically represented in Figures –4.

Table 1

RMSD, SD, R, and MAE of M-O2 BDEs Calculated with Different DFs While Using LANL2DZ & 6-31G(d) Basis Setsa

classes of DFT	DFs	rmsd	SD	R	MAE
LDA	LSDA	9.99	8.17	0.57	–7.09
GGA	BP86	7.12	6.47	0.52	–0.37
meta-GGA	TPSSTPSS	7.21	7.41	0.52	–1.26
H-GGA	B3LYP	4.36	1.55	0.81	–1.29
	B3PW91	4.87	3.11	0.34	–1.08
	PBE0	6.67	6.22	0.45	–0.70
	MPWPW91	7.28	6.43	0.53	–0.10
	B97	12.20	9.15	0.54	3.41
GH meta-GGA	M05-2X	17.54	16.33	0.54	5.51
	M05	17.15	12.35	0.77	7.50
	M06-2X	17.50	13.82	0.56	6.37
	M06	5.74	5.48	0.38	–1.56
RS H-GGA	CAM-B3LYP	14.79	12.68	0.53	4.37
H-GGA-D	B3LYP-GD3BJ	4.12	1.23	0.88	–3.16

All values are given in kcal/mol, except R which is presented as fraction of 1.0.

Figure 2

RMSD of different DFs with two series of basis sets for M–O2 BDEs.

Figure 4

MAE of different DFs with two series of basis sets for M–O2 BDEs.

B3LYP-GD3BJ of H-GGA-D class with LANL2DZ & 6-31G(d) series of basis sets has shown the best performance for the BDE measurement of M–O2 bond. Mean absolute error (MAE) for this functional (B3LYP-GD3BJ) is minimized to −3.16 kcal/mol. The observed standard deviation (SD) and root-mean-square deviation (RMSD) are 4.12 and 1.23 kcal/mol, respectively (Figures –4). A high Pearson’s R (0.85) is also seen (Figure ). Although MAE is high, still other statistical parameters are low with good correlation. B3LYP functional is known to minimize the error in geometrical parameters and energetic calculations, as reported in the previous literature.[12] The modified form of B3LYP with dispersion correction having parameter for dispersion forces further minimized the geometric and thermodynamic errors of the B3LYP and give best results for structural and energetic property analyses.[13−18] Similar results are obtained in the current study where B3LYP-GD3BJ shows good results for BDEs of M–O2 bond. The functional gives best results with LANL2DZ and 6-31G(d) basis sets.

Figure 5

Pearson’s correlation (R) of B3LYP-GD3Bj with SDD & 6-31+G(d) basis set for BDE calculation of the M–O2 bond.

B3LYP functional of the generalized gradient approximation (GGA) class with LANL2DZ & 6-31G(d) got the second position in reproducing good results against the experimental BDE measurements of M–O2 bonds. A smaller error of −1.29 kcal/mol is observed for this functional (see Table ). Despite the SD is above 4 kcal/mol, the RMSD is 1.55 kcal/mol with lower R value of 0.69. B3PW91, B97, PBE0, and MPWPW91 functionals from H-GGA class have lower correlation with experimental data with higher deviation between 4.87 and 12.20 kcal/mol. As a result, this class is designated as a moderate performer for the desired data set. The results of statistical analyses indicate that the efficiency of these DFs further decreases with LANL2DZ & 6-31G(d) series except B3LYP, which gives good results for the reason mentioned above.

Among five functionals of GH meta-GGA (M05, M05-2X, M06, and M06-2X), M06 functional has shown good efficiency with LANL2DZ/6-31G (d) basis sets. The deviations and error are low for M06, but a very low R (0.38) is found for correction between the experimental and theoretical data. RMSD, SD, and MAE values are 5.74, 5.48, and −1.56 kcal/mol, respectively [smaller error is observed here with LANL2DZ/6-31G (d)]. Efficiency of the rest of the DFs (M05-2X, M05, and M06-2X) is further decreased due to more deviation, errors, and lower R compared to the experimental data. RMSD, SD, and MAE are in the range of 6.37–17.54 kcal/mol, and R is in the range of 0.54–0.77 (Figures –4). Correlation of M05 is 0.77, but deviations (RMSD and SD) and errors are high. Hence, good results of M06 reflect the average performer of the respective class for the desired data set. Zhao and Truhlar empirically fit the dispersion parameter in standard DFs so they can capture dispersion interactions. For example, they included a set of noncovalent interaction energies in the fitting of the Minnesota functionals (M06 and M06-L), which are recommended good performer for transition-metal complexes.[19] These functionals give good results for BDE of the M–O2 bond for water splitting.

On the other side, TPSSTPSS functional from meta-GGA class with LANL2DZ & 6-31G(d) series of basis sets shows exceptional behavior. MAE of TPSSTPSS is decreased up to −1.26 kcal/mol, but its RMSD and SD values are increased up to 7.41 kcal/mol (in comparison to the experimental data) along with a low R value (0.52) between the experimental and theoretical data. Hence, these statistical results illustrate the moderate performance of meta-GGA class, similar to GH meta-GGA class with LANL2DZ & 6-31G(d) series of basis sets. CAM-B3LYP functional of RS-HGGA is less efficient because the deviations (both RMSD and SD) are 14.79 and 12.68 kcal/mol, respectively. The R (0.53) is low, although the larger MAE of 4.37 kcal/mol is less compared to deviations (see Table ). Comparatively, there is drastic increase in the deviations of the CAM-B3LYP/LANL2DZ & 6-31G(d) method compared to already discussed good performer functionals of selected density functional theory (DFT) classes. LSDA functional of the local density approximation (LDA) class has the least efficiency for the BDE measurements of M–O2 bonds due to high deviations (up to 9 kcal/mol), errors (−7.09 kcal/mol), and a low R value of 0.57. The results reflect the worst performance of LDA class for the desired data set.

Overall, the B3LYP-GD3BJ functional from the H-GGA-D class has better performance among all selected DFs from selected DFT classes. The efficiency is more enhanced at LANL2DZ & 6-31G(d) basis sets in comparison to SDD & 6-31+G(d) basis sets (vide infra). Conclusively, B3LYP-GD3BJ functional with LANL2DZ & 6-31G(d) basis set is the best methodology for BDE measurements of M–O2 bonds for water splitting.

Evaluation of DFs with SDD & 6-31+G(d) Basis Sets

SDD & 6-31+G(d) series of mixed basis sets are selected with 14 DFs, and their results are statistically analyzed (see Figures –4).

B3LYP-GD3BJ functional from the H-GGA-D class sustains its good performance for the respective BDE study with SDD & 6-31+G(d) basis sets. RMSD, SD, R, and MAE of B3LYP-GD3BJ functional are 6.62, 2.47, 0.35, and −5.81 kcal/mol, respectively (Table ). Lower deviations and error justify the good efficiency of this functional but still these deviations and error are higher than their results at LANL2DZ & 6-31G(d) basis sets. The correlation of the B3LYP-GD3BJ is low compared to the experimental data. Statistical results suggest that the B3LYP-GD3BJ functional is a better choice for BDE measurements of the M–O2 bond, and the H-GGA-D class is observed as good one among all selected DFT classes. Lonsdale et al. analyzed that the inclusion of the dispersion parameter in the B3LYP method significantly describes loosely bonded electrons in the valence shell of the transition-metal complexes.[20]

Table 2

RMSD, SD, R, and MAE of M–O2 BDEs Calculated with Different DFs While Using SDD & 6-31+G(d) Basis Setsa

classes of DFT	DFs	rmsd	SD	R	MAE
LDA	LSDA	16.66	19.86	0.02	–12.45
GGA	BP86	14.36	27.08	0.48	2.37
meta-GGA	TPSSTPSS	7.85	18.97	0.03	–6.59
H-GGA	B3LYP	17.22	27.95	0.44	4.29
	B3PW91	14.92	24.03	0.42	3.75
	PBE0	6.97	18.78	0.02	–6.08
	MPWPW91	6.25	18.93	–0.02	–5.49
	B97	16.33	29.31	0.40	3.73
GH meta-GGA	M05-2X	5.06	28.22	0.04	–4.44
	M05	6.20	9.57	0.39	0.05
	M06-2X	18.65	16.65	0.82	6.49
	M06	3.98	23.99	0.13	–2.20
RS H-GGA	CAM-B3LYP	4.70	22.51	0.04	–4.08
H-GGA-D	B3LYP-GD3BJ	6.62	2.47	0.35	–5.81

All values are given in kcal/mol, except R which is presented as fraction of 1.0.

From the GH meta-GGA class, the best performance is observed for M05 functional along with SDD & 6-31+G(d) basis sets. Statistical results indicate that the deviations, RMSD and SD are 6.20 and 9.57 kcal/mol, respectively (Figures –4). MAE is 0.05 kcal/mol, and R between the experimental and theoretical data is 0.39. RMSD, SD, R, and MAE of M06 functional are 3.98, 23.99, 0.13, and −2.20 kcal/mol, respectively (Table ). These results illustrate the moderate efficiency of M06 with SDD & 6-31+G(d) basis sets for the M-O2 BDE measurement. M05-2X and M06-2X functionals have deviations between 5.06 and 28.22 kcal/mol, and R values are 0.04 and 0.82, respectively. Pearson’s correlation of M06-2X is high (0.82), whereas for M05-2X, it is low (0.04). Due to low RMSD and errors values, M05-2X is more efficient than the M06-2X functional. Overall, the M05 functional is an average performer for BDE measurements of the M–O2 bond.

CAM-B3LYP from the RSH-GGA class is less efficient for the BDE measurement of the M–O2 bond. RMSD, SD, R, and MAE of CAM-B3LYP/SDD & 6-31+G(d) are 4.70, 22.51, 0.04, and −4.08 kcal/mol, respectively (Table ). SD is high along with a low R value compared to CAM-B3LYP functional along with LANL2DZ & 6-31G(d) basis sets, but RMSD and error are decreased.

From the H-GGA class, B97, B3PW91, B3LYP, PBE0, and MPWPW91 functionals are selected for the current study. The deviations and errors of all these DFs are high, and on the other side, R values is low. Their RMSD, SD, and MAE values range from 3.73 to 29.31 kcal/mol. However, R value of these functionals ranges from 0.02 to 0.44, compared to the experimental data (Figures –6). The addition of exchange–correlation (XC) decreases the self-interaction error, but the static correlation error appeared due to XC inclusion. Due to this reason, the efficiency of this class is low for the treatment of transition-metal complexes.[21] These analyses reflect the less efficiency of the H-GGA class for required BDE measurement.

Figure 3

SD of different DFs with two series of basis sets for M–O2 BDEs.

Figure 6

Structures of transition-metal complexes having oxygen molecule with known experimental BDEs of M–O2 bond. In all the complexes, O2 show side on the binding model with selected transition-metal complexes.

TPSSTPSS from meta-GGA class has a high SD value of 18.97 kcal/mol, but RMSD and MAE are above 6 kcal/mol (Figures –4). Pearson’s correlation of 0.03 is also lower. On the basis of these results, TPSSTPSS functional is designated as a poor performer in the current study.

A lower efficiency is observed for each of the selected functional from GGA and LDA classes (BP86 from GGA and LSDA from LDA) when SDD & 6-31+G(d) basis sets are used. The deviations and errors are above 12 kcal/mol except the lower error observed for BP86/SDD & 6-31+G(d) level of theory (Table ). The low R value classifies these DFs as least efficient performers for the respective bond cleavage.

Among all selected DFs, B3LYP-GD3BJ of H-GGA-D class has the best performance for the BDEs of the M–O2 bond of the selected compounds. The reason is the inclusion of dispersion as largest correction, which tremendously influenced the magnitude of BDE. Previously, Hirao observed a similar effect during analysis of BDE of the methylcobalamin.[22] Jensen and co-workers also noticed that the dispersion corrected method are the best for determination of the BDE of the metal–phosphine bond, which is underestimated by other DFs. The dispersion corrected method generally gives much better agreement between the calculated and experimental data, as cleared from its higher correlation.[23]

Actually, weak hydrogen bonds and van der Waals interactions between different parts of organometallic complexes can be important for the reproduction of the observed structures and bond energies, but these forces are not completely described by standard DFT. These properties can be properly described by the dispersion-corrected DF. Dispersion can be described well by a simple potential function of the form C6R–6, where C6 is the dispersion coefficients and R is interatomic distances, which is implemented and parameterized by Grimme and co-workers.[17,24−26] In Grimme’s implementation, the use of such corrections is usually denoted by adding the extension “-D” to the standard abbreviations used for functionals, that is, B3LYP-GD3BJ adds these empirical dispersion corrections to the standard B3LYP functional as we used in the current study.[27]

In the current study, the highest efficiency of B3LYP-GD3BJ is observed with LANL2DZ & 6-31G(d) basis set for the BDE measurement of the M–O2 bond (see Table ). Due to all merits of the LANL2DZ & 6-31G(d) method in methodology, the B3LYP- GD3BJ performance with these basis sets is outstanding among all selected functionals.

Conclusions

For BDE measurements of M–O2 bonds in five metal complexes with oxygen molecule, 14 DFs from seven classes of DFT along with two series of mixed basis sets are selected. These series of basis sets include LAN2DZ & 6-31G(d) and SDD & 6-31+G(d) basis sets. The H-GGA-D class shows better performance. Among all selected DFs, the B3LYP-GD3BJ functional with the LAN2DZ and 6-31G(d) method shows outstanding results due to lower deviation, errors, and best R between the experimental and theoretical data. RMSD, SD, R, and MAE of the B3LYP-GD3BJ/LAN2DZ & 6-31G(d) method are 4.12, 1.23, 0.88, and −3.16 kcal/mol, respectively. This level of theory is considered as an excellent method for the BDE measurements of M–O2 bonds in metal complexes. LSDA functional of the LDA class is observed as the least efficient performer for desired data. The proficiency of selected DFT classes along with LANL2DZ & 6-31G(d) series of basis sets is described below

The H-GGA-D class of DFT sustains its better performance with other basis sets series [SDD & 6-31+G(d)]. The trend of selected basis sets for BDE measurements of M–O2 bonds is as follows

These theoretical benchmark studies not only justify the already reported experimental results but are also fruitful for experimentalists and theoreticians working on the reactivity of important M–O2 bonds and predicting new chemical pathways.

Computational Methodologies Adopted in the Current Work

Gaussian 09 software[28] is used for calculation of BDE. 14 DFs are chosen from seven classes of DFs with two series of mixed basis sets for BDE of M–O2 bonds. Selected DFs cover most of the DFT classes including local density approximation (LSDA),[29] GGA (BP86),[22] meta-GGA (TPSSTPSS),[30,31] hybrid GGA (B3LYP,[32] B3PW91,[33] B97,[34] MPW1PW91,[35] and PBE0),[36] global hybrid meta-GGA (M05,[19] M05-2X,[37] M06,[38] and M06-2X),[39] dispersion corrected hybrid GGA (B3LYP-GD3BJ),[40] and range separated hybrid GGA (CAM-B3LYP)[41] (see Table ). For the proper description of occupied orbitals, different classes of basis sets are selected. These basis sets include 6-31G(d) and 6-31+G(d) basis sets of Pople basis sets[42] and LANL2DZ and SDD of effective core potential (ECP) basis sets.[43,44] Hans Hellmann introduced pseudopotential or ECP approximation for the treatment of complex systems with simple description in 1934.[45] For light atoms of the periodic table (carbon, hydrogen, oxygen, etc.) Pople basis sets are used but, for heavy metal atoms, ECP basis sets are implemented. Pseudopotential explicitly treats only the valence electrons, whereas the core electrons are “frozen” and considered as rigid species as nuclei is taken. Reduction of core electrons decreases the computational cost by focusing on valence electrons, and as a result, basis set size is reduced. Moreover, ECPs can include relativistic effects, which are important in heavy elements. In selected ECP basis sets, Stuttgart–Dresden (SDD) basis set treats the inner core electrons with a constant pseudopotential and the valence electrons with triple zeta valence basis set. Therefore, its efficiency is more enhanced compared to LANL2DZ, which is double zeta valence basis set. LANL2DZ and SDD basis sets are used to describe transition metals and oxygen. The inclusion of the diffuse function in Pople basis sets exclusively explains the nature of loosely bonded electrons in outer shell orbitals of complexes and their resultant radicals. The literature reveals that Pople basis sets are observed best for BDE measurements of organic and inorganic molecules.[46]

Table 3

List of Selected DFs from Seven Classes of DFT

classes of DFT	DFs
local density approximation (LDA)	LSDA
generalized gradient approximation (GGA)	BP86
meta generalized gradient approximation (meta-GGA)	TPSSTPSS
hybrid generalized gradient approximation (H-GGA)	B3LYP
	B3PW91
	PBE0
	MPWPW91
	B97
global hybrid generalized gradient approximation (GH meta-GGA)	M05-2X
	M05
	M06-2X
	M06
range separated hybrid generalized gradient approximation (RS H-GGA)	CAM-B3LYP
dispersion corrected hybrid GGA (H-GGA-D)	B3LYP-GD3BJ

Five compounds are selected from the literature for the benchmark study of M–O2 BDE (for water splitting). The reason for the selection of these transition-metal complexes was their efficiencies for reversible binding with O2, which was experimentally proved. The experimental BDEs of M–O2 bond of transition-metal complexes with oxygen are already reported in the literature[47−49] (see Figure ), and their structural representations (1–5) are given in Figure .

For all selected molecules, optimization and frequency calculations are performed at the same level of theory to confirm all the structures as true minima. The zero-point corrected energy is taken for BDE of all selected bonds at 298 K and 1 atm, and the results are compared with the already reported experimental data.

For the M–O2 bond, all the complexes and their resultant radicals are studied up to four lowest spin states (see Figure for the dissociation pattern). For validation of theoretical methods with the experimental data, various statistical analysis tools including rmsd, SD, Pearson’s correlation (R), and MAE are used (the details of these parameters are given in the Supporting Information). These tools are well known in describing the best method for various data sets through benchmark studies.[50−55] Some studies recommend the use of MAE instead of rmsd because it possesses advantages of interpretability over rmsd. Moreover, MAE is the average absolute difference between two variables, and it is fundamentally easier to understand than the rmsd. In rmsd, square of difference is taken, which depends on prominent errors not on small errors, whereas the MAE counts the small errors. Beside this, each error contributing to MAE is proportional to the absolute value of the error, which is not the case for RMSD.[52] MAE measures the average magnitude of the errors in data set of predictions, without considering their direction. R depends both on the strength and the direction of relationship between two variables. In Pearson’s correlation, the relationship is linear when one variable is changed, and the other variable also changed with similar increasing or decreasing trend. In the current study, the method which has less deviations (RMSD and SD) and errors (MAE) compared to the experimental data along with a reasonable Pearson’s correlation (R) is considered as the method of interest.

Figure 7

Modal reaction for M–O2 bond dissociation of the rhenium-based complex with oxygen molecule. The rest of the complexes follow the same pattern.