Combined DFT and Kinetic Monte Carlo Study of UiO-66 Catalysts for γ-Valerolactone Production

Zr-based metal–organic frameworks (MOFs) are excellent heterogeneous porous catalysts due to their thermal stability. Their tunability via node and linker modifications makes them amenable for theoretical studies on catalyst design. However, detailed benchmarks on MOF-based reaction mechanisms combined with kinetics analysis are still scarce. Thus, we here evaluate different computational models and density functional theory (DFT) methods followed by kinetic Monte Carlo studies for a case reaction relevant in biomass upgrading, i.e., the conversion of methyl levulinate to γ-valerolactone catalyzed by UiO-66. We show the impact of cluster versus periodic models, the importance of the DF of choice, and the direct comparison to experimental data via simulated kinetics data. Overall, we found that Perdew–Burke–Ernzerhof (PBE), a widely employed method in plane-wave periodic calculations, greatly overestimates reaction rates, while M06 with cluster models better fits the available experimental data and is recommended whenever possible.


S1. Analytical Expression for the Reaction Time
The reaction mechanism under investigation can be conveniently expressed as: in order to estimate the reaction time (see main manuscript).In it, R and P represent the reactant (methyl levulinate) and the product (γ-valerolactone), respectively, whereas C is the catalyst and I is a reaction intermedium.Notice that where x is the concentration of a given species, indicated by a subscript, and the 0 superscript indicates the initial concentration.
From the previous mechanism, we can write: If we assume the steady state approximation, i.e. the intermediate is consumed as quickly as it is generated, in combination with eq.(S2), we get: Thus, integrating eq.(S3): with t α being the time at which x R = (1 − α) • x 0 R , lead us to: where α is the conversion.In eq (S7), φ 0 represents the ratio between the initial populations of the reactant and the catalyst:

S2. Calculation of thermal rate constants
To simulate the UiO-66 MOF's rigidity, the four para-carbon atoms of the benzoate groups were kept fixed during the calculations.As a result, the systems with frozen atoms (I0 to I10 and the transition states) will contribute with a total of 3 • (N − 4) vibrational modes to the partition function, where N represents the total number of atoms.Notably, the overall translations and rotations are excluded in the calculation of the partition function.It is essential to mention that the Pilgrim software used for these calculations was not originally designed for systems with frozen atoms.To overcome this limitation, a modified version of Pilgrim was employed for the calculation of the rate constants.
Rate constants were calculated using several DFT (Density Functional Theory) functionals.Hessian matrices were calculated using the Stuttgart effective core potentials (SDD) in combination with the def2-TZVP basis set for the Zr atoms and def2-SVP for O, C, and H atoms.A cutoff of 50 cm −1 was implemented for the resulting vibrational frequencies.The electronic energies were refined using single-point calculations using the def2-TZVP/SDD basis set for Zr and def2-TZVP for the O, C, and H atoms.
The calculated rate constants and ratios between the initial concentration of ML and the catalyst, φ 0 , can be found in Tables S1 and S2.We highlight that: • Multi-structural transition state theory (MS-TST) was used in the calculation of the rate constant of R HAT , in order to account for the conformational wealth of ML (see section S6): In the previous equation, k SS−TST is the single-structure TST rate constant and F MS,X [X = ‡ (the TS), or R (reactants)] accounts for the contribution of the conformations of ML, that is, it is the ratio between the MS and SS harmonicoscillator partition functions.The rate constant k SS−TST was calculated for each DFT functional, but F MS was exclusively calculated for M06 and assumed to be the same for M06L anf PBE-D3.The value of F MS at T = 413 K is 0.46.
• Forward and backward rate constants for R nuc were calculated using SS-TST.
• The rate constant for R elim (and consequently for R n+e ) was calculated using canonical variational transition state theory (CVT): The variational coefficient, Γ CVT , associated with TS6-7 is constant between 403.15 and 413.15K and its value is 0.87, 0.38 and 0.70 for M06, M06L and PBE-D3, respectively.Figure S1 shows the electronic energy profiles for cluster and periodic models computed at PBE level with Grimme dispersion D2 (while in Figure 3, cluster PBE-D3 was compared with periodic PBE-D2).The dispersion correction does not affect much the activation energies of the reaction.

S3.2. Cluster PBE-D3 versus cluster PBE-D3 with a diffuse function basis set for O atoms
To better represent electronegative atoms with the basis set, we employed diffuse functions for all oxygen (O) atoms and referred to this modified basis set as "basis set diffuse".The results, shown in Figure S2, indicate that the inclusion of diffuse functions does not significantly impact the activation energies of the reaction.

