Estimating Free Energy Barriers for Heterogeneous Catalytic Reactions with Machine Learning Potentials and Umbrella Integration

Predicting the rate constants of elementary reaction steps is key for the computational modeling of catalytic processes. Within transition state theory (TST), this requires an accurate estimation of the corresponding free energy barriers. While sophisticated methods for estimating free energy differences exist, these typically require extensive (biased) molecular dynamics simulations that are computationally prohibitive with the first-principles electronic structure methods that are typically used in catalysis research. In this contribution, we show that machine-learning (ML) interatomic potentials can be trained in an automated iterative workflow to perform such free energy calculations at a much reduced computational cost as compared to a direct density functional theory (DFT) based evaluation. For the decomposition of CHO on Rh(111), we find that thermal effects are substantial and lead to a decrease in the free energy barrier, which can be vanishingly small, depending on the DFT functional used. This is in stark contrast to previously reported estimates based on a harmonic TST approximation, which predicted an increase in the barrier at elevated temperatures. Since CHO is the reactant of the putative rate limiting reaction step in syngas conversion on Rh(111) and essential for the selectivity toward oxygenates containing multiple carbon atoms (C2+ oxygenates), our results call into question the reported mechanism established by microkinetic models.


Hyperparameters
Within the Gaussian approximation potential (GAP) framework, the total energy (E tot ) of an atomic configuration with positions X n is expressed in terms of local energy contributions fitted by a Gaussian Process.We use here two-body (2b) and many-body (mb) contributions, which are defined as sums of kernel functions (k(i, j)) In this context, i and j are atom indices, r ij ≡ |r i −r j | represents an interatomic distance between atom i and j within a cutoff radius (r cut,2b ) and p i is the smooth overlap of atomic position (SOAP) vector of atom i, which describes its atomic environment within the cutoff sphere r cut,SOAP .A cutoff transition width r ∆ defines a region in which the neighborhood density of the SOAP representation is smoothly damped to zero.δ 2b and δ mb define weighting factors of the different contributions and represent the standard deviations of the respective Gaussian Processes.c 2b,n and c mb,n are the fitting coefficients obtained by minimizing the loss function in the training, where y l and ỹl are the DFT calculated and predicted quantities (here energies or forces).M is the training set size and R a regularization term weighted by σ 2 .These σ 2 values (σ 2 E for energies and σ 2 F for forces) represent the assumed variance of the errors and can be set globally for the entire training set, per configuration or in case of fitting forces per atom.To accelerate the training procedure, only a set of N sparse representative training points for the 2-body and many-body terms are used, respectively.These points are selected uniformly for the 2-body and based on the CUR algorithm for the many-body contributions .
In our iterative training scheme (see Fig. S5), we heuristically adjust σ E over the iterations so that the potential is more stronglz regularized during earlier iterations where the training set is far from complete.To this end, we start with a value of 0.01 eV and adjust σ E to √ RMSE (for the validation set) if that value is smaller than 0.01 eV (upper threshold) or larger than 0.001 eV (lower threshold).To illustrate this approach, the different validation set √ RMSE of the energies and σ E are both plotted as a function of iterations in Fig. S1.We also use different σ F for each atom in the training set.These regularization parameters are adjusted to the DFT calculated forces according to the heuristic formula Here, the idea is to incorporate both the DFT force norm of each atom (F DFT norm ) and the average force norm in each configuration (F with w i = 1 σF,i .In addition, we only incorporate forces of the surface atoms into the training evaluation (the CHO molecule plus the upper layer of Rh) and mask the remaining Rh-atoms (three bottom layers).This helps to decrease the force MAEs of carbon, oxygen and hydrogen as the surface slab is dominated by Rh-atoms (Rh : C + H + O = 36 : 1 + 1 + 1). Figure S3 indicates the surface and masked atoms.To describe the short and medium range contributions in an atomic environment differently, we use a multi-SOAP approach for the light species (C, H, O).This means that two SOAP representations with different length scale parameters are applied to obtain a comprehensive description of the atoms that are directly involved in the reaction.In contrast, the Rh-atoms are described by a single SOAP lengthscale.In the following, we will distinguish the representations for the light elements by using the notation SOAP 1 as well as SOAP 2 and SOAP Rh for the Rh-atoms.
We use a Gaussian kernel function to express the similarity between the 2-body contributions and a polynomial kernel for the many-body contributions.To avoid unphysical atomic clashes in early training iterations (where the potential is not necessarily fully robust), we add simple diatomic baseline potentials for all element combinations, as reported in [1].The corresponding O − O potential is plotted in Fig. S4 as an example.These baseline potentials account for the repulsive nature at small interatomic distances and keep structures away from nonphysical configurations in dynamical simulations.If no baseline potential is used, configurations with small distances have to be provided in the training set.

Iterative Training
We start with an initial training set of 50 configurations including 10 structures of the stoichiometry Rh 36 CHO (namely the NEB trajectory in Fig. 2 main text) and 12 structures of stoichiometry Rh 36 with an optimized lattice constant of 3.85 Å.In addition, we add some configurations which are computational efficient to evaluate, namely the four isolated atoms (Rh, C, O, H) and 24 dimer combinations of the three light species in gas-phase.With this initial training set, we train a first GAP model (iteration 0).At this stage, geometry optimizations on biased potential energy surfaces in ten randomly selected umbrellas were performed.The energies and forces of these 10 generated structures were recalculated with DFT and the corresponding configurations added to the training set.In subsequent iterations, 2.5 ps of biased dynamics at 573 K (umbrella sampling with k = 15) were carried out, again FIG.S5: Iterative training workflow for fitting GAP models.Our approach is explained in the main text of the SI. in 10 randomly selected windows.From each of these trajectories, one configuration was selected by farthest point sampling (FPS) [2] and added to the training set.After 12 iterations, we additionally perform NEB calculations and further add these images (8 configurations between initial and final state) to the training set.Note, that we rejected structures with force components higher than 15 eV/ Å .The workflow is illustrated in S5.We stop our iterations after 17 iterations and validate the GAP model by monitoring the energy and force errors over the different iterations (see Fig. 3 main text) as well as by comparing NEB trajectories produced by the GAP model and DFT (see Fig. S6).The final training set contains 292 configurations and the potentials were trained on DFT forces as well as atomization energies (AE) where E slab represents the potential energy of the surface slab and E atom the energy of the isolated atom s.DFT@GAP GAP@GAP DFT@DFT GAP@DFT FIG.S6: Model validation by comparing NEB trajectories produced with the GAP model (GAP@GAP) and DFT (DFT@DFT).In addition, the energies of the NEB images generated with the GAP model have been recalculated with DFT (DFT@GAP) and vice versa (GAP@DFT).The BEEF-vdW model was trained on a subset of the configurations used for the revPBE models, where energies and forces where recomputed with the BEEF-vdW functional.Here, for technical reasons only the configurations inculding the full surface adsorbate system (with composition Rh 36 CHO) where used.This is because isolated atoms and dimer configurations are more challenging to converge with the plane-wave QE code.Since these configurations are mainly included to provide robust prior information about interatomic interactions in early iterations, they can be neglected when retraining on the full revPBE training set.For similar reasons, formation energies were used as the fitting target, defined as: with Identical hyperparameters as in the last iteration of the revPBE potential are used, merely adjusting σ F according to the BEEF-vdW forces.
FIG. S1: Energy regularization parameter (σ E ) and the square root RMSE of atomization energies per atom of the validation set as a function of iterations.Dotted lines represent the upper and lower threshold of σ E .
DFT norm ) into the regularization parameter.The parameters σ min , C and A are additional hyperparameters, which are set to σ min = C = √ σ E and A = 0.01, while different values are used for dimer configurations (σ min = C = 0.1).The force regularization parameters are plotted as a function of the force norms in Fig. S2.Note that the force MAE stated in the main text represent MAEs weighted by the individual σ F :

F
FIG. S2: Force regularization parameter (σ F ) as a function of the force norm on an atom (F DF T norm ) and the average force norm on a configuration (F DFT norm ) .
FIG. S3: Illustration of surface and masked atoms during training.
FIG. S4: Example O − O baseline potential as a function of the interatomic distance from [1].

TABLE S2 :
Vibrational frequencies of initial and transition states at the revPBE+vdW surf level in cm −1 .* margraf@fhi.mpg.de

TABLE S3 :
Vibrational frequencies of initial and transition states at the BEEF-vdW level in cm −1 .

TABLE S4 :
Energy barriers (in eV) computed with different codes and functionals using the BEEF-vdW and revPBE+vdW surf initial and transition state geometries.FIG.S7: Estimation of CHO half-lives from unbiased MD simulations.Statistical errors are estimated via bootstrapping.