The Use of Density Functional Theory to Decipher the Electrochemical Activity of Metal Clathrochelates with Regard to the Hydrogen Evolution Reaction in the Homogeneous Phase

The energetic needs of a rising human population have led to the search for alternative energy sources. A promising route for the large-scale storage of renewable energy is water electrolysis, which is performed with a proton-conducting polymer electrolyte. However, only platinum group metal electrocatalysts have the adequate properties to minimize the overvoltages associated with either hydrogen or oxygen evolution reactions. Alternative materials based on transition metals are scarce, but molecular electrochemistry offers some alternatives. In particular, transition metal clathrochelates exhibit an interesting activity with regard to the hydrogen evolution reaction (HER). However, such complexes form a vast family, and there is a need to implement screening approaches to identify the most performing ones. Theoretical studies on molecular electrocatalysts are adequate for this purpose, since density functional theory (DFT) has a strong predicting capability to provide clues for the improvement of practical devices. This chapter describes the most recent theoretical methods applied to several members of the clathrochelate family. We describe the computation of their common spectroscopic and electrochemical properties. In addition, DFT analysis is used to decipher the multistep reaction mechanism of a model Co clathrochelate with regard to the hydrogen evolution reaction in the homogeneous phase.


Introduction
Clathrochelate is a term used in coordination chemistry to denote the ligands that can encapsulate metal ions. The word is specifically used here as a generic term to designate a family of macrocyclic organometallic compounds that contain tris-glyoximate ligands, forming a cage that is strongly bound to a central metal atom (typically Fe [1] or Co [2]; (Figure 1). Interestingly, such clathrochelates form a wide family of compounds since their apical and ribbed substituents may be varied with ease. Electron donor (such as alkyl) or electron-withdrawing substituents (such as halogens) can be used in the structure to adjust the electronic properties of the metallic center. The XRD structure of the cobalt clathrochelates in their resting state (i.e., neutral molecule with the metal ion in a 2+ oxidation state) shows four N atoms of two glyoximate groups closer to Co, typically at a distance between 1.88 to 1.99 Å, and two other N atoms of the third glyoximate group more elongated at a distance of 2.10 to 2.17 Å [2,3], reflecting a Jahn-Teller distortion. Another structural parameter employed for their description is the dihedral angle (θ) between the atoms B'N'NB (where the apostrophe distinguishes the two apical extremes of the complex), which expresses the rotation between the two N′··N′··N′ and N··N··N triangles, as illustrated in Figure 1. Both Co-N distances and θ depend on the oxidation state of the molecule, reflecting structural variations occurring upon electron gain or loss.
The relevant property of these clathrochelates that motivates this chapter is their electrocatalytic activity for hydrogen production from protons in the homogeneous phase (i.e., when the compound is dissolved in solution) [2] known as the hydrogen evolution reaction (HER). The versatility of these molecules makes them further suitable for their incorporation in heterogeneous systems for light harvesting [4] or as precursors for the deposition of Co-containing nanoparticles [5], which are also active for the HER. Herein, we illustrate the use of DFT to study the electrochemical properties of these clathrochelates with regard to the HER. First, we provide a concise overview of the theoretical methods currently employed to address these complexes and the typical procedures employed to aid and stabilize self-consistent field (SCF) and geometry optimization convergence. Subsequently, we provide examples of the computation of common spectroscopies, such as infrared (IR) or ultraviolet-visible (UV-Vis). The core of this chapter is then devoted to the computation of electron transfer processes and to unravel the HER electrocatalytic mechanism with a model Co clathrochelate. Finally, we provide an overall conclusion of the state of the art regarding the application of DFT methods to clathrochelates, along with conceivable perspectives for future work on the topic.

Comments on the selection of the calculation methodology
The choice of the theoretical method is the cornerstone requirement for a successful computation leading to results of predictive value. However, there is no established and unique approach to select the appropriate method for solving a particular problem. Despite this, some general guidelines may be set in order to select the suitable tools.
For example, DFT studies of organic molecules have overwhelmingly been performed with the B3LYP exchange-correlation functional [6], which is currently a standard procedure joined to Pople's basis sets (the family of 6-31G and so on). At the present time, there are no huge problems with the theoretical treatment of relatively small organic molecules. Conversely, organometallic complexes in general, and clathrochelates in particular, are much more challenging, because of the plethora of functionals and basis sets available. Ultimately, it is the comparison with an experimental quantity that tells whether a theoretical method is satisfactory or not. Typically, the initial manner is to skim the literature to find similar systems which have been treated before and to use such methodology. In the case of Co clathrochelates, the literature is dominated by the use of the B3LYP exchange-correlation functional. We have successfully used the B3LYP functional for geometry optimizations, with excellent agreement to experiment [2]; however, all attempts to use other functionals (M06-L, TPSS, BP86, PBE, or PBE0) yielded unreasonable geometries [7]. Regarding the selection of the basis set, the LANL2DZ with effective core potentials for heavy atoms predicts well the geometry of Co clathrochelates; on the other hand, the use of a more complete and more expensive all-electron basis set cc-pVTZ for Co and cc-pVDZ for all other atoms did not improve significantly the results [7].
Regarding geometry optimization convergence, it is reasonable to be generous on the allowed iteration steps; for instance, we have always employed 100 iterations for these kind of problems. Another frequent issue is SCF convergence. We recommend the use of Pulay's direct inversion in the iterative subspace (DIIS) and the damping of the Fock matrix, with the purpose to accelerate and stabilize convergence; alternatively to DIIS, the second-order SCF orbital optimization (SOSCF) method may be used as well. Finally, the computation of SCF may be performed directly or not. This has an impact mostly on the calculation time, and the choice depends on the hardware available for work. In the case of direct SCF, the integrals are calculated each time whenever necessary; on the other hand, if direct SCF is not chosen, integrals are calculated only once and then stored in memory. For computers that have mechanical-based hard disks, it is advisable to use direct SCF, because the calculation of integrals by the processor is generally faster than their search in memory. The time spent for the calculation of clathrochelates may be high, of course depending on the functional, the basis set, and the substituents of the molecule. As an illustrative example, by running a parallelized version of GAMESS [8] on 29 processors, it took 1 week for a complete optimization and Hessian computation of complex CoBd 3 (B-nC 4 H 9 ) 2 which possesses n-butyl and phenyl substituents in the apical and ribbed positions, respectively, using the B3LYP functional and the LANL2DZ basis set with effective core potentials for Co.

Computation of typical spectroscopic properties of a model clathrochelate
A major goal of theoretical chemistry is to be able to predict measurable properties and therefore to provide insights into the behavior of molecular systems. Two of the most ubiquitous techniques available are infrared (IR) and ultraviolet-visible (UV-Vis) spectroscopies, which provide information on the normal modes of vibration (IR) and electronic excitations (UV-Vis). At the end of any geometry optimization, a Hessian calculation should be performed, giving access to the normal modes of the molecule; indeed, all calculated frequencies must be positive in order to verify that the identified stationary point is a true minimum. A typical example of Hessian calculation for the Co(Cl 2 Gm) 3 (B-PhF 5 ) 2 molecule (where Gm and PhF 5 stand for the glyoximate and the pentafluorophenyl moieties, respectively) is depicted in Figure 2a, along with a comparison to the experimental IR spectrum. The calculation was performed with the B3LYP exchange-correlation functional and using the LANL2DZ basis set for all atoms, including effective core potentials for Co and Cl. Note the excellent agreement between experiment and calculation, showing that most signals observed in the experimental spectrum are not pure normal modes but a complex convolution of multiple vibrations that are close in frequency. In particular, the absorption band highlighted with a (1) is essentially composed by a single normal mode. However, such mode is not located in a small part of the molecule but comprises the whole molecular structure, with coupled vibrations from the ligand cage, and the pentafluorophenyl rings, as depicted on the right side of Figure 2b.
Another important property of these molecules that may be described is their optical absorption. In order to calculate the UV-Vis characteristics, time-dependent density functional theory (TDDFT) formalism has to be employed. The use of the B3LYP functional for the calculation of these kind of spectra has been shown to be an adequate choice [3].
However, the basis function should be chosen to be more complete than for geometry optimizations and Hessian calculations. In particular, it is desirable to include diffuse functions in the basis set, because the optical excitations may occur to nonbonding, unfilled orbitals which develop farther away from the molecule. Previous predictions of optical spectroscopy of clathrochelates have employed a mixture of basis sets, such as the triple-ζ with a polarization function (TZVP) for Co and double-ζ split valence for all other atoms [3]. For our calculations (Figure 2c and d), we have employed the augmented version of the correlationconsistent polarized valence triple-ζ (aug-cc-pVTZ) basis set for Co and the correlation-consistent polarized valence double-ζ (cc-pVDZ) for all other atoms, and we have computed the first 40 excited states. It should be noted that, for the calculation of UV-Vis spectra of usual small organic molecules, it is often sufficient to calculate only four excited states [6], in open contrast to organometallic compounds, where a large number has to be computed in order to adequately reproduce the spectrum. As for IR spectra, these complex molecules exhibit complex UV-Vis spectra, and the assignment of the different contributions is not straightforward because each transition is composed of several contributions from different orbitals. In particular, two transitions occurring in the visible region are highlighted with numbers

Estimation of redox potentials with DFT
DFT calculations can also be used to analyze the redox properties of molecules. Electron transfer reactions may take place between chemical species dissolved in solution or between a molecule in solution and a solid-state electrode or photoelectrode. In particular, electrochemistry is the science that studies matter transformation upon the passage of electric current. Species that gain electrons decrease their oxidation state (and are said to be reduced), while species providing electrons to others increase their oxidation sate (and are said to be oxidized). These oxidation-reduction reactions, which necessarily involve the exchange of electrons but may involve chemical steps as well, are referred with the shorthand term "redox." Certain molecules may exist in two (or more) oxidation states, and these pairs are known as redox couples, represented as O|R, where O stands for the oxidized form, and R represents the reduced form. The redox reaction O + e − = R occurs at a certain electrode potential (E°), which is mainly observable in electrochemistry. The redox potential is important because the more positive it is, the more oxidizing is O; likewise, the more negative it is, the more reducing is R. Essentially, the redox potential is a thermodynamic quantity that allows to predict the feasibility of redox reactions. Absolute redox potentials for molecules are not measurable, and therefore a proper reference electrode is employed. Arbitrarily, the standard hydrogen electrode (SHE) has been adopted as zero in the potential scale [9]. This relative scale is quantitatively related to the absolute energy scale used in physics, i.e., the electron in a vacuum [10].
DFT may provide an estimation of the absolute redox potential of molecules, following two main methodologies. The first method consists in the calculation of a Born-Haber process, by optimizing all relevant species both in gas phase and in the solvent. This procedure is relatively lengthy and therefore shall not receive further attention in this chapter, but the interested reader is referred to the specialized literature [11][12][13][14]. Conversely, the full calculation of the Born-Haber process may be bypassed, and the optimization and Hessian calculation of the relevant redox species, i.e., O and R, in the desired solvent may be undertaken directly [7,15]. The calculation of the absolute redox potential of the species is found according to Eq. (1): where F and n represent Faraday's constant and the number of exchanged electrons, respectively, and ΔG abs 0 is the difference in free energy of the reaction O + e − = R. The value of ΔG abs 0 contains (i) the total electronic energy of the molecule (i.e., electron-electron and nucleus-nucleus repulsions plus electron-nucleus attraction), (ii) zero-point energy (ZPE) corrections, (iii) vibrational corrections to the free energy, and (iv) thermal corrections. We note that the electronic energy of the molecule is obtained by a single-point calculation, but all other thermodynamic corrections are only accessible after a Hessian calculation.
This calculated absolute potential may be referred to an estimated value of the absolute potential of the SHE. Unfortunately, the literature is plenty of different values, for example, 4.2 ± 0.4 V [16] to 4.29 ± 0.02 V [17], 4.05 V [18], and 4.44 V [19], which makes it difficult to assess precision in the calculation, given that all authors provide arguments for their absolute values of the SHE. In order to overcome this drawback, the choice of a reference redox reaction is a safer option. The most commonly used reference reaction is that of the Fc 1+ |Fc (i.e., ferrocenium|ferrocene) redox couple. When this procedure is used, the optimization and Hessian calculation of both Fc 1+ and Fc at exactly the same theoretical level than the O|R couple under study are performed (Figure 3a). Then, the absolute potential of Fc 1+ |Fc is obtained using Eq. (1). Obviously, this is a sort of "redox calibration" of the computational method. The theoretical estimation of the redox potential of the O|R couple vs. the Fc 1+ |Fc couple is just the difference between the absolute potentials these two redox pairs. Current DFT methods provide estimates of redox potentials with an error of at most 200 mV [11,12,20], which is nowadays tacitly considered as the limit for a reasonable estimation. There are several causes for such discrepancies with experimental observations, such as the fact that (i) the exact functional linking the electron density and the electronic energy is not known, (ii) the use of an incomplete basis set, (iii) the consideration of only one conformation of the molecule for the computation, (iv) the treatment of the solvent as an electrostatic continuum, and (v) the harmonic oscillator approximation to estimate thermodynamic quantities. Despite these sources of errors, DFT is powerful in predicting reactivity trends while giving reasonable estimations for observable properties.  The comparison between theoretical and experimental redox potentials has another subtlety. Indeed, when an aqueous reference electrode, such as the saturated calomel electrode (SCE) or the silver-silver chloride (Ag|AgCl|Cl − ), is used in a nonaqueous environment, there is an inherent liquid junction potential (LJP) [21] which depends on electrolyte composition and biases the measurement. Hence, in order to compare theoretical results with experimental values, the first issue to pay attention to is that an aqueous reference had not been used to make the measurement in a nonaqueous environment. If an aqueous reference has been used in a nonaqueous solvent the Co atom moves to one side, it is still possible to make the comparison with theory, by finding a literature value for the Fc 1+ |Fc couple and converting the potential scale of the aqueous reference (which was used in a nonaqueous solvent) to the Fc 1+ |Fc scale. Literature values for Fc 1+ |Fc have been reported as +0.50 V vs. Ag|AgCl|Cl − [22] or + 0.46 V vs. SCE [2], comprising the LJP.
Finally, it is worth to stress that the redox potential and the energy of molecular orbitals are different things and cannot be equated. Certainly, when representing electron transfer reactions, the educational literature often depicts electrons as being transferred from frontier orbitals [23]. However, such representations might induce to think that the energy of orbitals is somehow linked to redox potentials. As a matter of fact, molecular orbitals are mathematical functions that describe the electron motion in a molecule and are certainly not observable [24], in contrast to the redox potential which is an observable thermodynamic quantity. Despite this, the energetic trend in orbital energies may be commonly paralleled to the oxidizing/reducing power of molecules. For example, with all other things being equal, the more electronegative the substituents in the ligand of a complex, the lower the energy of the molecular orbitals. Hence, the complex will be more likely to gain an electron in the lowest unoccupied molecular orbital (LUMO) than to give it from the highest occupied molecular orbital (HOMO); this is equivalent to say that the complex shall be a better oxidant and its redox potential should be more positive. This kind of trends, linking the energy of molecular orbitals and redox potentials, has been presented several times in the literature [25,26]. Though we have shown that the Co III |Co II couple in a series of clathrochelates complies with it, the Co II |Co I redox pair does not follow the orbital trends. Further details may be found in Ref. [2].

Mechanistic investigation of the electrocatalytic hydrogen evolution reaction and ligand noninnocence
The origin of the catalytic activity of metal-containing clathrochelates is currently under debate [1]. It is likely that the role of the ligand is essential to electrocatalysis, since it may stabilize low oxidation states in the metal to enable electron transfer or may participate in protonation reactions that may help bring protons together and evolve H 2 . That is why the term "ligand non-innocence" has been coined to refer to organometallic complexes in which the organic ligand has a major participation in its redox behavior to account for high reactivity involving small species, such as CO 2 or H + [27][28][29].
For example, upon the one-electron reduction process depicted in Figure 3b, we have predicted interesting structural modifications within the molecule. The asymmetry inside the Co-N cage is calculated to be reinforced, with four Co-N distances being at around 1.88 Å and the two other N lying farther apart at 2.45 Å; consequently, the Co atom moves to one side of the cage, thus decoordinating from two N atoms and resulting in a square planar geometry for the reduced complex; furthermore, there is a slight shortening of N-C(sp 2 ) bonds from 1.30 to 1.28 Å.
Ligand non-innocence in the case of clathrochelates may be further investigated because their chemical structure invites to think on multiple sites as putative candidates for protonation. Previous theoretical studies have addressed only protonation on the iminic C(sp 2 ) [30], but we have extended such approach by further exploring other protonation sites. Indeed, the basic character of N and the wide evidence of metal hydrides as intermediates in the HER have led us to explore such alternative protonation sites.
Until now, the most complete theoretical study of hydrogen evolution involving clathrochelates is found in Ref. [7]. The work was performed on a clathrochelate having chlorine and methyl moieties as the ribbed and apical substituents, respectively. Henceforward, this epigraph presents those results in a concise manner.
Optimized structures for three protonation possibilities are illustrated in Figure 4; the investigated protonation sites are highlighted with a red square (C), blue circle (N), and an orange triangle (Co). The insets on the right side of Figure 4 show a magnification of the coordination position of H + . Protonation of the iminic C(sp 2 ) brings a hybridization change to (sp 3 ) with the concomitant change in the geometry of such C from planar to tetrahedral. Protonation of N has a similar effect on the local geometry of the ligand, which is modified from planar to tetrahedral as well. Direct protonation of Co might seem unlikely at the first sight because the metal appears buried inside the organic cage. However, the movement of Co inside the  ligand cage upon the one-electron reductive activation and the associated decoordination from two N makes the metal accessible for protonation as well. The calculated X-H distances (where X stands for Co, N, or Csp 2 ) for these three species are 1.09 Å for Csp 2 -H, 1.45 Å for Co-H, and 1.03 Å for N-H. Such bond distances show that proton binding to N or C results in the formation of a strong bond, in contrast to the binding to the metal center, where a much more elongated Co(III)-H bond is created, indicating a weaker binding and therefore a more labile H + .
An energetic analysis of this initial part of the mechanism is illustrated in Figure 5a. Indeed, there is an energetic gain upon the reductive activation of the [ 2 CoL] 0 to yield the catalytically active [CoL] 1− intermediate. Then, protonation of N or C(sp 2 ) is likely, but the most favored protonation site was predicted to be the Co central ion.
We have also addressed the possibility of double protonation in two adjacent N; in a C(sp 2 ) and a N; in two adjacent C(sp 2 ); in Co and N; and in Co, for both high-spin and low-spin species. However, most optimized structures were largely distorted due to heavy structural changes caused by the sp 2 to sp 3 hybridization changes that should occur after protonation. As a consequence, the calculated energy of doubly protonated intermediates was prohibitively high, and these species are unlikely to be formed. The most stable doubly protonated species was the one with protons coordinated in two adjacent N atoms. However, it seems unlikely that the mechanism of hydrogen evolution involved initially a mono-protonated species that is a cobalt hydride and that the subsequent species in the mechanism is a double-protonated intermediate with protons in adjacent N atoms.
The Co(III)-H intermediate is interesting and deserved closer attention. In particular, such species may undergo a further one-electron reduction to yield Co(II)-H, which may also have a role in the mechanism. Certainly, we have predicted that such reduction brings an important energetic gain, but the search for doubly protonated intermediates and the search for a transition state to explain hydrogen evolution were unsuccessful. We were able to locate a transition state of the η 2 -Co(III)-(HH) type, but again this was far too energetic to be a plausible candidate to explain hydrogen evolution (Figure 6b). Considering that the hydrogen production step is of key importance, we have devoted wide efforts to explain how it happens. When the calculations on a mechanism yield high energy intermediates, it is possible that a concerted electron-proton transfer process occurred, in order to avoid energy-demanding paths [15,20]. Thus, potential energy surface (PES) scans were performed with the Co(II)-(HH) and Co(III)-HH) species at different (fixed) Co-H distances, while the rest of the structure was allowed to optimize freely. These results are summarized in Figure 6a and b. Interestingly, in the case of Co(II)-(HH), the increase in the Co-H bond distance was linked to a decrease in the H-H distance, approaching the equilibrium bond distance for H 2 , which is 0.74 Å. At the same time, the energy of the system gradually decreased to approach the final energy level of the reaction (observe the blue scale on the right side of Figure 6a), which is depicted in Figure 5b as [ 2 CoL] 0 + H 2 . These calculations show a nice picture of two protons that come together in the nearby of Co(II) and bind to one another to form H 2 while going away from the metal, with the concomitant decrease in total energy. Conversely, the same calculations on the Co(III) species gave the same bond-breaking (Co-H) bond-formation (H-H) sequence as just described for Co(II). However, at the same Co-H distance, the H-H separation was consistently higher in the case of Co(III), thus indicating that Co(III) was less effective than Co(II) to promote the reduction and binding of two protons. Moreover, the total energy of the system for the Co(III) ion did not decrease smoothly as in the case of Co(II) (Figure 6b). Instead, it was excessively high, and therefore such intermediate was definitively ruled out as a possibility to explain the catalytic activity of this clathrochelate.

