Alternative Pathway of CO2 Hydrogenation by Lewis-Pair-Functionalized UiO-66 MOF Revealed by Metadynamics Simulations.

The reaction between H2 and CO2 catalyzed by an intramolecular frustrated Lewis pair, which is covalently bonded to a UiO-66 metal-organic framework (MOF), is considered in this work. Free energy surfaces (FESs) for this reaction are generated throughout finite-temperature density functional theory (DFT) metadynamics (MD) simulations. The simulated FESs indicate an alternative stepwise pathway for the hydrogenation of CO2. Furthermore, indications of weaker binding of CO2 than H2 to the Lewis pair centers have been observed via metadynamics simulations. These findings were unknown from the results of static-DFT calculations, which proposed a concerted reduction of CO2. The results of the present work may influence the design of new efficient heterogeneous Lewis pair (LP)-functionalized MOFs to catalyze capture and conversion of CO2 to high-value chemicals.

Introduction

The removal of the byproducts of traditional oil-based fuel usage, such as CO2, from the environment, is one of the several challenging problems that require new technological solutions. One approach to reducing both carbon emission and consumption of fossil fuels is CO2 capture coupled with chemical recycling of CO2 to renewable fuels and valuable chemicals, e.g., CH3OH and CH2O.[1−4]

For a long period of time, it was thought that transition metals were the only viable materials for the adsorption of gases and catalysis of their conversion reactions. This paradigm changed with a 2006 report by Stephan et al. of metal-free molecules that reversibly activated small molecules by systems described as “frustrated Lewis pairs” (FLPs).[5] The combination of a Lewis acid (LA) and a Lewis base (LB) in such a way that they are prevented from forming a dative bond and a classical LA–LB adduct by means of steric hindrance or thermal hindrance is called an FLP. The unique reactivity of FLPs continues to find applications in new areas, such as transition-metal chemistry, bioinorganic chemistry, and materials science.[6] In the past decade, commercial interest in FLP reactivity has been triggered by its lower toxicity, distinct functional group tolerances, and reduced catalyst and product purification costs.

However, the implementation of such metal-free catalysts faces the same concerns as other homogeneous catalysts: stability, recyclability, and catalyst–product separation. It is for this reason that the notion of exploiting solid/heterogeneous FLPs to support FLP catalysts is an attractive topic that has been proposed and examined recently.[7,8] To this end, employing Lewis-pair-functionalized metal–organic frameworks (MOFs), in which the Lewis pairs are covalently bonded to the MOF as a solid platform to fix Lewis pair centers, has been recently proposed and experimentally utilized to capture, convert, and reduce multiple bonds of small molecules.[9−14]

In 2015, Lewis pair (LP)-functionalized UiO-66 was proposed by Johnson and co-workers to develop novel nanoporous materials for CO2 capture and conversion from flue gas. Their static-DFT calculations provided design strategies for efficient catalysts for CO2 reduction. The resulting functionalized MOFs are very stable and can have enhanced catalytic activity due to the stabilization or protection of the catalytic complexes.[9] However, a lack of computed free energy surfaces (FESs) of the reaction between H2 and CO2 to understand the details of the dynamics of the reaction path motivated us for the present work.

In the present work, we investigate the reaction of CO2 + H2 ⇌ HCOOH, reaction 1, inside the cavity of a Lewis pair (LP)-functionalized UiO-66 MOF. The recent findings obtained by advanced ab initio molecular dynamics (AIMD) simulations and the possibility of alternative pathways for the hydrogenation of carbonyl compounds suggest that the nature of complex reactions, e.g., FLP-mediated molecular activation, could be challenged in a dynamic picture.[15−26] This motivated us to investigate reaction 1 using the AIMD methodology, for the first time.

By employing metadynamics simulations, we calculate the FES of the reaction 1 at finite temperature (325 K). Simulation of a free energy landscape rather than a zero kelvin potential energy surface allows us to incorporate entropic contributions due to the flexibility of motion of free small gaseous molecules in the presence of LP-functionalized MOF complexes.

Figure shows the optimized structure of the unit cell of UiO-66 with the singled out benzene dicarboxylate (BDC) linker (I) and its LP-functionalized primitive (II). The structure of intramolecular LP (III) with N and B centers that is connected with a methyl group to the BDC linker is shown in Figure . The role of the selected MOF is to act as a scaffold and fix the Lewis pairs inside the periodic MOF structure. We note that the MOF atoms do not interfere in the reaction mechanism. It is important to mention that the geometric constraint of the LPs being covalently bonded to the MOF at specific sites provides the Lewis pairs without steric hindrance and prevents the mutual quenching of LP moieties.

