Free energy calculations with enhanced sampling simulations are less computationally demanding than using unbiased AIMD but can require a large number of DFT calculations to converge sampling statistics and quantify statistical uncertainties. Therefore, we adopt a computationally efficient DFT model, consisting of a 2 × 2 unit cell of an N 2 molecule on a twolayer Ru(0001) surface (Figure ) (see Computational Methods). In CFF-AIMD simulations, we use two-dimensional (2D) collective variables (CVs) to describe the bond breaking/formation and adsorption/desorption processes for the dissociation of N 2 on Ru: the N-N distance (𝑟) and the molecule-surface distance (𝑍) (Figure ). An ensemble of final dissociated states (2N*) can be sampled because there are two stable N* adsorption sites: hexagonal close-packed (HCP) hollow and face-centered cubic (FCC) hollow sites, which differ in the presence of a ruthenium atom in the layer below the HCP hollow site (Figure ). In total, we performed ~500 ps of enhanced sampling simulations using multiple walkers or replicas, amounting to ~17 million timesteps (see Table for a summary of simulations performed in this work). As demonstrated in Figures , CFF-AIMD successfully samples the configurational space describing the nitrogen dissociation reaction. Using 16 walkers at 1400 K, CFF-AIMD identifies the initial (N 2 *) and final states (2N*) of the surface reaction, respectively. By applying the nudged elastic band (NEB) method to the BRANN-fitted FES (Supporting Information), we find the minimum free energy path (MFEP) between the adsorbed molecular state (N 2 *) minimum at (𝑟, 𝑍) = (1.1, 2.7) Å, where the molecular axis is nearly perpendicular to the surface (the on-top configuration), and the dissociated atomic state (2N*) minimum at (𝑟, 𝑍) = (2.7, 1.1) Å, where nitrogen atoms adsorb to predominantly HCP hollow sites. This MFEP has a sharp curvature with a narrow passage connecting N 2 * and the transition state (TS) (Figure ). To examine how neural networks accelerate free energy calculations for surface reactions, we also perform simulations using an enhanced sampling method without neural networks: ABF-AIMD. ABF is chosen to mimic the adaptive sampling strategy in CFF without the improvements brought forth by neural networks (see Supporting Information on the CFF method description). CFF-AIMD converges to the underlying FES faster than ABF-AIMD (Figures ). ABF-AIMD does not converge after extensive simulations (~102 ps) (see Figure ). Regions proximal to the N 2 * state, which contains sharp free energy differences (~2 eV) over a narrow range of configurational space (𝑟 = 1.1-1.5 Å , 𝑍 = 2-3 Å), are undersampled (see Figure ). Specifically, ABF relies on determining the force a molecular system experiences along an order parameter and adds an opposing force (the bias) in an iterative fashion, such that the total force is zero. This iterative process is particularly slow across large free energy barriers (e.g., the configuration space proximal to N 2 *), where sampling is too sparse to obtain accurate estimates of the forces that arise from the underlying FES (see the Supporting Information on ABF-AIMD). In contrast, our CFF-AIMD method directly approximates the global FES using neural networks based on the potential energy and forces computed from on-the-fly DFT calculations. CFF-AIMD arrives at more accurate estimates of the biasing forces required by interpolating the free energies between both sparsely-sampled and well-sampled configurations at each iteration (see also the Supporting Information Sec. 5.1). These combined effects lead to a faster convergence as well as a more uniform sampling with CFF-AIMD compared to ABF-AIMD (Figures ), thereby highlighting the ability of the CFF-AIMD method to efficiently and accurately sample free energy landscapes containing narrow valleys and sharp curvatures, which are common in molecular dissociation on surfaces .