Conclusions and perspectives
The clathrochelate family has received attention due to its multiple applications. In this chapter, we have illustrated the use of density functional theory to describe several properties of clathrochelates. We began the presentation by providing tips on the choice of the calculation methodology to address these complex molecules. Then, the spectroscopic properties (IR and UV-Vis) have been reproduced. The possibility to accurately predict spectroscopic properties gives an adequate degree of confidence on the theoretical procedures employed. Subsequently, we have illustrated the procedures to predict the redox potentials of molecules, and we have highlighted some structural changes underwent during the electron transfer process. This chapter has introduced the concept of "ligand non-innocence," which is relevant in complex electrocatalytic reactions involving complexes of transition metals, where the ligand has a major role to play, not only to stabilize unusually low oxidation states but also to act cooperatively to bind small species, such as H + . Finally, a methodology has been presented to address the complex mechanism of hydrogen evolution. Certainly, the hydrogen production step was particularly intriguing for the Co(Cl 2 Gm) 3 (B-CH 3 ) 2 complex, given that neither doubly protonated intermediates nor plausible transition states could be found. Therefore, a concerted proton-electron transfer step was envisaged, with the Co(II)-(HH) hydride intermediate playing the most relevant role. The calculation of potential energy surface scans fixing the Co(II)-(HH) bond distance showed that the simultaneous elongation of Co-H bond leads to the decrease in H-H distance and the smooth decrease in energy, which nicely explained hydrogen evolution from this complex in the homogeneous phase. Overall, results reported in this chapter could contribute to stimulate theoretical and experimentalist chemists to explore the influence of different substituents in the organic cage on the HER activity, thereby providing new insights to researchers that try to develop and optimize alternative electrocatalysts for the HER. Clearly, the formal theoretical understanding of clathrochelate's electrochemistry is still in its early stages, and the missing piece in this scientific challenge is to apply the theoretical methodology outlined in this chapter to other clathrochelates in order to justify the distinct reactivity observed experimentally among different complexes. Moreover, the approach proposed here should provide guidelines to synthetic chemists in the search for more active compounds. Thus, we foresee that the future theoretical research shall be performed in such direction. This is a stimulating moment to conduct research aiming to find novel and efficient catalysts for hydrogen production, and transition metal clathrochelates are perhaps the family of complexes that has received the least attention from the chemical community.