Figure 1

Optimized structure of the unit cell of UiO-66 with the benzene dicarboxylate (BDC) linker singled out with a circle (I) and its LP-functionalized primitive (II). (III) Structure of an intramolecular Lewis pair covalently bonded to the BDC via a CH2 group.

Utilizing MOFs as a solid platform for Lewis pairs is a practical way, according to experiments, to use the catalytic functionality of LP centers in the solid/heterogeneous state to overcome the drawbacks of homogeneous FLP catalysis, i.e., stability, recyclability, and catalyst–product separation. Hence, binding the LP centers to the BDC linkers of UiO-66 by a methyl group prevents the migration and association of the LP centers. In this work, we have selected nitrogen as the Lewis base (LB) and boron as the Lewis acid (LA) (Figure ). First, we investigate the energetics of adsorption/desorption of H2 and CO2 on the LP centers and then analyze the energetics of reduction of CO2 by the dissociated H···H on the LP centers. The lattice constants of the optimized primitive cell are a = b = c = 14.788 Å and α = β = γ = 60°. As indicated in previous studies by Johnson and co-workers, fully relaxing the geometry and cell parameters of LP-functionalized UiO-66 gave lattice constants that were almost identical to the relaxed UiO-66 values. Furthermore, optimizing the structure with chemisorbed CO2 or H2 in UiO-66-LP also perturbed the lattice constants and energies by a very minor amount.[9,27−30] Therefore, the lattice constants were fixed at the ground-state UiO-66 values to save the computational time. Our DFT-MD simulations indicate that the porous framework remains stable after the functionalization and chemisorption of H2 and CO2 and catalytic moieties can be covalently bonded into UiO-based MOFs without losing their catalytic activity. More importantly, the possibility of an alternative pathway for the hydrogenation of CO2 starting from dissociated H···H on the LP centers is indicated. This work presents a step toward the discovery of practical heterogeneous metal-free catalysis for CO2 reduction and conversion.

Computational Details

All periodic DFT calculations—including geometry optimizations and ab initio molecular dynamics (AIMD) simulations—were performed using the CP2K program[31] with the Gaussian and planewave (GPW) method. The valence orbitals were expanded in the DZVP-MOLOPT Gaussian basis set in combination with Goedecker, Teter, and Hutter pseudopotentials with a planewave cutoff energy of 280 Ry. We used the PBE density functional[32]augmented with the Grimme D3 dispersion correction.[33] The criterion for the self-consistent field convergence was set to 5.0 × 10–7. The unit cell of UiO-66 was taken from the Cambridge Crystallographic Data Center, and the calculations were performed on the primitive cell obtained by the RASPA program to save computational efforts.[34] This means that the P1 cell was obtained from the CIF file listing only unique atoms with symmetry operations, and the full unit cell was reduced to the primitive cell using RASPA. The AIMD simulations were performed in the NVT ensemble, with the temperature controlled by a CSVR (canonical sampling through velocity rescaling) thermostat, set at a temperature of 325 K and a period of 500 fs. The MD time step was 0.5 fs. To investigate the reaction mechanism at finite temperature, characterize the reaction pathway, and identify the transition-state (TS) region between the reactant and product states, we performed metadynamics simulations using the PLUMED plug-in in combination with CP2K.[35−37] The collective variables (CVs) were the relevant distances in each step of the reaction, i.e., B···H, N···H, and H···H for H2 dissociation/formation; B···O(CO2) and N···C(CO2) for CO2 adsorption/desorption; and B···H–, H–···C(CO2), and H+···O(CO2) for CO2 hydrogenation. The simulations were initiated from the reactant states, i.e., free H2 and CO2, and conducted until several transition-state recrossing events could be observed. To prevent sampling from the unnecessary regions of the FES, harmonic walls with a force constant K = 250 kcal mol–1 Å–2 were used for CVs. The Gaussian bias potentials were added every 25 steps until the first conversion happens in each reaction step and then decreased to every 100 steps till 20 ps simulations. Multiple test simulations were run to set the computational parameters, including the number and type of CVs, height and width of Gaussian bias potentials, and quadratic walls for each part of the reaction mechanism. The accuracy of the calculated free energy values is 1.0 × 10–9 kcal/mol. The computational details of each reaction step are explained in the Supporting Information (SI). Using the trace_irc program structures of the minima, the transition states were refined.[38]

