Modeling the potential energy surface by force fields for heterogeneous classification, applications, and challenges

Authors: Chenglong Qiu, Tore Brinck, Jiacheng Wang

Abstract

The concept of the potential energy surface (PES) in computational simulations is essential for studying material properties and heterogeneous catalytic processes. However, constructing the PES using quantum mechanical methods is computationally expensive and typically limited to small systems. Force field methods, which rely on quantum mechanical data, use simple functional relationships to establish a mapping between system energy and atomic positions or charges. Force field methods are more efficient for handling large-scale systems, such as catalyst structures, adsorption and diffusion of reaction molecules, and heterogeneous catalytic processes. To further promote in-depth research in this field, this review introduces the classification, development, and characteristics of various force field methods including classical force fields, reactive force fields, and machine learning force fields. It summarizes the forms, fitting methods, and distinct periods of these force field methods. Additionally, these force field approaches are compared in terms of their applicability, accuracy, efficiency, and fitting methods. Finally, the optimization and challenges of force field methods in constructing the PES are discussed. It is expected that this review will assist researchers in selecting and applying different force field methods more effectively to promote in-depth understanding of catalytic reaction mechanisms and the efficient design of catalysts.

Introduction

Computational simulation has emerged as a powerful tool for investigating chemical reaction processes and the physical and chemical properties of materials.^1–6^ The potential energy surface (PES), based on the assumption that electron and nuclear motions can be separated, is a crucial concept in computational simulations, representing the total energy of a system in a certain state.^7,8^ The PES is widely applied in fields such as physics and chemistry, particularly in their theoretical subfields. From a geometric perspective, the energy landscape is the plot of the energy function over the configuration space of the system. It is used to explore the properties of atomic structures, such as determining the minimum energy configuration of a molecule or calculating reaction rates.^9,10^ Additionally, in dynamic simulations based on Newton's laws, the force Fi exerted on each atom must be known at each time step for numerical integration of the equations of motion.^11^ This force can be derived from the PES by using the relation Fi = ∂E/∂ri, where the force is the negative gradient of the potential energy E with respect to the atomic position ri. Forces are also used in geometric optimization to identify the special structure of the system that corresponds to the critical point on the PES.^12,13^ For instance, a saddle point represents a transition state—the peak energy point along the reaction coordinate that determines the most energetically favorable path between reactants and products. The magnitude of the reaction energy barrier can be calculated as the energy difference between the saddle point and the two energy minima it connects. Therefore, the PES is an essential tool for analyzing reaction processes and predicting system evolution. The primary challenge lies in constructing the PES both efficiently and accurately.

Quantum mechanics (QM) and the force field method are the primary methods for constructing PESs (Fig. 1a). For smaller systems, PESs constructed by QM can accurately describe molecular properties, crystal structures and microscopic reactions. In recent years, QM-based simulation methods have gained increasing popularity, partly due to the development of some software packages that facilitate the generation of PESs. For example, periodic density functional theory (DFT) codes such as VASP and CP2K are flexible for extended systems, albeit requiring external automation or scripting for reaction-path sampling.^14,15^ And codes such as Q-Chem and Gaussian offer built-in tools for PES scans along reaction coordinates.^16,17^ This increased availability has proven particularly useful in material design, where QM often serves as a theoretical guide and screening tool.^18,19^ However, the computational cost inherent to QM-level calculations severely limits the scalability of simulations.^20^ In QM, the electronic structure and energy of a system are determined by solving the Schrödinger equation (Fig. 1b), with the analytical solution applicable only to two-body systems, such as hydrogen atoms. For multi-atom systems, several methods (e.g., semiempirical wavefunction,^21^ density functional theory,^22^ and CCSD(T)^23^) have been developed to approximately solve the Schrödinger equation. Despite these approximations, obtaining a precise numerical solution remains a computationally demanding task. For instance, CCSD(T), a high-precision ab initio method for electron correlation, scales ∝ N,7 with N being the number of atoms.^24^ Consequently, constructing a PES using QM to model the dynamic evolution of large chemical systems containing diverse molecules is impractical.

In contrast to QM, the force field method uses a simple functional relationship to establish the mapping between the system's energy and the atomic positions and charges (Fig. 1b). Compared to solving the Schrödinger equation, calculating system energy using the force field method is significantly less complex, allowing it to handle large-scale systems (e.g., polymers, biomolecules, and heterogeneous systems) more efficiently. The force field method dates back to 1924,^25^ when Jones proposed a molecular model involving a repulsive force λnr^−n^ and an attractive force λ3r^−3^. Building on this, the Lennard-Jones potential function was developed to describe interactions between non-bonded atoms or molecules.^26,27^ With the advancement of computational simulations, increasingly accurate forms of potential functions have been introduced. Based on their forms and the types of systems they apply to, current force fields are categorized into three classical force fields, reactive force fields, and machine learning (ML) force fields. The construction of a force field-based PES primarily relies on energy values calculated from discrete geometric configurations by QM, followed by fitting a PES using these discrete data points. Thus, the accuracy of the force field method is influenced by the quality of the QM calculations. Additionally, due to errors inherent in the fitting process, the force field method cannot achieve the precision of QM (Fig. 1c). Consequently, force field methods often trade computational cost for accuracy, enabling simulations at scales that are orders of magnitude beyond the reach of QM.^28^

Various types of force fields and their applications have been extensively reviewed in the literature. For examples, Wang et al. summarized the application of different force field methods in mechanism exploration and performance prediction of electrocatalysis.^29^ Thomas and Han et al. provided an overview of the development and application of the ReaxFF reactive force field.^30,31^ Oliver et al. presented a detailed mathematical and conceptual framework of ML force fields, along with their applications and the chemical insights they offer.^28^ Additionally, Cheng et al. elaborated on the principles and application of ML force fields, in conjunction with global optimization algorithms, to identify in situ active sites in heterogeneous catalysis.^32^ However, most reviews focus on the application of a single force field,^28,30,33–38^ and lack the comparison of different force fields, particularly in terms of fitting methods and their applicability to various systems.^39–41^ Furthermore, the principle and development of the PES constructed using force field methods are also inadequately summarized. Therefore, a comprehensive review is necessary. Summarizing the development trends of force field methods and the characteristics of different force fields can not only assist researchers in selecting the most suitable force fields for specific applications and provide valuable theoretical guidance for experimental work, but also foster innovation in new methods to accelerate progress in this field.

To this end, this review systematically summarizes the development and application of these three force field methods (classical force fields, reactive force fields and ML force fields), and compares them in terms of applicable systems, computational accuracy and cost, fitting methods, and other relevant aspects. For each force field, we discuss its general form in detail and outline the corresponding parameterization strategies. Examples of various applications are provided to illustrate the range of systems that can be modelled with the respective force field. Finally, we discuss future directions for improvements and potential challenges in constructing PESs using force field methods.

Classification and comparison of force fields (classical force fields, reactive force fields, and machine learning force fields)

Introduction to various force fields

According to the Born–Oppenheimer approximation, the energy of a molecule can be expressed as a function of the spatial coordinates of the individual atoms, each characterized by distinct atomic properties or parameters.^7,42^ Consequently, a molecular force field consists of two a functional form that describes the interactions between atoms and the force field parameters specific to each atom. Based on their functional forms and applicable systems, current force fields are categorized into three classical force fields, reactive force fields and ML force fields.

The classical force fields calculate a system's energy using simplified interatomic potential functions. This approximation is well-suited for modeling nonreactive interactions, such as bond stretching and angle bending (represented by harmonic functions), dispersion force (represented by the Lennard-Jones potential), and electrostatic interactions (represented by atomic charges). Currently, a variety of classical force fields with different simplified formulae have been developed for different types of molecular systems. Classical force fields typically contain between 10 and 100 parameters, which often possess clear physical meanings and are relatively easy to interpret (Table 1). However, due to the simplicity of their formulae, such descriptions are inadequate for modeling changes in atom connectivity such as bond breaking and formation during reactions, which also lead to a reduction in calculation accuracy. Despite this limitation, this method significantly accelerates computations (Fig. 1c). The simulation length scale can reach 10–100 nm for extended systems, with time scales ranging from tens to hundreds of nanoseconds, and occasionally extending to the microsecond range on modern hardware. Therefore, the classical force fields are particularly suitable for describing the motions of atoms or molecules driven by their interactions. Thermodynamic or kinetic properties such as adsorption, diffusion, dissolution, separation and stress–strain can be further studied through statistical analysis of the motion behavior of particles.^43–45^ For example, the calculation of the diffusion coefficient from molecular dynamics (MD) simulations involves simulating particle motion, extracting the time series of particle positions, and calculating the root mean square displacement (MSD) over . The diffusion coefficient is then extracted using Einstein's .

The limitations of the classical force fields in describing the reactive processes motivate the inclusion of connection-dependent terms, resulting in the development of reactive force fields. Bond order is a key concept in reactive force fields, describing the strength and properties of bonds between atoms. The bond order value is calculated using specific equations, and it is dynamically adjusted based on the relative positions of atoms and their local environment. Therefore, reactive force fields can describe the breaking and formation of chemical bonds, as well as the conversion between reactants and products, enabling the modeling of reaction processes on a large scale. Their parameters are derived from a combination of physical principles and empirical insights, often involving some degree of abstraction. The number of parameters typically ranges from 100 to 500. And reactive force fields can incur a relatively high computational cost (approximately 10–100× that of classical force fields). The typical simulation length scale for condensed-phase reactive systems is 5–20 nm, with accessible time scales of 1–10 nanoseconds for large systems. In particular, ReaxFF, the most widely used reactive force field, can simulate reaction events at the interfaces of solid, liquid, and gas phases because the parameters for each element in the force field can be transferred across different phases. For example, in the simulation of an oxygen evolution reaction (OER) catalyzed by a metal oxide, the oxygen atoms use the same parameters, whether they are in the gas phase as O2, in the liquid phase within an H2O molecule, or bound in a solid metal oxide. Additionally, since the catalytic process involves not only the reaction but also the migration of molecules, such transferability allows ReaxFF to simulate the influence of dynamic factors, such as solubility and diffusivity, on the catalytic process. As a result, ReaxFF can simulate systems involving multiphase complex processes.