S6.1. Reference Z-matrix
Here, we show the Z-matrix used in the conformational search with TorsiFlex.The value associated to each internal coordinate is listed below.Bond distances are given in Å; bond and dihedral angles, in degrees.# t e m p e r a t u r e s (K) f o r p a r t .f u n c t i o n s f r e q s c a l L L 1.000 # f r e q .s c a l i n g f a c t o r ( LL ) f r e q s c a l H L 1.000 # f r e q .s c a l i n g f a c t o r (HL) sigmamj 0.02 # max v a l u e f o r sigma ( Mj ) ; >= 0.01

S6.4. Conformers
Relative energy for the conformers for methyl levulinate, calculated at M06/Def2SVP, is shown in Table S4.The two reference energies, one for the total electronic energy (U ) and another for the total electronic energy plus the vibrational zero-point energy (U +ZPE), are:

Figure S2 :
Figure S2: Electronic energy profiles (in kcal•mol −1 ) of defective UiO-66 using cluster at the PBE-D3 level and at the PBE-D3 level with a diffuse function basis set for O atoms.R = (CH 2 ) 2 CO 2 Me.

Figure S3 :
Figure S3: Relative Gibbs energies calculated at the PBE-D3 level for the regeneration routes at 403.15 K (in kcal•mol −1 ).
3 # c o n t r o l s t h e c o n n e c t i v i t y c r i t e r i u m #−−−−−−−−−−−−−−−−−−−−−−−# # S t o r a g e # #−−−−−−−−−−−−−−−−−−−−−−−# d i r l l f i l e s L L / # f o l d e r t o s t o r e LL conformers d i r h l f i l e s M 0 6 / # f o l d e r t o s t o r e HL conformers t m p l l / s c r a t c h / LL ml / # f o l d e r f o r LL temporal f i l e s tmphl / s c r a t c h /M06 ml/ # f o l d e r f o r HL temporal f i l e s #−− Guess geom ( Conn , S i m i l , Hard , S o f t ) t e s t s O 1 1 1 1 # f o r Opt geom ( Conn , Redun , Hard , S o f t ) d i s t 1 D 7 # domain s i z e about each p o i n t ( d e g r e e s ) epsdeg 2 # max d i f f between two i d e n t i c a l a n g l e s ( d e g r e e s ) #−−−−−−−−−−−−−−−−−−−−−−−# # Gaussian c a l c u l a t i o n s # #−−−−−−−−−−−−−−−−−−−−−−−# optmode 1 # 0 : opt ( z−m a t r i x ) , 1 : opt ( modredundant ) f c c a r d s y e s # Use LL H e s s i a n i n HL opt ( y e s /no ) n p r o c l l 8 # Number o f t h r e a d s ( low−l e v e l ) n p r o c h l 8 # Number o f t h r e a d s ( high−l e v e l ) memll 16GB # dynamic memory ( low−l e v e l ) memhl 16GB # dynamic memory ( high−l e v e l ) l o w l e v e l HF 3−21G # low−l e v e l o f c a l c u l a t i o n h i g h l e v e l M06 Def2SVP i n t=u l t r a f i n e # high−l e v e l o f c a l c u l a t i o n #−−−−−−−−−−−−−−−−−−−−−−−# # P a r t i t i o n f u n c t i o n s # #−−−−−−−−−−−−−−−−−−−−−−−# temps 403.15 # t e m p e r a t u r e s (K) f o r p a r t .f u n c t i o n s temps 413.15

Table S1 :
Rate constants calculated by different electronic structure methods for the reactions indicated in the "Reaction kinetics simulations and comparison with experiments" section.Rate constants are in cm 3 • molecule −1 • s −1 for R HAT and in s −1 for the remaining reactions.Forward and backward rate constants are indicated with (fw) and (bw), respectively.

Table S2 :
Experimental (Exp)and calculated values by different electronic structure methods for the ratio between the initial concentrations of ML and the catalyst (φ 0 ).The value of the product conversion (α) for each experiment is also shown.

Table S3 :
Computed geometries of intermediates and transition states across different DFs.Selected distances in Å.

Table S4 :
: -459.57003HartreeNotice that each conformer in TableS4presents a torsional enantiomer.Therefore, the total number of located conformers is twice the number shown in the table.Relative total energy, U , and relative total energy plus zero-point energy, U + ZPE, for the located conformers of methyl methyl levulinate (kcal/mol).Dihedral angles are in degrees.In this section, the Cartesian coordinates (in Å) of each conformer optimized at the M06/Def2SVP level are listed. ref