Results and Discussion

According to the static-DFT calculations reported in the literature, the reduction of CO2 takes place at the heterolytically dissociated H···H that are bonded at the LP centers, which is found to be energetically lower than the alternative pathway of chemisorbed CO2 at the LP centers followed by H2 dissociation. Hence, in the present work, we focus on the former pathway. As a matter of fact, binding the CO2 to the LP centers before H2 deactivates the catalytic sites and disables them for H2 dissociation. Due to the structural properties of LP covalently bonded to the BDC linkers, the LA–LB distance is fixed. We simulate the FES of the CO2 hydrogenation starting from free H2 dissociation on the LP centers and in each step analyze the dynamic behavior of the free or adsorbed molecules and energetics. For the sake of energetic comparison, the CO2 adsorption on the LP centers is also investigated.

H2 Adsorption–Desorption on the LP Centers

We started the ab initio, DFT-based molecular dynamics (DFT-MD) simulations with the free H2 molecule in front of the Lewis pair (LP) centers, i.e., N and B atoms, covalently bonded to the UiO-66 BDC linker. During the first test runs of dissociated H2 (i.e., N–H and B–H bonds are already formed) on the LP centers, we observed that the HNBH dihedral angle is very flexible and this means that the H–B bond freely rotates around the plane of NNH atoms (Figure ). By this rotation, the conformation around this dihedral varies from initial cis (N–H and B–H in the same direction) to trans (N–H and B–H in the opposite direction). In the trans conformation, H···H atoms are not close enough to form the H2 molecule and release it. Hence, in metadynamics simulations of molecular H2 dissociation/formation in addition to N···H and B···H distances, we used a third CV that is the H···H distance and put a wall at 2.5 Å to maintain the two H atoms close enough during simulations. We have run 20 ps metadynamics simulations, where Figure shows the variation of N–H, B–H, and H–H distances (Å), which are the three CVs in this metadynamics simulations, versus simulation time (ps). In Figure , the entire trajectory has been divided into seven representative regions denoted I, II, III, IV, V, VI, and VII. We discuss hereafter each of the regions step by step along the simulation time. Region I is basically the motion of the free H2 close to the LP centers at room temperature. There are some fluctuations and variations in the distances that cause minor elongation in the H–H bond length and decrease N···H and B···H distances to 1.5 Å but no bond formation or breaking can be observed yet since the H2 is either close to N or B. At the end of region I and before H2 dissociation, the B···H has larger fluctuations than N···H and is shorter.

Figure 2

Two conformations around the HNBH dihedral. In the trans conformation, the H···H distance is farther than what is required for H2 formation. For a clearer picture, we have removed the UiO-66 primitive in the figure.

Figure 3

Variation of N–H, B–H, and H–H distances (Å) that are three CVs in metadynamics simulations of H2 dissociation versus simulation time (ps).

In region II, the H2 bond dissociates and the N–H and B–H bonds form. The N–H and B–H bonds vibrate around their equilibrium values, i.e., 1.0 and 1.2 Å, respectively. However, the B–H bond vibrates with a larger amplitude than the N–H bond. We note that the MTD parameters were the same for both N···H and B···H. In region III, 5.0–6.5 ps, both molecular and dissociated H2 can be observed. This region indicates an equilibrium between free and dissociated H2, i.e., a rather large elongation in the N–H and B–H distances leads to H2 formation. Region III is an interface between the free H2 molecule and the dissociated one in which N and B strongly interact with H–H. Region IV is a long interval between 6.5 and ca. 12 ps, in which the H2 molecule forms and vibrates around the equilibrium bond distance, 0.74 Å. There are large fluctuations in N···H and B···H distances within 1.3 and 2.5 Å, and at ca. 11 ps, the strong interactions between N and B with H–H cause the dissociation of the H–H bond. However, N···H and B···H distances behave often oppositely (one elongated and the other shortened) in region IV and consequently H2 stays in the molecular form. In most of the states, H–H is closer to B. This indicates an asynchronous interaction (formation) of B–H and N–H bonds. In region V, the strong interaction between LP centers and H–H leads to the dissociation of H2 molecules (B–H and N–H bonds are formed) for a period of 1 ps. In region VI for a long period between ca. 13 and 18.8 ps, free H2 forms but still there is a strong interaction between LP centers and H–H. In region VII, H2 dissociates and N–H and B–H bonds form; however, there are large fluctuations in bond distances, which indicate an equilibrium state between free and dissociated H2; the formation of H2 at some instances can be observed.