Classical force fields and reactive force fields rely on predefined mathematical functions to describe atomic interactions, often lacking the flexibility and accuracy to model complex chemical environments. In contrast, ML force fields represent an emerging class of computational models that use ML algorithms to construct a system's PES. These data-driven approaches model the PES directly from data, enabling them to capture intricate interactions and chemical behaviors with high accuracy. The number of parameters typically ranges from 10^3^ to 10^6^ (for the emerging foundational/universal ML force fields, the number of parameters can even exceed 10^6^), and they are represented through mathematical or ML models without fixed physical forms. ML force fields can achieve quantum-level accuracy within the training domain, but their performance is constrained by extrapolation limitations and computational scaling (often 10–100× the cost of reactive force fields). This high precision comes at the cost of greater computational expense compared to the other two force fields. To reduce the fitting cost, researchers apply transfer learning to emerging foundational or universal ML force fields by leveraging pre-trained models, such as those from the Open Catalyst Project (OCP) and Materials Project (MP). The simulation length scale depends strongly on descriptor efficiency, generally reaching 2–10 nm in current implementations. The typical time scale is 0.1–1 nanoseconds for extended systems, although linear-scaling ML force fields may extend this range. Nevertheless, ML force fields are more efficient than QM and are thus applicable to large-scale systems (Fig. 1c). For example, water and aqueous systems often involve complex interactions, such as hydrogen bonding, that classical force fields may not accurately capture. As an efficient alternative to QM, DeepMD employs neural networks to model complex atomic interactions and construct PESs to predict hydrogen bond networks, solvation effects, and the diffusion behavior of water.^46,47^ Due to their accuracy and efficiency, ML algorithms like neural networks, Gaussian processes, and ridge regression are increasingly applied in areas such as catalysis, materials science, and drug design.^48–50^ However, ML force fields also have limitations, including high fitting costs, poor model generalization, and lack of interpretability.

Considering the current state of force field development, we believe that no single method currently bridges all relevant catalytic scales (from atomic events at ∼ fs, Å to mesoscale ones at ∼ ms, μm). Therefore, researchers often adopt hybrid simulation strategies. For example, the QM/MM approach applies ML force fields or reactive force fields to the key reactive region, while describing the remaining parts with classical force fields.^51^ Another strategy is coarse-graining,^41^ in which coarse-grained parameters are derived from atomic-scale ML force fields or reactive force field trajectories, thereby extending simulations to the μs–ms regime or beyond. However, ML force fields augmented with active learning and hybrid coupling show the most promise for unifying length- and timescales in heterogeneous catalysis simulations.^33,52^

Development periods for force fields

After establishing a fundamental understanding of force field concepts and classifications, we will further explore the development of force field methods. This involves reviewing their evolution process, analyzing key advancements and breakthroughs across different periods, and identifying the underlying driving forces behind their continuous development to help us better understand the characteristics and applicability of various current force fields. Lennard-Jones, ReaxFF, and neural networks are selected as representatives of classical force fields, reactive force fields, and ML force fields, respectively. Fig. 2 shows the number of publications in molecular simulations that have utilized these approaches. Based on the trends in the number of publications over time, it can be found that there are three key nodes (2003, 2011 and 2017) in the development of force fields.

Before the 21st century, force fields were typically based on simplified potential energy functions to describe atomic interactions. The classical force fields, represented by the Lennard-Jones potential, were introduced early on and proved suitable for most non-bonding interaction due to their simplicity. Building on this foundation, force fields were further developed to include additional physical details, such as bond stretching and bond angle bending, making them more applicable to molecular systems. Several classical force field models have been developed in this period, such as AMBER,^53^ CHARMM,^54^ OPLS,^55^etc. However, during this early stage, the body of related research was limited. With advancements in computational power and molecular simulation methods in the 21st century, classical force fields began to account for more complex multi-body effects and environmental factors, leading to their broader adoption. These force fields were further refined to include polarization effects between atoms, enhancing both their accuracy and applicability. Additionally, to better describe reactive processes, van Duin introduced the ReaxFF force field, based on the concept of bond order in 2001.^56^ Since then, with ongoing improvements in ReaxFF parameters, its transferability has become increasingly evident, making it widely used in fields such as heterogeneous catalysis, atomic deposition, organic molecule combustion, and liquid-phase chemistry after 2011.

The advent of powerful ML algorithms significantly enhanced the accuracy of PES construction. Although Blank et al. proposed the use of neural networks to represent PESs for molecule-surface scattering in 1995,^57^ the development of ML force fields progressed slowly due to computing power limitations. However, in recent years, with advances in technology, neural network potentials and other ML force fields have seen rapid growth in their application across chemistry, physics, and materials science. The evolution of these force fields reflects a shift from simplified models to more precise, comprehensive, and adaptable approaches. These advancements have not only driven research in computational chemistry and molecular simulations but also facilitated the application of these methods in fields such as chemical engineering and materials science.

Fitting process of force field parameters

The development and evolution of force fields are closely tied to the parameter fitting method, and accurate fitting techniques are crucial for the precision of the PES constructed using the force field method. Currently, many universal force field parameters exist, particularly in the classical force field, where various combinations of simplified formulae and corresponding parameters have been developed to suit different systems. However, for new materials and complex systems, the accuracy of universal parameters is often limited. In such cases, it becomes necessary to derive the force field parameters specific to the system through a fitting approach. The process of fitting force field parameters generally follows a standard workflow, which can be summarized in three main (1) construction of the initial database, (2) selection of the force field form, and (3) training of the regression model, as shown in Fig. 3. By following these steps, the force field method can be applied to construct the PES for systems ranging from small to large scales.

In the database construction phase, the structures and corresponding energies of the systems are obtained using an appropriate QM method. Since the simulation of PESs is both high-dimensional and complex, and the cost of acquiring high-precision data from QM is substantial, it is important to build the database with broad coverage while keeping the data set as representative as possible. Sometimes, experimental data are also included, such as vibrational frequency data obtained from infrared spectroscopy, which can be used to fit bond and angle parameters within the molecule. Additionally, experimentally measured system properties, such as solubility, elastic modulus, and crystallographic data, can be employed to calibrate the force field parameters. Regarding data volume, ML force fields typically require the largest amount of data, while classical force fields require the least. In addition to manually constructing databases, various training databases have been established during the development of force fields, providing high-quality reference data for accurately fitting the PES. Widely used examples include the QM9 dataset for small organic molecules, ANI datasets (ANI-1, ANI-2x) for organic chemistry coverage, and MD17 for molecular dynamics trajectories of small molecules. In materials science and catalysis, databases such as the MP, OCP, and NOMAD repository offer large-scale atomistic data for diverse material classes and catalytic systems. Additionally, specialized datasets such as the COMP6 benchmark set and the Alexandria database provide rigorous testbeds for evaluating model performance. These curated datasets form the foundation for developing transferable, data-driven force fields capable of capturing complex chemical and physical interactions across broad domains.

The choice of force field form depends on the simulation reactive force fields are preferred for simulating reaction processes, ML force fields are selected for high structural accuracy, and classical force fields are chosen for efficiency in simulations. The regression model training step aims to minimize the difference between the predicted and true values of energy and atomic forces, which is crucial for accurately fitting the force field parameters. The fitting quality is commonly assessed using metrics such as mean absolute error (MAE) or root mean square error (RMSE).

A growing number of benchmark studies have quantitatively compared classical, reactive, and ML force fields in terms of accuracy, generalizability, and computational efficiency. These studies reveal a clear trade-off between accuracy and interpretability, as well as between computational cost and flexibility. Classical force fields, such as Lennard-Jones or harmonic-bond models, typically yield MAEs of 0.5–2.0 eV in energies and >1 eV Å^−1^ in forces when benchmarked against QM references for reactive systems or surface chemistry.^58^ A notable exception is the recently developed GFN force field, which achieves energy MAEs below 0.2 eV.^59^ Reactive force fields, such as ReaxFF, generally outperform most classical models for reactive processes, with reported force MAEs in the range of 0.3–1.0 eV Å^−1^ and energy MAEs around 0.1–0.5 eV, depending on the quality and scope of training.^56^ However, their accuracy varies significantly across chemical systems, and parameter transfer is nontrivial. For example, benchmark studies such as that of Senftle et al. highlight ReaxFF's success in hydrocarbon combustion but also point to its failure in accurately capturing certain transition-metal surface reactions.^30^ ML force fields demonstrate superior performance across a range of systems. Widely benchmarked ANI-1x^60^ and sGDML^61^ models report energy MAEs below 20 meV per atom and force MAEs below 0.1 eV Å^−1^ for equilibrium and near-equilibrium configurations. More recent equivariant models like NequIP^62^ and Allegro^63^ have reduced force MAEs to below 0.03 eV Å^−1^ and energy MAEs <10 meV per atom on diverse datasets, including bulk phases, surfaces, and transition states. These models often approach DFT-level accuracy with orders-of-magnitude speedup. Nonetheless, these models depend heavily on the quality and diversity of the training data, and extrapolation to out-of-distribution chemical environments remains a challenge. Additionally, the OCP provides a large-scale, standardized benchmark for catalyst surfaces and reaction intermediates.^52^ In OCP leaderboard results, ML Force Fields such as GemNet,^64^ Allegro,^63^ and MACE^65^ outperform both classical and reactive models by a substantial margin, achieving force MAEs as low as 0.035 eV Å^−1^ on surface configurations relevant to catalysis, compared to >0.2 eV Å^−1^ for traditional force fields.

The accuracy and efficiency of different regression models vary. For classical force fields, which typically involve a small number of parameters, gradient descent methods can quickly fit the parameters. In contrast, reactive force fields, which involve a larger set of parameters, are more suited to genetic algorithms. ML force fields, using algorithms like neural networks, can approximate any continuous function with arbitrary accuracy. Their strong fitting capability ensures high precision in constructing PESs.

