For R2, the US/MLIP-MD free energy profile was computed using r 1 -r 2 as the reaction coordinate, split into 30 equally-spaced windows with a force constant of 20 eV Å -2 . The simulation ran for 40 ps at 365.6 K in every window and the first 10 ps in each window were used for equilibration and discarded from the analysis. Free energy was calculated by the weighted histogram analysis method (WHAM) as implemented in mlp-train. WTMetaD was performed along the same reaction coordinate, r 1 -r 2 , as the US. The width of the deposited Gaussians was 0.07 Å with the initial height 1 2 k B T = 0.363 kcal mol -1 at 365.6 K to be consistent with the experimental condition and the bias factor of 50. Gaussians were deposited each 100 fs. WTMetaD was run in 10 independent replicas for 500 ps each, starting from the same configuration but with different velocities initialised using the Maxwell-Boltzmann distribution. WTMetaD+IB was performed using the inherited bias from the 16 th AL iteration as an initial bias. The parameters of the deposited Gaussians were the same, with the bias factor increased to 80. Ten independent WTMetaD+IB runs were performed with a simulation time of 250 ps each. The free energy profile was reconstructed by reweighting based on the final bias potential, 90 using kernel density estimation as implemented in PLUMED with a bandwidth of 0.02 to compute the histograms. For R3, the free energy surface was investigated using WTMetaD/MLIP-MD. The simulation system comprised the substrate in 86 DCM molecules within a box with the size of 22 Å using periodic boundary conditions (PBC) to maintain the density of DCM (1.33 g/cm 3 ).