For a closer insight into the behavior of N···H and B···H versus H···H distances, in Figure , upper panel, we have shown the variation of H···H (vertical axis) versus N···H and B···H (horizontal axis) distances. Figure , upper panel, shows that both N···H and B···H have similar behavior versus H···H, i.e., for both N···H and B···H distances within an interval of 2.6–1.5 Å, H···H oscillates around the equilibrium H–H distance (0.74) and for N···H and B···H distances less than 1.5 Å, the H···H sharply increases. Figure , middle panel, shows the variation of the B···H versus N···H distances. In addition to the reactant (free H2) and product (formed N–H and B–H bonds) area, three specific regions with a higher population of states can be detected. These three regions have been singled out with the blue and red rectangles and a gray circle. In the blue rectangle, the B···H distance mainly fluctuates around 2.5 Å and the N···H distance decreases up to 1.5 Å. On the other hand, in the red rectangle, the N···H distance fluctuates around 2.5 Å and the B···H decreases up to 1.5 Å. Hence, the hydrogen molecule is initially polarized either by N or B and then moves toward the other center. As soon as both distances are shortened enough, i.e., 1.6–1.8 Å, singled out with a gray circle in Figure , middle panel, the H···H passes through the transition-state region. This region consists of structures with both N···H and B···H distances shortened, i.e., prior to the transition-state area. The reaction proceeds in such a way that the polarized H···H (with an elongated distance) is either close to the N center or nearby the B center, the regions singled out by blue and red boxes, respectively. Two representative structures of each region are shown in the lower panel of Figure . The side-on arrangement close by a boron center and the end-on arrangement nearby the N center can be observed.

Figure 4

Upper panel shows the variation of H···H (vertical) versus N···H and B···H (horizontal) distances. The middle panel shows B···H versus N···H. The lower panel shows two selected structures that the polarized H···H is close to the B center (side on) in the red box and the polarized H···H is close to the N center (end on) in the blue box.

Figure , upper part, shows the variation of free energy versus N–H, B–H, and H–H distances. The three-dimensional free energy plot obtained from metadynamics simulations shows two minima: the two components complex H2 + LP-MOF (A in Figure ) and the hydrogenated N–H and B–H bonds covalently bonded to MOF (C in Figure ) and the transition-state structure (B in Figure ). In Figure , upper part, the x, y, and z axes are N···H, B···H, and H···H distances (CVs) in Å, respectively. The four layers with different colors in the free energy plot show the increase in free energy.

Figure 5

Upper part shows a three-dimensional (3D) variation of free energy versus N–H, B–H, and H–H distances (CVs), in xyz axes, respectively. The energy values (numbers) are in kcal/mol. The four layers with different colors in the free energy plot show an increase in free energy. The structures of reactant, transition state, and product complexes with selected distances are shown in the lower part. For a clearer picture, we have removed the UiO-66 primitive in the figure.

The barrier of the H2 dissociation catalyzed by LP-functionalized UiO-66, in this case, is 18.8 kcal/mol, and the dissociation is exergonic with a ΔG of reaction equal to −6 kcal/mol. The barrier of recombination of the dissociated H···H is 24.7 kcal/mol. It is worth noting that the barrier of dissociation of the free H2 molecule by the static-DFT method was calculated as 11.6 kcal/mol, according to the literature.[9] However, metadynamics simulations in this work predict a ca. 6.0 kcal/mol higher barrier for H2 dissociation. This may be due to the movements of the free H2 molecule in the proximity of the LP centers that increases the entropic contributions and decreases the ΔG of the reactant complex, which was unseen in the static-DFT calculations.

The structures of reactant, transition state, and product complexes, refined by the trace_irc program, with selected distances are shown in Figure , lower part. At the transition state (TS), the elongated H–H has a side-on arrangement to the boron with a BHH angle of ca. 100° and a NHH angle of ca. 140° with respect to the nitrogen center, which is similar to the arrangement of the adsorbed H2 molecule on a solid surface.[18] The transition-state area is a typical saddle-shaped surface with no local minima in the proximity of TS, which is consistent with a single-step H2 cleavage.

CO2 Adsorption–Desorption