Furthermore, it may be possible to train each one of them on several out-of-equilibrium configurations to make them more robust when performing a more global exploration of the PES. In principle, each can be augmented with out-of-equilibrium training data; however, their ability to incorporate and benefit from out of equilibrium configurations differs. For classical force fields, parameters are usually fitted to equilibrium or near-equilibrium data. The functional forms are rigid and chemically specific, limiting the ability to generalize to far from equilibrium configurations. Extending coverage requires manual re-parameterization. For reactive force fields, parameters can be optimized against broader training sets, including strained geometries and transition states. However, the complexity of cross term interactions can introduce competing minima in parameter space, making it difficult to incorporate high energy configurations without sacrificing low energy accuracy. Achieving robustness away from equilibrium often requires multi-objective fitting and careful weighting of reactive versus non-reactive configurations. For ML force fields, the flexible functional form enables accurate interpolation across diverse regions of configuration space if trained on representative data. Active learning strategies and on-the-fly uncertainty detection can systematically identify and incorporate under-represented high-energy or rare configurations. The latest ML frameworks (e.g., equivariant graph neural networks^62^) can incorporate large, heterogeneous datasets without hard-coded chemistry assumptions, making them particularly suited to integrating off-equilibrium data without manual re-design of the potential.

Applicability and limitations of classical, reactive, and machine learning force fields in heterogeneous catalysis

In heterogeneous catalysis, the accurate modeling of active sites, surface reconstruction, charge transfer, and dynamic electronic effects remains a significant challenge. Different PES modeling approaches (classical, reactive, and ML force fields) offer varying capabilities in addressing these catalysis-specific complexities.

Classical force fields generally rely on fixed atom types, rigid bonding environments, and static charge assignments. While computationally efficient, they are not suited to capture phenomena such as bond rearrangement, changes in oxidation states, or surface reconstruction. So, their applicability is typically limited to non-reactive processes or systems where local electronic structure remains largely invariant.

Reactive force fields, exemplified by ReaxFF, incorporate bond-order-dependent interactions and dynamic charge equilibration schemes (e.g., QEq^66^ or ACKS2),^67^ enabling them to describe bond breaking and formation, as well as approximate charge redistribution. These features allow reactive force fields to model catalytic reactions, including proton/electron transfer, albeit with limited accuracy in systems involving subtle electronic effects or complex charge delocalization. Furthermore, their performance is often constrained by the parametrization strategy and the empirical nature of their functional forms.

ML force fields, trained on high-level QM data, exhibit improved flexibility and generality. When trained on sufficiently diverse datasets including charged systems, surface reconstructions, and transition-state geometries, ML force fields can implicitly learn to model charge transfer, hybridization, and active-site-specific reactivity. Architectures such as message-passing neural networks (e.g., DimeNet++,^68^ Allegro^63^) or equivariant graph neural networks (e.g., NequIP,^62^ MACE^65^) have demonstrated success in modeling chemically complex environments with high fidelity. Some models incorporate explicit charge awareness (e.g., PhysNet,^69^ SpookyNet^70^), while others rely on large training sets to implicitly encode charge redistribution. Despite their promise, many ML force fields still face limitations in modeling long-range electrostatics or variable charge states in heterogeneous systems, and ongoing efforts aim to combine ML potentials with electrostatic or embedding schemes to address these gaps.

Collectively, while reactive and ML force fields represent significant progress toward modeling catalytic complexity, no current approach offers a fully transferable solution across all catalytic scenarios. Future developments will likely benefit from hybrid frameworks that combine data-driven potentials with physical constraints and electrostatic corrections, tailored for the unique demands of catalysis.

Furthermore, in heterogeneous catalysis, the explicit modeling of surfaces introduces several unique challenges beyond those encountered in bulk or molecular systems. Long-range effects including dispersion forces, electrostatics, and image charge interactions play a crucial role in determining adsorption geometries, surface reconstructions, and reaction barriers. Accurately capturing dispersion interactions often requires empirical or semi-empirical corrections, such as the DFT-D scheme of Grimme,^71^ while classical force fields incorporate such effects via parameterized functional forms.^72^ For ML force fields, recent developments in E(3)-equivariant architectures have enabled more explicit incorporation of long-range information into the learned potential.^62^ Another layer of complexity arises from open-shell systems, which are common in transition-metal surfaces, magnetic oxides, and adsorbed radical intermediates. These systems require spin-polarized QM references because their electronic and magnetic states can change along reaction coordinates. Mapping spin-dependent interactions into reactive or ML force fields remains challenging; however, recent ML architectures have begun to address electronic degrees of freedom and spin explicitly. For example, SpookyNet introduces electronic-degree-of-freedom and nonlocal corrections that improve the description of charge and spin dependent effects.^70^ More recent works extend ML potentials to spin and magnetization dynamics, enabling ML modeling of magnetic materials and nonequilibrium spin forces.^73,74^ Incorporating these developments into ML force field workflows is important for accurate simulations of catalytic surfaces that exhibit open-shell character.

To provide readers with a deeper understanding of force fields, the following sections offer a detailed overview of three types, focusing on their functional forms, parameter fitting methods, and applications.

Classical force fields

Forms of classical force fields

For classical force fields, the functional forms typically include the bond term (bond stretching energy, angle bending energy, dihedral angle torsion energy, inversion energy) and the non-bond term (van der Waals term, Coulomb term). An example of a universal force field function is as 1

The first term, bond stretching energy, describes the relationship between two bonded atoms. As atoms oscillate around the equilibrium bond length, a repulsive force is generated when they approach each other, and the energy increases as they move apart, eventually leading to bond dissociation when there is no longer any interaction between the atoms. The second term, bond bending energy, refers to the tension in bond angles and represents the relationship between three bonded atoms. Similar to bond stretching, the bond angle oscillates around an equilibrium angle. The third term, bond torsion energy, describes the interaction between four bonded atoms, specifically the function curve of the dihedral angle and its associated energy. Unlike bond stretching, dihedral angle torsion cannot be represented by simple harmonic functions for most cases. Instead, these periodic changes are typically modeled using Fourier series, which are combinations of trigonometric functions. The fourth term, the van der Waals energy, describes the interaction between non-bonded atoms, often represented by the Lennard-Jones potential. The 12th power term corresponds to repulsion (positive energy), while the 6th power term represents attraction (negative energy). At large distances, the 6th term dominates, creating attractive forces, while at short distances, the 12th term prevails, leading to repulsion. The fifth term, Coulomb energy, models the electrostatic interaction between charged atoms using the classical Coulomb formula, which is inversely proportional to the distance between atoms. For specialized systems, additional terms such as hydrogen bonding, torsional terms, and cooperative terms (such as the Urey–Bradley term in the CHARMM^75^ force field) may be incorporated to improve the precision of the PES.

Different force fields utilize various function forms but can generally be classified into three main periodic function-based, harmonic potential-based, and tabulated potential-based. Over time, numerous force field forms have been developed and applied to different systems (Fig. 4a). Examples include the AMBER^76^ and CHARMM^75^ force fields, which are designed for simulating biological systems; the OPLS force field,^58^ which focuses on simulating condensed phase properties; and the COMPASS force field,^77^ which can predict properties of both gaseous and condensed phases. The DREIDING^78^ and UFF^72^ force fields describe interactions for most elements in the periodic table, although with limited accuracy. Additionally, the newly developed GFN force field strikes a balance between computational efficiency and accuracy, enabling precise large-scale simulations of biomolecular and material systems.^59^

Fitting methods of classical force fields

There are three main sources of classical force field (1) wave number measured experimentally, from which force constants can be calculated using the wave number formula; (2) structural, energy, and frequency data obtained through high-precision calculations, which are then used for parameter derivation and fitting; and (3) force constants derived from the first two methods, which can be further refined using Badger's rule and linear regression. This paper primarily focuses on the second method for obtaining force field parameters.

The most widely used non-fitting method is the Seminario method (Fig. 4b), developed by Jorge Seminario in 1996. This method aims to derive the intramolecular force field parameters, such as the harmonic functions for bond stretching and angle bending, directly from the Hessian matrix obtained by QM calculations.^79^ Currently, several auxiliary tools support the use of this method to fit force field parameters, including the AMBER and AuToFF programs. Below is a brief overview of the method.

For a system with N atoms, the force δF of the 3N component due to a displacement δx can be expressed as a second-order Taylor series 2δF = −[k]δxwhere tensor [k] is the 3N × 3N Hessian matrix. And the [k] used here is defined in unweighted Cartesian coordinates, without mass–weighting. This is because force field fitting aims to reproduce the second derivatives of the PES, which are independent of atomic masses. In contrast, vibrational frequencies and normal modes require mass–weighted and projected Hessians, which are not directly involved in the force field parameterization. Then, eqn (2) can be written as.3

The tensor [k] represents the intramolecular force field up to second order for small displacements δx, and it can be obtained from QM calculation. The eigenvalues λi of [k] correspond to the 3N force constants, which are associated with the 3 translational, 3 rotational and 3N-6 vibrational modes of the molecule. By analogy with eqn (2), the force δFA = (δF~xA~~, δFyA~~,δFz~A~~) on atom A due to a displacement δrB = (δxB, δyB, δzB) of atom B is given by δFA = −[kAB]δrB. The tensor [kAB] can be written as4

The bond stretching constant can be obtained from the projections of the eigenvectors onto the unit vector, û^AB^:5where λ^AB^i and v̂^AB^i are eigenvalues and eigenvectors, respectively.

For the derivation of angle bending constant kθ, first, a unit vector ûN is defined as perpendicular to the plane ABC:6

The unit vectors perpendicular to the bonds AB and CB on the plane ABC can be 7û^PA^ = ûN × û^AB^8û^PC^ = û^CB^ × ûN

The force constants kPA and kPC are then defined as the corresponding constants derived by projecting the eigenvectors of the partial Hessian matrix onto these two 910

And the angle bending constant is then given by11where rAB and rCB are the two bond lengths.

Consider the dihedral defined by the atoms A, B, C, and D which are linked by bonds AB, BC, and CD. The approach used to determine the dihedral force constant is similar to that for bond and angle. The dihedral force constant is given by12

The premise of the Seminario method for deriving force field parameters is that the change in energy associated with the displacement of atom A along the direction û^PA^ will only affect the A–B–C angle. However, in a complex system, neighboring angles may also be altered due to changes of the A–B–C angle. This often leads to an overestimation of bond angle parameters. To address this issue, Allen et al. modified the Seminario method by rescaling the value of kPA with a factor that accounts for the geometry of the molecule.^80^ In addition to the Seminario method, Hirao et al. proposed other rapid force field parameterization schemes, such as partial Hessian fitting (PHF), full Hessian fitting (FHF) and internal Hessian fitting (IHF). These methods are based on the idea of minimizing the Hessian matrix of molecular force fields and quantifying the differences in the Hessian matrix.^81,82^

Numerical fitting methods are commonly applied to fit van der Waals parameters. The process begins by scanning the PES generated from changes in the interatomic distance based on QM calculations. An appropriate van der Waals interaction expression (e.g., Lennard-Jones or Morse) is then selected and used to fit the PES. The fitting process involves iteratively adjusting parameters such as the potential well depth and zero-energy distance until the sum of the squared differences between the two PESs is minimized (eqn (13)). Once the optimization is complete, the van der Waals interaction parameters can be obtained.13

Various methods can be used to adjust the force field parameters to minimize J(θ), with one of the most straightforward approaches being gradient descent. The process primarily involves two primary first, an initial set of force field parameters is provided, which may be randomly assigned; second, the parameters are iteratively adjusted in the direction of the negative gradient until the loss function converges to a predefined value. As illustrated in Fig. 4c, the red areas represent higher values of J(θ), and as the parameters are updated, the value of J(θ) decreases, eventually reaching the deep blue regions. The core idea behind gradient descent is to move in the direction of the negative gradient from the current position. To further accelerate convergence, the method was further extended into the steepest descent method, where the step size for parameter adjustment reduces as the gradient diminishes. This leads to slower progress as the optimization nears the target value.

However, the final point reached by gradient descent may not always correspond to the global minimum, but rather to a local minimum, especially when dealing with force fields that involve a large number of parameters. To overcome this challenge, a common strategy is to initialize multiple sets of parameters and perform gradient descent along several different paths. The optimal parameters are then selected from these paths, allowing for a more robust approach to global optimization (Fig. 4c).

Atomic charge is a critical component of classical force field parameters, and various methods have been developed to fit atomic charges, including Merz-Kollman,^83^ CHELPG,^84^ and RESP.^85^ Among these, the RESP method, proposed by Kollman et al., effectively addresses issues such as conformational dependence, numerical instability, and atomic equivalence of internal atomic charges, making RESP charges widely used. The principle behind RESP involves iteratively adjusting atomic charges using the least squares method until the error between the calculated classical potential and the QM potential is minimized. To avoid unreasonable charge distributions, the RESP method introduces constraints, such as the neutral constraint (which ensures the sum of molecular charges is zero) and a penalty term that controls the charge distribution. Programs such as AmberTools^86^ and Multiwfn^87^ make it easy to fit RESP charges.

Applications of classical force fields

In this section, we highlight selected studies that utilize classical force field methods, focusing on their application in heterogeneous catalysis rather than providing a comprehensive literature review. Specifically, we explore how classical force fields are used to study catalyst interface properties, particularly in areas like metal particle morphology and molecular diffusion. These applications demonstrate the strengths of classical force fields in efficiently simulating the dynamics of heterogeneous interfaces.

The morphology of metal catalysts (specifically the size, shape, and distribution of metal particles) plays a crucial role in determining catalytic performance.^20,88–91^ The size and shape directly influence the number of active sites on the catalyst surface. For example, Hu et al. investigated the size effect on the atomic structure of amorphous systems originating from a Cu64Zr36 particle (containing 50–5000 atoms) using MD simulations with an embedded atom method (EAM) potential.^92^ Their findings show that particle size strongly impacts the local atomic structure, with Cu64Zr36 particles exhibiting core–shell structures. The shell component of the particle has a lower average coordination number, shorter bond lengths, higher ordering, and lower packing density compared to the core, due to Cu segregation on the shell. Given that metal nanoparticles often exhibit high surface energy, they tend to aggregate during catalytic processes, which can lead to catalyst deactivation. To mitigate this, metal nanoparticles are commonly dispersed on support surfaces to enhance stability.^93^ In a study by Wang et al., the effect of metal–support interactions on the stability of metal nanoparticles was examined using a Morse force field for parameterization.^94,95^ Their simulations revealed that the melting point of Pt nanoparticles supported on a substrate was significantly higher than that of unsupported nanoparticles, and the melting point increased as the metal–support interaction strengthened. In addition to improving stability, the support material also affects the morphology of the metal nanoparticles. The simulations indicated that increasing the metal–support interaction leads to a higher number of unsaturated coordination atoms in the nanoparticle. Notably, when Pt nanoparticles were supported on bare MXene, a film-like structure formed on the substrate surface. Wang et al. further explored this by tuning the nanoparticle structure through variations in the surface functional groups of MXene.^96^ Their results showed that when the surface functional groups on Nb2C MXene transitioned from –Cl, –Br, and –O to partial –O, the supported Pd561 nanoparticle exhibited a distinct morphological shift from 3D to 2D (Fig. 5a), consistent with electron microscopy observations (Fig. 5b). Initially, a monolayer of Pd preferentially formed on the exposed Nb sites, followed by the creation of a second Pd layer upon encountering oxygen functional groups, ultimately exposing the (111) facet. Based on these findings, they constructed Pd metalenes supported on Nb2C MXenes, which demonstrated efficient selective hydrogenation of phenylacetylene at room temperature.

Classical force fields are also commonly employed to simulate the dynamic behavior of molecular systems, including the diffusion, adsorption, and desorption processes of reactants, intermediates, and products, and their impact on catalyst performance.^97–100^ To improve the catalytic performance of N2 electroreduction (NRR) on metal surfaces, hexanethiol (HEX) was selected as a modifier to inhibit the competitive adsorption of water molecules in the hydrogen evolution reaction (HER). Chen et al. studied the effects of HEX on the diffusion and adsorption behaviors of H2O and N2 molecules on the Cu surface through MD simulations.^101^ The study involved 100 nitrogen molecules and 5000 water molecules, based on the solubility of N2 in water. The results revealed that on the HEX-modified Cu surface, N2 molecules could diffuse to the catalyst surface, while H2O molecules were blocked by the HEX layer (Fig. 5c). In contrast, the non-HEX Cu surface is predominantly covered by H2O molecules. After HEX modification, the number of N2 molecules adsorbed on the Cu surface increased significantly from 24 to 48, while the number of H2O molecules adsorbed decreased from 700 to 23 (Fig. 5d). Additionally, the calculated potential of mean force indicated that the presence of HEX notably raised the equilibrium position of H2O (from 0.18 nm to 0.54 nm). After HEX modification, the Faraday efficiency of various metal catalysts (Cu, Au, Pd, Pt, and Ni) prepared by experimental collaborators was significantly enhanced, confirming the hydrophobicity of HEX on the metal surface and its effective role in promoting N2 adsorption.

For the diffusion properties of molecules in porous materials, Xu et al. performed MD simulation on the gas diffusion in the interlayer gallery of graphene and graphene oxide (GO) membranes, and elucidated the mechanisms of gas selective separation (Fig. 5e). They found that both the layer spacing and the chemical modification of the membrane surface significantly influenced the selective gas penetration.^43^ For example, for GO membranes, the He/CO2 selectivity can reach as high as 30, compared to 4.5 for CH4/CO2 and between 2 and 3 for CO/CO2, N2/CO2, and O2/CO2. In contrast, for graphene layers, the selectivities for N2/CO2, CO/CO2, and O2/CO2 are relatively low, ranging from 1 to 3 (Fig. 5f). To further explore the correlation between pore properties and molecular diffusion coefficients, Wang et al. studied the interlayer diffusion of CO in graphene nanosilts with different Pt loading. They found that the diffusion energy barrier is related to the distribution of CO molecules in the system, which in turn is influenced by factors such as temperature, pressure, interlayer distance, and the properties of the supported metal. High temperature, low pressure, and fewer surface atoms were found to facilitate the diffusion of gas molecules. To quantify the relationship between the diffusion coefficient and the environmental and structural properties of the system, a generalized formula for confined diffusion of CO in the supported system is derived based on the simulation data.^97^ Further simulations of interlayer diffusion of CO2 were conducted to study the effects of the metal–molecule interaction on the adsorption and diffusion properties of molecules.^44^ The results showed that the difference in the adsorption amounts of CO and CO2 on Pt nanoparticles correlates with the adsorption site and exhibits a volcanic trend with respect to temperature. The metal–gas interaction and the surface atomic number of metal nanoparticles have great influence on the diffusion of gas molecules, especially at low temperatures. The simulation data further validate the accuracy and generality of the generalized formula.

Compared to the above classical force fields, the newly developed GFN-FF force field demonstrates superior generality and accuracy, approaching the precision of QM methods in certain cases. For instance, metal–organic frameworks (MOFs), widely used in heterogeneous catalysis, typically have few available force fields (such as UFF), which are often inadequate for describing conjugated systems and metal coordination. In contrast, GFN-FF enables accurate geometric optimization of MOF structures. For example, for a Goldberg polyhedra composed of 46 Pd^2+^ ions and 96 organic ligands, with a total of 3888 atoms, the RMSD of the heavy atoms in the GFN-FF optimized structure is only 0.75 Å, which is in excellent agreement with experimental results. Similar accuracy can be achieved for other metal–organic polyhedra and MOFs.^59^ Additionally, water, as the most common liquid, is notoriously difficult to describe accurately, leading to the development of various force field models, including rigid, flexible, and polarizable force fields. Molecular dynamics simulation of water using the GFN-FF force field shows that, while the calculated self–diffusion coefficient (2.06 × 10^−5^ cm^2^ s^−1^) is remarkably close to the experimental value (2.35 × 10^−5^ cm^2^ s^−1^), its density is significantly overestimated (1.23 g cm^−3^ compared to the experimental value).^102^ Therefore, while GFN-FF offers a fast and accurate classical force field, like all universal force fields, it may not always provide the most accurate results.