To investigate the dynamics and energetics of CO2 adsorption/desorption on the LP centers and compare it with the dynamics of H2 dissociation, we have run 20 ps metadynamics simulations starting from free CO2 in front of the Lewis pair (LP) centers. The CVs in this run were the distances between the Lewis pair centers and carbon and one oxygen of CO2, i.e., N···C(CO2) and B···O(CO2) distances. The variation of the CVs versus simulation time is shown in Figure .

Figure 6

Variation of the C···N and B···O distances (CVs) for CO2 adsorption/desorption on the LP centers versus simulation time.

In Figure , the entire trajectory has been divided into seven representative regions denoted I, II, III, IV, V, VI, and VII. In region I, the motion of free CO2 near the LP centers with a small variation in N···C(CO2) and B···O(CO2) distances (2.4–3.0 Å) can be seen. In region II, the LP centers strongly interact with C(CO2) and O(CO2) and both N–C and B–O bonds form; however, there are oscillations in these bonds (stronger for B–O) that end in breaking of both bonds and desorption of CO2 in region III (6–7 ps). In region IV, the strong interaction between N and C(CO2) forms the N–C bond; however, there are large fluctuations in the B···O(CO2) distance that prohibit the B–O bond formation and only elongated B–O bond (ca. 1.8 Å) can be observed in this region. In region V, free CO2 is formed and adsorbed again on the LP centers, and in region VI, CO2 is desorbed and the interactions between N···C(CO2) and B···O(CO2) can be observed that often behave oppositely, i.e., one elongates and the other decreases and vice versa, and the motion is diffusive and unbound. In region VII, the B–O bond forms and there are strong interactions between N and C(CO2) without any real N–C bond formation.

The calculated FES of CO2 adsorption/desorption against the LP centers is shown in Figure . The barrier of CO2 adsorption on the LP centers is 10 kcal/mol, and the adsorption reaction is exergonic with a ΔG of −11.1 kcal/mol. The barrier for CO2 desorption is 21.1 kcal/mol. As indicated in Figure , the FES is dependent on both N···C and B···O distances almost equally. The structures of reactant, transition state, and product complexes, that have been refined by the trace_irc program, with selected distances are shown in the lower part. For a clearer picture, we have removed the UiO-66 primitive in the figure. A comparison between the results of 20 ps metadynamics simulations for the dissociation–formation of H2 and CO2 adsorption–desorption indicates that the states with dissociated H2 (stable N–H and B–H bonds) are more dominant than the adsorbed CO2. In the case of CO2, the states with weakly bonded C or O of CO2 to the LP centers are dominant. Moreover, the ΔG‡ of H2 release is 3.6 kcal/mol higher than that of the CO2 desorption, i.e., 24.7 vs 21.1 kcal/mol for backward reactions.

Figure 7

FES of CO2 adsorption/desorption. The energy values (numbers) are in kcal/mol. The contour lines are spaced at 0.7 kcal/mol intervals. The structures of reactant, transition state, and product complexes with the selected distances are shown in the lower part. For a clearer picture, we have removed the UiO-66 primitive in the figure.

The barrier height depends on the chemical nature, structural properties, and interactions between various fragments at the structure of the transition state. Regarding that, the analysis of these interactions is beyond the scope of the present work and will be addressed in a future study.

The three-dimensional flexibility of small free gaseous molecules such as H2 and CO2 in front of the LP centers (reactant complex), which is absent at the transition state due to the presence of stronger bonding interactions, increases the degree of freedom in the reactant state. Higher freedom of the reactant state increases the entropic contribution and lowers the ΔG of the reactant complex vs transition state. The notion that the three-dimensional flexibility of molecular pairing is connected with entropic gain has been considered by Kim and Rhee in FLP chemistry.[39]

Dynamics of the Molecular Complex with Dissociated H2-Formed N–H and B–H Bonds

Considering the flexibility of the HNBH dihedral, we found that trans-like conformations are also viable for dissociated H2, Figure . The H···H distance in the trans-like conformation is too long for H2 formation, and the dissociated H atoms cannot easily form a bond to reproduce a H2 molecule. In addition to the higher barrier of H2 recombination than CO2 desorption, the flexibility around HNBH and consequent the probability of trans-like conformations influence the reaction mechanism. This means that in the case of trans-like conformations the H···H atoms are far from each other and the reaction mechanism of CO2 hydrogenation proceeds through either a hydride transfer to C(CO2) followed by proton transfer or vice versa, leading to a stepwise mechanism. However, in the case of cis conformation, the proton and hydride are almost in the same orientation and the hydride and proton transfer can take place simultaneously, leading to a concerted hydrogenation mechanism. In the next sections, we have discussed in detail the CO2 hydrogenation mechanisms starting from both trans and cis conformations around the HNBH dihedral.