From the above cases, it is evident that classical force fields remain valuable tools for modeling heterogeneous catalysis, particularly when large system sizes and long simulation times are required. Their computational efficiency enables exploration of phenomena such as adsorbate diffusion, nanoparticle sintering, and solvent-mediated surface restructuring over length and time scales inaccessible to QM methods. In supported catalyst systems, classical force fields have been successfully applied to investigate nanoparticle morphology evolution, interfacial interactions between the support and active phase, and adsorption–desorption equilibria under realistic reaction conditions.

However, the inherent limitations of classical force fields arise from their fixed topologies, predefined interaction functional forms, and parameter sets that are typically tuned for narrow chemical spaces. These constraints hinder their ability to capture bond formation/breaking, charge transfer, and polarization effects, which are central to catalytic processes. Furthermore, the lack of explicit electronic degrees of freedom makes it challenging to describe surface reactions involving variable oxidation states, adsorbate-induced reconstruction, or metal-support charge redistribution. Transferability across different surfaces, phases, and chemical environments is often limited, necessitating reparameterization for each new system.

Future developments should aim to improve the accuracy and transferability of classical force fields for catalysis by integrating more flexible interaction terms (e.g., polarizable force fields), incorporating long-range electrostatics more rigorously, and systematically coupling classical descriptions with reactive or ML force fields in hybrid schemes. Such approaches could preserve the efficiency of classical force fields while extending their applicability to the complex, reactive, and dynamically evolving environments characteristic of heterogeneous catalysis.

Reactive force fields

Forms of reactive force fields

The reactive force field is typically based on bond order to construct the system's PES. In the early stage of development, its form was designed specifically for the bonding situation in a specific reaction system, making it difficult to generalize. For example, based on the relationship between Pauling bond length and bond order, Johnston developed the BEBO force field to study the reactive PES of the H + H2 system.^103^ Besides the BEBO force field, several other reactive force fields have been developed to address bond formation and breaking in specific systems. The AIREBO (Adaptive Intermolecular Reactive Empirical Bond Order) potential,^104^ an extension of the REBO model,^105^ is designed for hydrocarbons and carbon-based materials and has been widely used in simulations involving graphene, CNTs, and organic reactions. COMB3 (Charge-Optimized Many-Body) potentials^106^ enable the simulation of metal–ceramic interfaces and allow for dynamic charge transfer using a charge equilibration scheme, making them suitable for oxide-based catalysts and interfacial systems. The Environment-Dependent Interatomic Potential (EDIP)^107^ is tailored for covalent materials, particularly silicon, and can model defect formation and surface reconstructions with relatively low computational cost. While these models have not yet seen widespread application in complex catalytic reactions, they offer valuable tools for studying specific materials and interfaces under reactive conditions.

Currently, ReaxFF, developed by van Duin et al. in 2001,^56^ is the most widely used reactive force field due to its parameter transferability. The initial version of ReaxFF (2001) was focused on hydrocarbons and used the same dissociation energy for C–C single, double, and triple bonds, which worked well for hydrocarbons but was limited for more complex systems. In 2003, the ReaxFF functional form was extended to systems containing Si, O, and H, with separate parameters for single, double, and triple bond dissociation energies. This extension also introduced a lone-pair energy term to handle the formation and dissociation of oxygen lone pairs. The 2003 extension was further improved to handle more complex group chemistries, such as the conjugation term for –NO2 group chemistry in nitramines and a triple-bond stabilization term for better describing terminal triple bonds.^108^ By 2005, the ReaxFF functional form stabilized,^109^ and the general form now includes terms 14Esystem = Ebond + Eover + Eunder + Eval + Etors + Econj + Evdw + ECoulomb + EH–bond + Erestwhere the terms of bonding interactions include bond, valence angle, lone pair, conjugation, and torsion angle, and the terms of nonbonding interactions include van der Waals, Coulomb interaction, and hydrogen bond. In addition, to deal with special systems, additional terms will be introduced into the formula, such as angle bending terms for Mg–Mg–H zero-degree angles^110^ and double-well angular terms for aqueous transition metal ions.^111^

As mentioned above, the ReaxFF is categorized into bond-order-dependent and bond-order-independent contributions. Bond order (BO) a function of interatomic distance rij, which is divided into the contribution of single bond (BO^σ^ij), double bond (BO^π^ij), and triple bond (BO^ππ^ij):15

Based on the bond order formula, interactions such as bond and valence angle terms can be derived through a series of transformations; for details, see the work of van Duin et al.^109^ Additionally, this bond order formula accounts for long-distance covalent interactions in the transition state, enabling the ReaxFF to accurately predict reaction barriers.

For the Coulomb interaction, the QEq method is used to calculate and adjust the charge distribution of each atom in the system.^66^ This method achieves charge distribution through electronegativity equilibration, accounting for both atomic interactions and the dynamic changes occurring during chemical reactions. During a reaction, the electronegativity of each atom varies as its local environment changes. Atoms with higher electronegativity tend to attract electrons, while those with lower electronegativity tend to lose electrons.

Initially, the ReaxFF parameters were developed exclusively for organic systems.^56^ Subsequently, the parameters gradually incorporated metallic elements^110,112,113^ and other non-metallic elements,^108,114^ further expanding the applicability of ReaxFF to a wider range of systems. Currently, the ReaxFF method has been successfully applied to simulate various reaction dynamics, including hydrocarbon organic small molecule systems,^56^ polymer systems,^115,116^ high energy material systems,^117,118^ metal oxide systems^119,120^ and transition metal catalyst systems.^113,121^ Rapid reaction processes, such as explosion and combustion, can also be simulated using ReaxFF parameters.^117,118,122^ Additionally, the ReaxFF-based Monte Carlo reaction kinetics method has been employed to investigate experimental structures that are difficult to resolve experimentally,^119,120^ catalytic reactions in fuel cell electrode materials,^123,124^ and catalytic processes in porous materials.^114,125^ To date, ReaxFF parameters for various elements in the periodic table have been developed, as shown in Fig. 6a.^30^ However, in the specific study, the transferability of these parameters is limited, and it is impossible to simply combine these parameters to get satisfactory results. Depending on the different O/H atomic and bond parameters, there are currently two main categories of ReaxFF parameter sets that are intra-transferable with one the combustion group and the aqueous group. Additionally, there are several independent groups whose parameters often require more extensive refitting. For specific details on grouping, refer to the summary by Thomas et al.^30^

Fitting methods of ReaxFF parameters

Compared to the classical force field, the ReaxFF function is more complex and has a greater number of parameters, making it more challenging to develop. During the development of ReaxFF, it is often not possible to fit all parameters. Some validated force field parameters can be obtained from existing parameter libraries or theoretical literature, and relevant parameters can be selected for fitting based on the system under study. Two main methods are commonly used for the fitting successive one-parameter search and genetic algorithm.

Successive one-parameter search is the earliest method used for ReaxFF parameter fitting.^126^ The method assumes a parabolic relationship between a single parameter value and the total fitting error, determining the optimal parameter value by calculating the total error at three distinct parameter values. However, because most parameters in ReaxFF are interdependent, a change in one parameter will cause the optimal values of others to shift. Therefore, the parameters must be optimized iteratively to ensure that the fitting error is minimized. This method offers certain advantages for fitting processes with many parameters, as the fitting process can be interrupted if a parameter appears to have an unrealistic value. However, its fitting efficiency is relatively low.

For ReaxFF with multiple parameters, in addition to the fitting efficiency issue, the fitting process is a high-dimensional, non-separable optimization problem with multiple minima. As a result, the fitting results often converge to local optima. While deterministic global optimization techniques are available, they encounter significant practical challenges in high-dimensional search spaces and with computationally expensive objective functions, both of which are typical of the present problem. Viable alternatives include nondeterministic search heuristics, such as genetic algorithms (GA).^127^

GA are optimization techniques inspired by the process of natural evolution. The two primary concepts underlying genetic algorithms are natural selection and genetic dynamics. Natural selection involves choosing individuals with greater fitness to pass their traits to the next generation, based on superior performance. The goal of the genetic algorithm is to maximize fitness as the generations progress. A simplified depiction of the genetic algorithm, illustrating the key steps, is shown in Fig. 6b. First, based on parameters from the literature, multiple sets of parameters are generated through mutation operations to form the initial population, with each set representing an individual. The fitness function is then used to assess the quality of each individual, typically defined by the difference between the system energy calculated from the parameters and the system energy obtained from QM calculations. Next, a selection operation is applied to the current population, with methods such as rank-based selection used to retain individuals with higher fitness. These selected high-fitness individuals undergo crossover operations simulating genetic recombination to produce offspring. To enhance population diversity and prevent premature convergence to a local optimum, mutation operations are employed to randomly alter specific force field parameters, thereby generating new individuals. By applying selection, crossover, and mutation operations, a new generation is created to replace the old one, followed by a fitness evaluation of the new population. After each generation, a check is made to determine whether the termination criteria have been met. Common termination criteria include reaching the maximum number of iterations or achieving a predefined fitness threshold. If the termination criteria are satisfied, the individual with the best fitness in the current population is selected as the final solution.

To enable efficient and rapid parallel optimization of parameters, Goddard et al. developed GARFfield (genetic algorithm-based reactive force field optimizer method), a hybrid multi-objective Pareto-optimal parameter development scheme combining genetic algorithms, hill-climbing routines and conjugate-gradient minimization.^128^ The software architecture is illustrated in Fig. 6c. The core of the GA search mechanism is stochastic gradient search, which requires the evaluation of the objective function's fitness for each generation. A drawback of this method is the significant time cost associated with the repetitive search and fitness evaluation. To address this, a hybrid algorithm combining artificial neural networks (ANN) and GA has been proposed and accepted for fitting the ReaxFF parameters.^129,130^ The ANN read and analyze the data and total error values generated during the GA process, thereby enhancing the efficiency and reducing the time cost of the heuristic search.

Applications of ReaxFF

Modeling catalytic processes requires consideration of bond breaking and formation, yet the time and length scales inherent in nanoscale interface reactions cannot be addressed with QM. The bond order concept, coupled with the low computational cost, enables ReaxFF to bridge the gap between QM and non-reactive force fields. In this context, we focus on the application of ReaxFF to metal catalysts and their oxides to highlight the methodological strength of ReaxFF: modeling reactive chemistry at heterogeneous interfaces, rather than providing a comprehensive review of the literature.

Transition metals are widely used in heterogeneous catalysis reactions due to their unique electronic structures, with their geometric structures closely influencing their properties.^131–134^ The simulation scale achieved by ReaxFF is well-suited for modeling catalytic processes on metal surfaces and clusters, particularly at defect sites or unsaturated sites. Iron-based catalysts play a significant role in the water–gas shift and Fischer–Tropsch synthesis reactions. Based on ReaxFF fitted by genetic algorithm, Wen et al. studied the structure–activity relationship of Fe nanoparticles in CO activation.^135^ The results showed that CO dissociation can be effectively promoted by introducing line dislocation and vacancies on Fe nanoparticles. Furthermore, four mechanisms of CO2 formation catalyzed by Fe nanoparticles were analyzed through adsorption and activation of surface carbon in MD simulation trajectories (Fig. 7a–d). They found that, at the initial stage of the reaction, CO molecules adsorb on the surface of Fe nanoparticles and dissociate. At this point, the oxygen concentration on the surface is low, and the probability of CO2 formation via the Langmuir–Hinshelwood (LH) mechanism (where adsorbed O atoms react with gaseous CO molecules) is minimal. However, some CO molecules adsorbed on the Fe surface without dissociation are more likely to react with gaseous CO molecules, forming CO2 through the Eley–Rideal (ER) mechanism. As the reaction progresses, O atoms accumulate on the surface of the nanoparticles, and the surface also contains undissociated CO molecules, which can form CO2 through the LH mechanism. Additionally, adsorbed undissociated CO molecules can also react with each other to form CO2. According to the above analysis, it can be concluded that both the ER and LH mechanisms play crucial roles at different stages of CO2 formation catalyzed by Fe, with the ER mechanism primarily dominating in the early reaction stages, while the LH mechanism becomes more prominent in the later stages. For metal-catalyzed carbon nanotube (CNT) growth, the process involves the dissolution and migration of carbon atoms within the metal. CNT growth is initiated by the adsorption and dissociation of hydrocarbons on the metal surface. To investigate the factors influencing this process, Mueller et al..developed a ReaxFF reactive force field to describe the interaction between hydrocarbons and the Ni surface.^121^ They found that surface defects play a key role in the decomposition rate of CH3, particularly in the final step, where CH decomposes into C and H. Furthermore, due to the transferability of ReaxFF parameters, the Ni/C/H parameter set was utilized by Neyts et al. to elucidate the influence of ion bombardment on CNT formation.^136^ The simulation results showed that the energy of the impacting ion can be adjusted to break low-coordination C–C bonds, leading to the formation of new bonds in the network, which facilitates guided growth.

The ReaxFF formula combines bond order and charge transfer formalisms, making it particularly well-suited for describing the evolution of the reaction process in metal oxide materials, which exhibit both covalent and ionic interactions. For example, Chenoweth et al. developed V/O/C/H ReaxFF parameters aimed at describing the interaction of hydrocarbons with vanadium oxide.^119,120,137^ Based on the MD simulation with this ReaxFF parameter set, they investigated the dissociation process of methanol on the surface of V2O5. On a defect-free surface, C–H dissociates more easily than O–H. On the defect surface, the oxygen atom in methanol binds to the reduced vanadium defect site, lowering the energy barrier of O–H dissociation. Additionally, they observed the desorption of a water molecule formed on the hydroxylated V surface, which leads to interlayer bonding between two metal atoms of different oxidation states, transitioning from V^V^ and V^III^ sites to two V^IV^ sites. Metal oxide materials typically exhibit partial, mixed, and irregular metal occupations at various crystallographic sites. To obtain the lowest-energy structure, researchers employ the combined Monte Carlo method to explore possible metal oxide configurations. For example, Janik et al. used ReaxFF-based Grand Canonical Monte Carlo (GCMC) simulations to construct thermodynamically stable Pd/CeO2 surface models under reaction conditions.^138^ Based on the equilibrium configuration, they further simulated methane activation using reactive MD to evaluate the catalytic performance of different interface morphologies. By counting the number of reactive species, it was found that a sharp decrease in the number of gaseous CH4 molecules was accompanied by an increase in adsorbed H atoms at approximately 1440 K, indicating that methane was activated (Fig. 7e). By comparing methane activation on Pd and PdO~x~ clusters that are both supported on and embedded in the CeO2 lattice (Fig. 7f), they found that the supported metal clusters facilitate methane activation more efficiently, with C–H bond breaking occurring at the unsaturated coordination sites on the cluster edges. But in the embedded model, methane activation is slower because CH4 has limited access to these sites. In contrast, PdO~x~ clusters in the embedded structure are more likely to activate methane. This can be attributed to the reduced exposure of edge sites in the embedded clusters, while PdO~x/CeO2~ interface mixing generates unique active sites (Fig. 7g). These studies demonstrate the utility of ReaxFF in modeling catalytic processes on metal and metal oxides, as it can simulate large-scale reaction processes that are not achievable through QM.

In summary, reactive force fields such as ReaxFF, provide a powerful framework for simulating bond-breaking and bond-forming events at scales relevant to heterogeneous catalysis. They have been widely applied to investigate surface reaction mechanisms, catalyst activation and deactivation, support-metal interfacial chemistry, and degradation pathways under realistic temperature and pressure conditions. By explicitly accounting for variable bond orders and dynamic charge equilibration, reactive force fields enable the exploration of catalytic cycles beyond the reach of classical force fields, while maintaining computational tractability for systems containing thousands of atoms.

Despite these strengths, several limitations hinder their predictive reliability. The accuracy of reactive force fields strongly depends on the breadth and quality of the parameterization dataset, which must adequately sample relevant reactive configurations, surface reconstructions, and intermediates. Transferability remains a challenge as parameters optimized for one catalytic material or reaction type may not generalize to others without significant re-fitting. Current charge equilibration schemes (e.g., QEq) often fail to capture non-local charge transfer, polarization under strong electric fields, or complex redox processes, limiting their applicability in electrocatalysis and photocatalysis. Furthermore, the high dimensionality of parameter space makes systematic optimization difficult, and the lack of rigorous error estimation complicates the assessment of model reliability.

Future directions for reactive force fields in catalysis include the development of next-generation reactive potentials with improved electrostatic and polarization models, the incorporation of machine learning-assisted parameter optimization to accelerate and improve transferability, and the construction of large, diverse, and surface-specific training datasets. Hybrid simulation schemes that couple reactive force fields with on-the-fly QM calculations or ML force fields may offer a promising strategy to combine chemical reactivity with long-timescale dynamics, thereby bridging the accuracy-efficiency gap for realistic catalytic environments.

Machine learning force fields

Forms of machine learning force fields

The application of ML methods in the construction of PESs dates back to 1992, when Bobby et al. innovatively introduced neural networks to model and predict the energy of molecular systems, thereby describing the behavior of polymer molecules in vibrational modes.^139^ Subsequently, in 1995, Blank et al. used neural networks to approximate the PES for hydrogen formation reactions on silicon surfaces.^57^ This work marked a significant milestone in PES modeling and paved the way for the application of ML methods to model complex chemical reactions.^140–142^

Unlike classical force fields and reactive force fields, which rely on specific formulae, ML force fields do not have a predefined functional form. Instead, ML force fields use atomic local environments as descriptors and are trained using ML algorithms, such as neural networks or Gaussian regression, to construct the PES of a system. By moving away from intuitive, physics-based expressions, ML force fields can offer a more accurate PES.

Designing appropriate descriptors to represent the local environment of atoms is crucial for the effectiveness of ML force fields. The local chemical environment refers to the interactions between an atom and its neighboring atoms within a specified cutoff radius rcut, as illustrated in Fig. 8a. Behler et al. introduced the concept of locality approximation,^33,143^ which asserts that atomic interactions are primarily determined by the local chemical environment of the central atom, rather than by the influences of all other atoms in the system. This approximation reduces computational costs while maintaining the model's transferability.

Building upon the locality approximation, the ML force field further assumes that the total potential energy E of the system can be expressed as the sum of the energies Êi of individual central atoms, as shown in eqn (16). The energy Êi is determined by the local chemical environment of atom i.16

In a neural network force field, each atom is represented by a neural network (Fig. 8b), with the input corresponding to the information about its local chemical environment and the output representing the atomic energy term Êi. The total potential energy E of the system is obtained by summing the outputs of the neural networks for all atoms. This approach constructs the model at the atomic level rather than for the entire system, ensuring its scalability. It is important to emphasize that the atomic potential energy Êi is not a predefined label but rather an intermediate variable introduced based on the potential energy decomposition assumption. The total potential energy E is the direct target of model training.

Although the types and coordinates of atoms in the system can fully describe the PES, directly using these as inputs to the ML force fields is not appropriate. On the one hand, atoms do not have absolute coordinates, and the overall translation and rotation of the system, while altering the coordinates of individual atoms, do not affect the system's potential energy. On the other hand, regarding atom types, exchanging positions of identical atoms in the system has no effect on its properties, meaning that changing the atomic index order does not influence the system's behavior. Therefore, it is necessary to transform the atomic coordinates and type information into suitable descriptors for ML, which can then be used as inputs to the ML force field model.