Trans Conformation—Hydride Transfer to C(CO2)

We started 20 ps metadynamics simulations with an optimized trans conformer with an HNBH dihedral angle of −120°. The H···H distance in the initial structure is 3.4 Å. Considering the long distance between two H atoms, and in particular between H+ and O(CO2), the CO2 molecule can be either close to the hydride or nearby the proton; hence, the possibility for a concerted hydride and proton transfer to the CO2 molecule in a single step is minor. So, we consider each case (proton or hydride transfer) separately throughout a stepwise mechanism. CVs in this simulation are the B···H– and H–···C(CO2) distances. In Figure , we have shown the variation of these collective variables. The bond formation between H+ and O(CO2) takes place after hydride transfer to C(CO2). Figure indicates that H–···C(CO2) (in black) fluctuates before the C–H bond completely forms at ca. 2.5 ps. The distance between each of the oxygen atoms of CO2 with B during the simulation at the shortest is 2.0 Å. This means that after hydride transfer to C(CO2) there is no evidence of B–O(CO2) bond formation, which indicates a rather fast proton transfer to O(CO2), after ca. 400 fs. During the entire 20 ps metadynamics simulations, the H···H distance remained at the shortest 2.3 Å. Hydride transfer to C(CO2) happens ca. 400 fs earlier than proton transfer to O(CO2). This means that hydride transfer catalyzes the forthcoming proton transfer to O(CO2) (Figure ). The variation of the H+···O(CO2) distance, that is not a CV in the metadynamics simulation but significantly changes during the simulation, is shown in Figure , lower panel. The variation of HNBH during 20 ps simulations shows a transformation from the trans to cis conformation at ca. 2.5 ps, i.e., HNBH changes from −120° in the initial molecular complex to 15°, Figure , structure E. The formation of the CO2 in the backward reaction is triggered by proton transfer to the N atom (LB center) and then within a short time interval hydride leaves the carbon atom.

Figure 8

Upper panel: variation of CVs, B···H–, and H–···C(CO2) distances. Lower panel: variation of the O···H+ distance and HNBH dihedral versus simulation time for the stepwise hydrogenation of CO2, starting from the trans conformation around HNBH.

Figure 9

FES of the stepwise mechanism of the hydrogenation of CO2 versus CVs, i.e., hydride transfer catalyzes proton transfer. The energy values (numbers) are in kcal/mol. The contour lines are spaced at 1.5 kcal/mol intervals. For a clearer picture, we have removed the UiO-66 primitive in the molecular structures.

Figure shows the FES of the stepwise mechanism of CO2 hydrogenation, i.e., hydride transfer catalyzing proton transfer. The FES is constructed on the hydride transfer from B Lewis acidic center to C(CO2). The horizontal axis is the B···H distance, and the vertical is the C···H distance. As Figure demonstrates the CO2 hydrogenation, the reaction starts by decreasing the C···H distance at an almost constant B···H bond length (1.2 Å). Around the transition-state region, the B···H distance is elongated accordingly. This means that the CVs behave rather asynchronously. Following the formation of the C–H bond due to the excess of negative charge on the OC(H)−O anion, it captures the proton and forms the product of hydrogenation, i.e., HCOOH. The stepwise hydrogenation of CO2, in which hydride transfer proceeds with a barrierless proton transfer, is an endergonic reaction with a ΔG of 6.0 kcal/mol and a barrier of 20.9 kcal/mol.

Trans Conformation—Proton Transfer to O(CO2)

For an investigation of the possibility of a stepwise mechanism starting from the trans conformation first with proton transfer to O(CO2) before hydride transfer to C(CO2), we have run 10 ps metadynamics simulations. For this simulation, the CVs were the N···H+ and H+···O(CO2) distances. It turned out that this is not an alternative to the stepwise mechanism since the barrier for proton transfer followed by hydride transfer is 30 kcal/mol, which is 10 kcal/mol higher than the reaction starting with hydride transfer and proceeding by proton transfer. The analysis of the atomic charges localized on the B, H–, H+, and N atoms of the B–H and N–H bonds in the dissociated form of H2 indicates that the B–H bond is significantly polarized vs the N–H bond. This causes easier hydride donation than proton donation. Furthermore, the highest occupied molecular orbital–lowest unoccupied molecular orbital (HOMO–LUMO) energy gap between N–H σ* and the lone pair of O(CO2) is 1 eV higher than the HOMO–LUMO gap between the B–H σ and π* of CO2. Proton and hydride migrations are the results of the mentioned HOMO–LUMO interactions. The higher barrier of proton transfer to O(CO2) is due to the less favorable proton release from the N–H bond than hydride release from the B–H bond. The Mulliken charges and frontier molecular orbitals (FMOs) are depicted in Figures S1 and S2 in the SI, respectively.