Therefore, the descriptors should adhere to the principle of symmetry, meaning they must be invariant under translation, rotation, and permutation of atomic indices. More precisely, the descriptors should have a bijective relationship with the atomic structure and type information, ensuring that for any given system configuration, there is a unique corresponding descriptor, and different structures have distinct descriptors. In addition to fulfilling these requirements, a good descriptor should also exhibit continuity, low computational cost, and high representational efficiency. Several methods for constructing descriptors have been proposed, including the Smooth Overlap of Atomic Positions (SOAP),^144^ Coulomb Matrix (CM),^145^ Atom-Centered Symmetry Functions (ACSF),^143^ Spherical Harmonics,^146^ Many-Body Tensor Representation,^147^ Bispectrum features,^148,149^ and others.

In addition to atomic environment descriptors, another crucial aspect of the ML force fields is the ML algorithms themselves. With advancements in ML techniques and the widespread use of GPU computing, numerous ML force field models have been developed, as shown in Fig. 8c. Early ML force fields primarily relied on traditional ML techniques, such as Gaussian Approximation Potentials (GAP)^148,149^ and fully connected neural network architectures like Atom-Centered Symmetry Functions (ACSF).^143^ As algorithms have advanced, more sophisticated ML force fields have emerged, including GDML,^61^ ANI,^150,151^ DPMD^46^ and others.^152–158^ Recently, ML force fields based on graph neural network (GNN) architectures have gained traction, with models like NequIP^62^ and Allegro.^63^ These models do not require explicit descriptor construction; instead, they represent the entire system as a graph and incorporate the effects of the chemical environment through a message-passing mechanism. Given the wide variety of ML force fields, this review focuses on introducing the more widely used neural network force field. This emphasis does not imply that other models are less effective. For example, GAP based on Gaussian process regression, sGDML rooted in nuclear ridge regression, and other linear and kernel-based methods such as Moment Tensor Potential (MTP)^159^ and Atomic Cluster Expansion (ACE),^160^ offer interpretable and systematically improvable representations of the PES. These approaches typically require less training data than deep neural networks and permit rigorous error estimation, making them particularly appealing for catalytic studies where generating high-fidelity training datasets is computationally expensive.

Fitting methods of neural network force fields

Neural networks are mathematical models that mimic the structure and functioning of neurons in the human brain. Typically, these models consist of an input layer, hidden layers, and an output layer, with each hidden layer containing a set of interconnected neurons. Fig. 8b illustrates a simple neural network model, where each input is progressively transformed through neurons in the hidden layers before reaching the output. The conversion between neurons is achieved by the following 17where x represents the input from the preceding layer, ω denotes the weights, b is the bias term, and f signifies the activation function.

In the neural network force field, the input consists of the localized chemical environment information of each atom, and the corresponding atomic energy is obtained through the mapping between neurons, as shown in eqn (18)–(20).181920

A notable approximation in this method is the restriction of atomic interactions within a cutoff sphere. The resulting short-range potential is well-suited for describing local bonding, even in complex atomic environments. However, for many systems, long-range interactions (e.g., electrostatic and dispersion forces) are also important. To address this issue, like the environment-dependent atomic energies, the atomic charges also depend on vectors of ACSFs that describe the atomic environments, and are obtained as outputs of the atomic charge neural networks. These charges are then used to calculate the long-range electrostatic energy using standard methods, such as Coulomb's law or the Ewald summation. PhysNet^69^ is a prototypical example of this class of methods. Although this method incorporates electrostatic interaction calculations, it relies on atomic charges determined solely by the local chemical environment, which limits its accuracy in describing non-local effects. To overcome this limitation, methods such as the charge equilibration neural network technique (CENT)^161^ have been developed to determine atomic charges based on the global environment, effectively addressing non-local interaction challenges. This method also uses ACSF vectors to describe the atomic environment as input to the neural network. The key difference lies in its output, which is the electronegativity values that are subsequently used in a charge equilibration scheme to determine atomic charges.

Compared to traditional neural network force fields that rely on manually designed fixed descriptors, the graph neural network force fields developed recently demonstrate superior generalization capabilities and dynamically adapt to varying topological structures. NequIP exemplifies this approach by employing an E(3)-equivariant architecture to accurately model atomic interactions.^62^ In this framework, atoms serve as nodes, while edges represent interactions between a central atom and its neighbors within a cutoff radius, forming an atomic graph. NequIP operates directly on these graphs, with features that transform equivariantly under three-dimensional rotations, translations, and reflections. This geometric equivariance allows the network to inherently respect physical symmetries, enhancing data efficiency and prediction accuracy. Through multiple layers of equivariant message passing that update node and edge features, NequIP effectively captures complex many-body interactions and anisotropic effects essential for realistic molecular and materials simulations. Similar to NequIP, methods such as MACE^65^ and Allegro^63^ are also graph-based architectures that model interatomic interactions through message passing schemes, enabling the explicit treatment of many-body effects.

During the training process of a neural network, the weights Wi and biases bi are initialized randomly and optimized to minimize a loss function that quantifies the difference between the predicted values and the reference data. The primary focus is typically on two the system's potential energy and the atomic forces. Commonly used loss functions include the MAE and RMSE. The smaller the MAE or RMSE on the test set, the more accurate the model's predictions.

In some cases, additional terms are included based on the system's properties; for example, for solid material systems, a stress-related loss function term may be added. In addition to the weights Wi and biases bi that are determined during training for a given data set, model training also involves selecting hyperparameters, such as the number of hidden layers and neurons per layer. In theory, if the number of hidden layers and neurons per layer is sufficiently large, the model can output energy with arbitrary precision. However, too many hidden layers and neurons can significantly increase the computational cost. Therefore, selecting appropriate hyperparameters is crucial for the accuracy and efficiency of the model. Hyperparameters are typically optimized using exhaustive search schemes like grid search or random search, often combined with informed guesses for suitable search ranges. Currently, for many hyperparameters, model performance remains fairly robust for small changes, and defaults perform well across different data sets.

Before optimizing any hyperparameters, the test set must be separated from the available reference data, and the remaining data is then divided into training sets and validation sets. This separation is essential because, in force field applications, the structures are not included in the training data, so the model must be trained to be capable of predicting unseen data. Therefore, for each trial combination of hyperparameters, the model is trained on the training data, and its performance is evaluated on the validation set to estimate the generalization error. The best-performing model is then selected. To enhance the model's generalization ability, a k-fold cross-validation method is also used. For more complex neural network models, additional considerations are required, and further details can be found in the summary by Tokita et al.^162^

Training ML force fields often relies on large amounts of high-quality data, which typically demands expensive QM calculations. To reduce training costs, several methods for optimizing data acquisition have been developed. One approach is to generate initial data using less expensive classical molecular dynamics simulations, followed by selecting a minimal subset for more computationally expensive ab initio calculations.^163^ Additionally, strategies such as active learning and enhanced sampling can further improve efficiency.^164^ Active learning significantly reduces data requirements by intelligently selecting the most informative data points for labeling, while enhanced sampling techniques (e.g., meta-dynamics, umbrella sampling) can effectively capture rare events and improve the coverage of reaction paths. By combining these methods, data-efficient ML force fields can be developed, enabling construction of high-precision PESs with limited data conditions.

The ML force field models typically rely on local structural information, such as interatomic distances and angles, to predict the PES. However, the traditional model tends to overlook long-range electrostatic interactions between charged particles in the system. As a result, traditional models based on local structural features still faces challenges when studying systems (such as clusters, interfaces, and gas-phase systems) that are significantly influenced by long-range electrostatic interactions. To address this, researchers have incorporated a long-range electrostatic interaction term into the traditional model, learning the hidden variables in the local atom descriptor and applying Ewald summation to account for these interactions.^165–167^ This approach has been successfully validated in systems such as charged polar molecules and biological macromolecules. Additionally, another challenge with ML force fields is that the system's structural features can change under different environmental conditions. For example, variations in coordination atoms, temperature, pressure, or solvents can significantly affect the structure and interactions of molecules, making it challenging for a model trained under one set of conditions to transfer efficiently to another. In contrast to the parameter transferability observed in classical force fields and the ReaxFF force field for similar systems, ML force fields, particularly those based on neural networks, generally lack parameter transferability. As a result, the model needs to be retrained for new systems.

To improve learning efficiency and accuracy on new systems or tasks while reducing reliance on expensive high-quality training data, transfer learning strategies have proven effective. This process mainly includes pretraining followed by fine-tuning, incremental learning, and active learning–assisted transfer. Typically, foundation models, such as MACE-MP-0,^168^ are first pretrained on large-scale general datasets, like the open catalyst project, to capture universal physicochemical features. They are then fine-tuned to rapidly adapt to specific systems, thereby reducing the demand for high-quality training data. Additionally, incremental learning enables continuous model updates as new data become available, preventing the forgetting of previously learned knowledge. Active learning incorporates uncertainty evaluation to selectively sample high-value data, enhancing training efficiency.^169^ These strategies, when combined, effectively facilitate the rapid deployment and efficient application of ML force fields across diverse catalytic and materials systems.

Applications of machine learning force fields

The advantage of ML force fields lies in their ability to achieve accuracy comparable to that of QM, while being significantly faster than QM methods. These simulations have been widely applied to a variety of systems, including molecular clusters,^170^ solid materials,^171,172^ solutions,^173,174^ and biomolecules.^69^ Additionally, ML force fields are particularly well-suited for simulating reactive dynamic processes. To illustrate the potential of ML force fields in catalysis, this paper will focus on three primary application the stability of supported metals, catalytic reaction processes, and material structure prediction.

The surface energy of metal nanoparticles in the heterogeneous catalytic reaction is high, making sintering more likely, which subsequently leads to catalyst deactivation. A common solution is to disperse the nanoparticles on the surface of a substrate. The construction of PESs by neural networks offers comparable accuracy to QM but is several orders of magnitude faster, making it suitable for long-term simulations of nanoparticle sintering processes. To explore the influence of substrate on sintering behavior, Li et al. simulated the agglomeration process of gold nanoparticles on silica and ceria surfaces using neural network-based deep potential molecular dynamics (DPMD).^175^ The DPMD simulation results (Fig. 9a) showed that small gold nanoparticles are more likely to migrate on the surface of silica and rapidly merge with large nanoparticles. Similarly, on the flat CeO2−x (111) surface, small gold nanoparticles also migrate and agglomerate. In contrast, a gold nanoparticle at the step site of CeO2−x (111) is highly stable, with no significant migration during the simulation period. Further analysis revealed that the migration of the small particles is affected by the interaction between the metal and the support. Similar to hydrophilic behavior, gold nanoparticles on ceria exhibit smaller contact angles compared to silica, especially on the step site of the surface. These simulation results were also verified by experimental characterization. The rapid migration of gold nanoparticles on the surface of silica was observed using sequential high-resolution TEM (HRTEM), where all nanoparticles merged into a single particle. Conversely, Au nanoparticles on the surface of ceria were stabilized at the step site.

The precision and speed with which PESs are constructed using neural networks enable the simulation of more MD trajectories for the reaction process. For example, Jiang et al. investigated the equilibration dynamics of hot oxygen atoms following the dissociation of O2 on Pd(100) and Pd(111) surfaces using MD simulations based on a scalable neural network force field.^176^ By analyzing hundreds of trajectories, they found that, on both surfaces, oxygen atoms produced by oxygen dissociation tend to neighboring sites and perform a random-walk-type motion. This mechanism results in a finite distance distribution of equilibrium atoms, which is consistent with experimental observations. And the initial molecular orientation and surface thermal fluctuations have significant effects on the overall dissociation kinetics. Similarly, a study on the decomposition dynamics of formate (HCO2) on Cu surfaces demonstrated the application of the neural network force field in catalytic processes.^177^ Based on globally accurate high-dimensional PESs fitted with density functional theory data, the study predicted the mean translational energy distribution and angular distribution of desorbed CO2 on Cu(111) and Cu(100) surfaces. Additionally, the decomposition of HCOOH on different Cu surfaces is structurally sensitive due to different surface repulsions. These studies represent a significant advancement in modeling surface reactions using ML force fields, providing a more comprehensive understanding of the state of reactive species on catalyst surfaces and improving theory-experiment agreement.

ML force fields can also be used to predict the composition and structure of materials. For example, the thermodynamic phase diagram of metal oxide alloys remains largely unknown due to their compositive variability and associated atomic structural complexity, as does the catalytic kinetics on different surfaces of different compositions. By combining a global neural network potential with stochastic surface walking global optimization, Liu et al. investigated the relationship between the structure of ternary zinc-chromium oxide (ZnCrO) catalysts and the performance of syngas (CO/H2) conversion.^178^ They explored the PES structure under different components (Zn~xCryO) and constructed the ternary Zn–Cr–O phase diagram (Fig. 9d). Further calculation of the formation energy revealed the presence of a small, stable composition island in the phase diagram, where the oxide alloy crystallizes primarily in a spinel phase (Fig. 9e and f). Two representative crystal phases (ZnCr2O4~ and Zn3Cr3O8) were selected for further analysis of the syngas conversion mechanism. Reaction kinetics results indicated that a planar [CrO4] site, dynamically formed under reaction conditions, is the active center for methanol production. This planar [CrO4] is only present when Zn : Cr exceeds 1 : 2, after the appearance of [ZnO6], demonstrating that ZnCrO catalysis is highly sensitive to the Zn : Cr ratios. This work highlights the role of ML force fields in predicting the composition and structure of materials, providing valuable insights into the connection between atomic structure and its properties.

Overall, ML force fields have rapidly emerged as a transformative tool for heterogeneous catalysis, offering near-quantum accuracy at a lower computational cost. By learning from large and diverse QM datasets, ML force fields can capture complex chemical environments, including surface reconstructions, multi-element active sites, adsorbate-induced electronic effects, and reaction pathways, while enabling simulations at scales inaccessible to direct QM methods. Recent advances in equivariant graph neural networks (e.g., MACE, NequIP, GemNet) and databases (e.g., MP, OCP) have further improved data efficiency, accuracy, and transferability across different catalyst surfaces and reaction intermediates. Notably, MACE-OFF^179^ represents a significant advancement by extending the MACE architecture to a large-scale pretraining paradigm, analogous to foundation models in natural language processing. MACE-OFF is trained on millions of atomic environments spanning diverse elements, bonding motifs, and structural phases. By leveraging this pretraining, the model can be adapted to new catalytic systems with minimal fine-tuning or, in some cases, without any additional retraining. And MACE-OFF has demonstrated near-DFT accuracy across a range of downstream tasks, including molecular dynamics, adsorption energy predictions, and surface reaction energetics, with orders-of-magnitude lower data requirements compared to conventional task-specific ML force fields. Moreover, SO3LR^180^ further pushes the frontier by combining SO(3)-equivariance with explicit long-range interaction modeling. This architecture allows the model to capture electrostatics, polarization, and other nonlocal effects that are critical in catalytic environments but often neglected in standard ML force fields. Unlike most equivariant graph neural networks, which focus primarily on short-range many-body interactions, SO3LR incorporates efficient formulations of long-range physics directly into the network, enabling accurate treatment of charged surfaces, polar adsorbates, and extended catalytic interfaces. Importantly, SO3LR achieves these capabilities while remaining computationally scalable to large systems, positioning it as a promising candidate for modeling realistic catalytic reactors or electrochemical interfaces.

Despite their promise, ML force fields face significant challenges. Their predictive power is fundamentally limited by the quality, diversity, and representativeness of the training data, making out-of-distribution generalization a key bottleneck. Transferability across phases, surface terminations, and adsorbate coverages often requires retraining or fine-tuning, and constructing high-quality datasets for reactive, charged, or open-shell catalytic systems can be computationally prohibitive. Moreover, most current ML force fields lack explicit treatment of long-range electrostatics, charge transfer, and excited-state effects, which are critical in electrocatalysis, photocatalysis, and plasmon-assisted catalysis. Benchmarking across relevant catalyst classes and establishing rigorous uncertainty quantification protocols remain underdeveloped.

Looking forward, the integration of active learning, transfer learning, and foundation models (e.g., those developed in the OCP) offers a promising route to enhance data efficiency and generalization. Hybrid simulations that couple ML force fields with reactive force fields or on-the-fly quantum refinement can provide a balanced description of reactivity and long-timescale dynamics. Furthermore, expanding benchmark datasets to include realistic catalytic conditions such as solvent effects, applied potentials, and high coverage regimes will be essential to closing the gap between ML force field predictions and experimental observables, thereby enabling their widespread adoption in predictive catalyst design.

Summary and outlook

In summary, this review provides an introduction to the principles of force field-based PESs, highlighting three main types of force fields, their corresponding fitting methods and application. The evolution of force field-based PESs reflects a progression from simplified models to more precise, comprehensive, and adaptable approaches. Furthermore, this review offers a detailed comparison of the strengths and limitations of each type, alongside a discussion of their applicability to different systems, with a particular focus on heterogeneous catalysis. These advancements encompass studies on supported metal nanoparticle morphology, interfacial molecule distribution and diffusion, catalytic reaction processes, and the prediction of catalytic materials, though they do not represent an exhaustive list of all possible examples. By leveraging the cost-effectiveness of force field methods, researchers have overcome the limitations of QM, enabling simulations of previously inaccessible phenomena and significantly advancing scientific understanding.

Although force field methods have seen widespread application, they still face several challenges that lead to the inconsistency between the experimental results and the computational simulations. The key to solving this issue lies in improving the precision of constructing PESs using force fields. The first challenge is in database construction. Since experimental data are often costly to obtain, force field fitting typically relies on more readily available computational data. However, computational data often involves too much simplification, which can compromise the accuracy of the PES. To address this, a hybrid dataset can be built, where computational data are used to fit the force field, and the resulting model is validated and fine-tuned using experimental data, thereby reducing the gap between simulation and experiment. The second challenge concerns the architecture of the force field model. QM methods achieve accuracy through electronic-level descriptions, while force field methods still have shortcomings in electrostatic interactions. For example, the QEq method, which is primarily used for static charge distribution estimation, cannot deal well with the redistribution and transfer of electrons in non-equilibrium systems such as chemical reactions. Therefore, future force field models should aim to improve charge equilibration methods, so as to explicitly consider electrostatic interactions and non-local effects. Furthermore, the tradeoff between the accuracy and computational cost in force field fitting needs to be considered. Force field fitting typically requires significant computational resources and data, especially for high-precision fitting. Traditional fitting methods may require thousands of simulations to adjust and optimize parameters, making the process time-consuming. As a solution, the concept of meta-learning could be introduced to fit force fields for new systems. By selecting optimal algorithm and parameter configurations based on knowledge acquired from historical fitting tasks, the meta-learning model can provide an effective fitting effect and strong generalization capabilities, even in the case of sparse data.

Beyond the general limitations discussed above, several catalysis-specific challenges remain urgent. First, the transferability of force fields, particularly ML force fields, across different phases, chemical compositions, and surface states is often limited by the diversity of their training data. Addressing this requires multi-fidelity and transfer learning strategies that leverage both large general-purpose datasets (e.g., MP, OCP, MD17) and system-specific fine-tuning. Second, long-range electrostatics and charge transfer remain difficult to capture accurately in heterogeneous environments, especially under dynamic or reactive conditions. Extensions such as polarizable models and explicit charge equilibration schemes within ML frameworks show promise. Third, electrochemical environments, involving applied bias or constant potential conditions, introduce additional complexity for reactive force fields and ML force fields, necessitating integration with constant potential molecular dynamics and implicit/explicit solvation models. Fourth, fitting procedures and data generation for high-dimensional potentials remain computationally expensive, and active learning, data selection, and meta-learning approaches offer routes to reduce cost without sacrificing accuracy. Finally, benchmark datasets such as the OCP, MD17, and MP are critical not only for training but also for enabling reproducible model comparisons, helping guide community progress toward more robust and transferable force fields.

Author contributions

Chenglong Qiu, Tore Brinck and Jiacheng Wang made substantial contributions to discussions of the content. Chenglong Qiu prepared the draft manuscript. Tore Brinck and Jiacheng Wang revised the manuscript before submission.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.