Cis Conformation—Concerted Hydride and Proton Transfer to CO2

We begin with an optimized cis conformation, in which the H–N–B–H dihedral angle is 25° and the H···H distance is 2.33 Å. To prohibit free rotation around the HNBH angle during the simulation and the formation of the trans conformation, we put a potential energy wall at 2.5 Å for the H···H distance. To obtain the FES, we have run 20 ps metadynamics simulations with two collective variables that are the H+···O(CO2) and H–···C(CO2) distances, which both are close to the CO2 molecule in the initial structure. Figure shows the variation of the H+···O(CO2) and H–···C(CO2) distances during the simulation time. As can be seen in Figure , within the first 7 ps of the simulation run, there are large amplitude fluctuations in both the H+···O(CO2) and H–···C(CO2) distances, which indicate a strong interaction between the hydride with C(CO2) and the proton with O(CO2). At 7.1 ps, where H–···C(CO2) and H+···O(CO2) are 2.4 and 3.1 Å, respectively, the hydride and proton transfer to CO2 and HCOOH forms. This points out that a decrease of H–···C(CO2) triggers H+···O(CO2) shortening. On the other hand, in the reverse reaction, which is the formation of CO2, the proton first leaves the O(CO2) and then the hydride transfers to the boron atom, Figure .

Figure 10

Variation of H+···O(CO2) and H–···C(CO2) distances (CVs) versus simulation time.

Figure shows the free energy surface of the concerted hydrogenation of CO2 starting from the cis conformation around the HNBH dihedral. As can be seen at the reactant-state minimum, the reaction proceeds by decreasing the H+···O(CO2) distance at almost constant H–···C(CO2), ca. 3.0 Å. At the H+···O(CO2) distance of 1.8 Å, the H–···C(CO2) distance starts decreasing and after passing through the transition-state region the HCOOH product forms. The reaction is endergonic with a ΔG of 6.9 kcal/mol, and the barrier for concerted CO2 hydrogenation is 21.4 kcal/mol, which is slightly higher than that of the stepwise mechanism starting from the trans conformation, ca. 0.5 kcal/mol.

Figure 11

FES of the concerted mechanism of the hydrogenation of CO2, starting from the cis conformation around HNBH. The energy values (numbers) are in kcal/mol. The contour lines are spaced at 1.5 kcal/mol intervals. For a clearer picture, we have removed the UiO-66 primitive in the figure.

The shape of the FES indicates that the reactant state is a flat surface for H+···O(CO2) within 3.0 to 1.8 Å with decreasing the H+···O(CO2) distance; then, the H–···C(CO2) distance decreases to form the transition state. This behavior illustrates that the concerted mechanism of CO2 is asynchronous, i.e., first H+···O(CO2) decreases and then H–···C(CO2) starts decreasing.

As shown in Figure , in the concerted path, many strong interactions between H+ and O(CO2) are involved for the elongation of the N–H bond before H+ completely leaves. The FES for the concerted mechanism shows these interactions up to O···H distance of 1.5 Å within an energetically flat region. On the other hand, as shown in Figure , lower panel, in the stepwise mechanism, the O···H distance continuously decreases within the first 2.5 ps and the O–H bond forms 400 fs after hydride transfer. Hence, an earlier hydride transfer catalyzes proton transfer.

The preferred conformation of the HCOOH is cis for both pathways. The analysis of the trajectories of the stepwise and concerted paths indicates for the stepwise path in the interval of 2.0–14.0 ps and for the concerted path within 7.0–9.5 ps, at which the HCOOH forms, the HOCH dihedral varies within ±0.78 rad (±45°). The preferred cis conformation of HCOOH has been also indicated through an equilibrium molecular dynamics simulation of the HCOOH inside the cavity of MOF (the variations of the HOCH dihedral are presented in Figure ). Figure S3 in the SI shows the variation of the HOCH dihedral in the HCOOH product molecule during 20 ps metadynamics simulations of both paths.

Figure 12

Evolution of H···H, HCOOH···O(MOF), and Zr···O(C=O of HCOOH) distances during 20 ps equilibrium dynamics simulations. The variation of the HOCH dihedral of the HCOOH (ca. ± 0.5 rad) indicates that the stable conformation of the product during 20 ps simulations is cis.

Dynamics of UiO-66 Alongside the Entire Reaction Mechanism

A computationally proposed and experimentally examined idea of employing MOFs as a solid platform to fix the Lewis pair center involves the retention of the catalytic functionality of the LP inside the cavity of MOF. Hence, MOF may have some effects in stabilizing the reactants and products but does not provide any catalytic action in the reaction. According to our metadymanics simulations and analysis of the trajectories, the oxygen and metal sites of the UiO-66 MOF do not participate in the reaction mechanism. However, to investigate the role of UiO-66 in the stabilization of HCOOH (product of the reaction between H2 and CO2), or the possibility of the interaction between LP centers and HCOOH, we have run 20 ps equilibrium MD simulations on the free HCOOH inside the cavity of the UiO-66 primitive.

Figure shows the evolution of the H···H, HCOOH···O(MOF), and Zr···O(C=O of HCOOH) distances during 20 ps equilibrium dynamics simulations. This equilibrium run indicates an interaction between the carbonyl oxygen of HCOOH and the Zr site of UiO-66. Furthermore, indications of H-bond interaction between HCOOH and O site of UiO-66 are observed, i.e., HCOOH···O(UiO-66) fluctuates around 2.0 Å after 7 ps of simulations. There is no evidence of interaction between the transferred hydride and proton with LP sites. It seems that the thermodynamic force for the HCOOH formation is the reasonably strong interaction between the product molecule (HCOOH) and the MOF, which results in the preferential adsorption of the HCOOH molecules and prohibits the reverse reaction. Hence, HCOOH is continuously removed from the process, rendering high overall conversion rates.

The variation of the HOCH dihedral of the HCOOH (ca. ±0.5 rad) indicates that the stable conformation of the product during 20 ps simulations is cis.

Conclusions

In the present study, the reaction between H2 and CO2 to form HCOOH catalyzed by Lewis pairs, B and N as the Lewis acid and base centers, respectively, covalently bonded to the benzene dicarboxylate (BDC) linker of the UiO-66 platform was investigated using DFT-based metadynamics simulations for the first time. The Lewis pair retains its chemical activity when bonded in the pore, being able to facilely bind H2 and CO2. We simulated free energy surfaces (FESs) for the entire reaction, and it resulted that the reaction mechanism is more eventful than the previously proposed by static-DFT calculations according to the literature. Specifically, the possibility of the trans conformation around the HNBH dihedral of the dissociated H2 molecule is observed with the present metadynamics simulations, which leads to an alternative stepwise mechanism for the hydrogenation of CO2. The stepwise mechanism starts with hydride transfer to the carbon of CO2 and is followed by proton transfer to the oxygen of CO2. Furthermore, according to our FESs simulated by metadynamics, the CO2 desorption from the Lewis pair centers has a lower barrier than H2 formation from the dissociated H···H on the Lewis pair. Metadynamics simulations reveal the possibility of two alternative mechanisms for the hydrogenation of CO2. A comparison of the energetics of the stepwise and concerted mechanisms indicates that the stepwise mechanism has a slightly lower barrier than concerted, which means that the two pathways are kinetically similar. The CO2 hydrogenation is endergonic in both mechanisms. In the stepwise mechanism, hydride transfer takes place earlier than proton transfer and catalyzes the ongoing proton transfer. The analysis of the structure of the TSs shows that the C···H– distance is shorter in the trans conformation than that in cis at the TS. On the other hand, the O···H+ distance is longer in the trans conformation than that in cis. Moreover, according to the simulated FESs for the stepwise and concerted mechanisms, it turned out that in the concerted mechanism the reactant state is constructed on a flat FES. The analysis of the trajectories at each step of the CO2 hydrogenation indicates that each point in the CV space is explored several times during the simulation (Figures , 6, 8, and 10). This leads to improved accuracy in the free energy estimate. Conceptually, the required protic and hydridic hydrogens for CO2 hydrogenation can be provided by Lewis bases and acids, respectively, and as a result, the CO2 activation and subsequent hydrogenation are a metal-free route for the production of hydrogen-rich C1 fuels.